Rehabilitated Tailing Piles in the Metropolitan Ruhr Area (Germany) Identified as Green Cooling Islands and Explained by K-Mean Cluster and Random Forest Regression Analyses

Abstract

Urban green spaces, such as parks, cemeteries, and allotment gardens provide important cooling functions for mitigating the urban heat island (UHI) effect. In the densely populated Ruhr Area (Germany), rehabilitated tailing piles (TPs), as relicts of the coal-mining history, are widespread hill-shaped landscape forms mainly used for local recreation. Their potential role as cooling islands has never been analyzed systematically. Therefore, this study aimed at investigating the TP surface cooling potential compared to other urban green spaces (UGSs). We analyzed the factors controlling the piles’ summer land surface temperature (LST) patterns using k-mean clustering and random forest regression modeling. Generally, mean LST values of the TPs were comparable to those of other UGSs in the region. Indices describing vegetation moisture (NDMI), vitality (NDVI), and height (VH) were found to control the LST pattern of the piles during summer. The index for soil moisture (TVDI) was directly related to VH, with the highest values on the north and northeast-facing slopes and lowest on slopes with south and southeast expositions. Terrain attributes such as altitude, slope, aspect, and curvature were of minor relevance in that context, except on TPs exceeding heights of 125 m. In conclusion, we advise urban planners to maintain and improve the benefit of tailing piles as green cooling islands for UHI mitigation. As one measure, the soil’s water-holding capacity could be increased through thicker soil covers or soil additives during mine tailing rehabilitation, especially on the piles’ south and southeast expositions.

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 improve 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) show more substantial cooling effects than smaller green patches [15,16]. Kuang et al. [17] even found an LST difference of over 6 °C between UGSs 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, and 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 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, and crown diameter and their effect on cooling [14,28]. For example, Rahman et al. [29] observed that fast-growing tree species such as Prunus umineko and P. calleryana generally show a higher stomatal conductivity, providing more cooling. Stumpe et al. [14] observed that tree age controls 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. [33] observed a significant negative correlation between tree canopy height and LST (r = −0.83).

Besides configuration and composition, other factors, such as water availability and terrain attributes, have been discussed in the context of the LST patterns. Concerning water availability, a strong negative correlation has often been found between the 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,34,35,36]. In this context, the Temperature Vegetation Dryness Index (TVDI) became a well-established parameter for monitoring surface soil moisture [37,38,39,40].

Concerning terrain attributes, several studies have 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. [41] 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. [42] 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.

Most of these investigations focused on UGSs such as urban parks, urban forests, and urban gardens, with parks being most relevant as urban cooling islands. Consequently, the term “park cooling island (PCI)” arose [18,20,43]. More recent, studies also focused on other UGSs, such as peri-urban greens and cemeteries [14,44]. However, in the Ruhr Metropolitan Area (Germany), one of the biggest urban agglomerations in Germany, numerous rehabilitated tailing piles from the former coal-mining industry characterize the landscape. The tailing piles are man-made hills of 50 to 200 m altitudes, consisting of rock debris from deep coal mining. Since they are often located in residential areas they were systematically revegetated for local recreation several decades ago. Unlike typical park cooling islands, piles form sloped green areas with high relief energies, which can, thus, provide additional cooling effects through cool air drainage 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 aimed to systematically characterize and understand 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 as follows:

(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 Germany’s largest urban economic zone with about 5.5 million inhabitants and a total area of about 4400 km² [46,47] (Figure 1).

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, and slate rocks [47]. Most of the mining waste material was dumped in tailing piles, using only smaller amounts as underground stowing materials. Today, the tailing piles are mostly found along a broad east–west band stretching through the central Ruhr Area (Figure 1), corresponding to the distribution of the former coal mines. In the early 1980’s, the regional mining company (Ruhrkohle AG, Essen, Germany) 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).

2.2. UGS Delineation

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), and cemeteries (C) from their surroundings, we used the “Authorative Topographic-Cartographic Information System” (ATKIS), as presented in

Table 1. Characteristics and sources of the data used in this study.

Data	Source	Spatial Resolution	Reference
Vector data delineating parks, cemeteries, and allotments	Authorative Topographic-Cartographic Information System (ATKIS)	not specified	[50,53]
Vector data delineating the tailing piles	Ruhr Regional Association (RVR)	not specified	[52]
Azimuth and altitude of sun light	SunCalc	not specified	[54]
Spectral Bands 4, 5, 6, 10 of Landsat 8 OLI	Earth Explorer	30 m × 30 m	[55]
Digital elevation model	Geobasis NRW	1 m × 1 m	[53]
Digital terrain model	Geobasis NRW	1 m × 1 m	[53]
Meteorological data	German Meteorological Service (DWD)	not specified	[56]

. It is the official German nationwide open source digital database for topographic spatial data and is updated annually [50]. However, within the ATKIS database, tailing piles (TPs) are not assigned as one specific ATKIS object group. Therefore, we used data from the digital web map of the Ruhr Regional Association [51], which was provided by Heiko Geyer [52] for pile delineation (Figure 1).

However, just like the ATKIS data, the tailing pile map is based on the German digital basic landscape model (Digitales Basis-Landschaftsmodell) [53]. Consequently, each of the P, A, C, and TP polygons is characterized with an accuracy of ±3 m.

2.3. Data Collection and Processing

Table 1 summarizes the data characteristics and sources that are further explained in detail in the following section. Our overall methodological approach is documented in the flowchart of Figure 2.

2.3.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 days with high irradiation from April to September (termed summer scenes), available from 2013 to 2023 [55]. The winter months with a low sun level (October to March from 2013 to 2023) were also selected separately. Within these selected time slots, all relevant Landsat 8 images were filtered for 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 also characterized with respect to the climate conditions during satellite overpassing in

Table A1. Characterization of the different selected Landsat scenes according to their cloud coverage and climatic condition before and while satellite overpassing. All Landsat 8 scenes are from row 197, path 24 with C2L1, Tier1 data quality.

Date	Season	Cloud Cover of the Total Satelite Scene	Precipitation 1 3 Weeks Before Satelite Overpassing	Air Temperatur 1 During Satellite Overpassing	Potential Evaporation 2 During Satellite Overpassing
		%	mm	°C	mm h−1
21 July 2013 (10:29)	summer scene	0.03	4.82 ± 2.52	28.0 ± 1.02	7.00 ± 0.01
6 June 2014 (10:27)	summer scene	2.39	48.94 ± 3.31	20.8 ± 0.53	5.10 ± 0.22
30 August 2016 (10:27)	summer scene	1.55	20.46 ± 6.42	21.3 ± 0.50	4.47 ± 0.33
14 June 2017 (10:27)	summer scene	0.74	28.88 ± 3.33	21.3 ± 0.71	5.33 ± 0.33
9 July 2023 (10:27)	summer scene	4.58	93.96 ± 21.92	27.6 ± 0.71	5.67 ± 0.90
22 February 2022 (10:27)	winter scene	0.13	96.00 ± 11.66	7.0 ± 0.40	2.50 ± 0.50

¹ Mean values of five meteorological stations across the Ruhr Area: Bochum (555), Dinslaken (989), Essen (1303), Gelsenkirchen (1595), Duisburg (13,670); ² Mean values of three meteorological stations across the Ruhr Area: Bochum (555), Essen (1303), Duisburg (13,670).

(Appendix A).

The LST was assessed of all selected Landsat 8 OLI/TIRS images, separately, based on spectral band 10 (C2L1 dataset, TIRS 1, 10.6–11.19 µm) using a common mono-window algorithm (MWA) according to previous studies [57,58,59], ensuring correct relative LST values that are most important in our research. Nevertheless, for retrieving the correct absolute land surface temperature, other algorithms including split-window algorithms are still in discussion [60,61]. However, the TIRs band was resampled from the original 100 m to 30 m resolution in the delivered data product [55].

In the first step, the digital number (DN) values of band 10 were converted into the top of atmospheric (TOA) spectral radiance (Lλ) as follows:

Lλ = M_L × DN + A_L

(1)

where M_L is the band-specific multiplicative rescaling factor, DN is the digital number of the band 10 image, and A_L is the band-specific additive rescaling factor. Values for M_L and A_L were taken from the corresponding metadata file of the satellite image.

In the second step, the spectral radiance Lλ was converted into the brightness temperature (BT):

$$\mathrm{B}\mathrm{T}={\displaystyle \frac{\mathrm{K}2}{\mathrm{l}\mathrm{n}[\left(\mathrm{K}1/\mathrm{L}\mathsf{\lambda }\right)+1]}}$$

(2)

where K1 and K2 are band-specific thermal conversion constants from the metadata file.

In the third step, the normalized differentiated vegetation index (NDVI, Equation (8)) was used for emissivity correction by calculating the proportion of vegetation P_v (Equation (3)) and emissivity (Equations (4) and (5)):

$${\mathrm{P}}_{\mathrm{v}}={\left({\displaystyle \frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{s}}}{{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{v}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{s}}}}\right)}^2$$

(3)

where NDVI_S and NDVI_V are global NDVI thresholds for soil and vegetation of 0.2 and 0.5, respectively.

Ɛ = Ɛ_v × P_v + Ɛ_s (1 − Pv) + C

(4)

where Ɛ_v and Ɛ_s are the vegetation’s and soil’s emissivity’s, respectively. C is the surface roughness which is a constant value of 0.005. The following equations can represent the conditions:

$${\mathsf{\epsilon }}_{\mathsf{\lambda }}=\left\{\begin{array}{c}{\mathsf{\epsilon }}_{\mathrm{s}\mathsf{\lambda },} \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}<{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{s}}\\ {\mathsf{\epsilon }}_{\mathrm{v}\mathsf{\lambda }}\ast {\mathrm{P}}_{\mathrm{v}}+{\mathsf{\epsilon }}_{\mathrm{s}\mathsf{\lambda }}\left(1-{\mathrm{P}}_{\mathrm{v}}\right)+\mathrm{C}, \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{s}}\le \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\le \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{v}}\\ {\mathsf{\epsilon }}_{\mathrm{s}\mathsf{\lambda }}+\mathrm{C} \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}>\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{v}}\end{array}\right.$$

(5)

Pixels with NDVI < 0 are classified as water with an emissivity value of 0.991. For 0 < NDVI < 0.2, bare soil is assumed with an emissivity value of 0.996. For 0.2 < NDVI < 0.5, it is assumed that a mixture of bare soil and vegetation covers the land. Pixels with NDVI > 0.5 are considered fully covered with vegetation with an emissivity value of 0.973.

In the last step, the emissivity corrected LST (in °C) was assessed as follows:

$$\mathrm{L}\mathrm{S}\mathrm{T}={\displaystyle \frac{\mathrm{B}\mathrm{T}}{\left\{1+\left[\left(\mathsf{\lambda }\mathrm{B}\mathrm{T}/\mathsf{\rho }\right)\mathrm{l}\mathrm{n}{\mathrm{Ɛ}}_{\mathsf{\lambda }}\right]\right\}}}-273.15 °\mathrm{C}$$

(6)

where LST is the land surface temperature in °C, λ is the mean wavelength of the band 10 emitted radiance of 10.895, Ɛ_λ is the emissivity calculated in Equation (5), and ρ is defined by

$$\mathsf{\rho }=\mathrm{h}{\displaystyle \frac{\mathrm{c}}{\mathsf{\sigma }}}=1.438\times {10}^{-2}\mathrm{m}\mathrm{K}$$

(7)

where $\mathsf{\sigma }$ is the Boltzmann constant (1.38 × 10⁻²³ JK⁻¹), h is the Planck’s constant (6.626 × 10⁻³⁴ Js), and c is the velocity of light (2.998 × 10⁸ ms⁻¹). Finally, the Kelvin value of the LST is transformed into Celsius by adding the absolute zero (approx. −273.15 °C).

All calculations were performed using the programming language R-4.2.2 [62], including the packages raster and sp.

2.3.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, NIR) and the red Landsat 8 OLI band 4 (0.64–0.67 μm, R) according to [62]. The calculation of the index is documented in Equation (8):

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}}$$

(8)

NDVI values range from −1 to 1, with positive values representing vegetated areas and negative values representing non-vegetated areas [63].

The NDMI index is used as a proxy for vegetation moisture [64], which is a highly important factor for the TP characterization due to the different shading regimes of the hill-shaded piles. This index was calculated from near-infrared 8 OLI band 5 (0.85–0.88 μm, NIR) and short-wave infrared Landsat 8 OLI band 6 (1.57–1.65 μm, SWIR), as documented in Equation (9).

$$\mathrm{N}\mathrm{D}\mathrm{M}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}}$$

(9)

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 [65].

Vegetation height (VH) was assessed based on a normalized digital surface model (nDSM) according to Yu et al. [66]. Since the nDSM is a derivate elevation product obtained by subtracting a digital terrain model (DTM) from a DSM, it can represent the relative height of features above the surrounding ground surface [67]. Since all UGS sites are vegetation-dominated, the nDSM are well suited to assess vegetation heights. The corresponding DSM was derived by airborne laser scanning (ALS) with horizontal and vertical accuracies of 1 m × 1 m and 0.1 m, respectively [53].

2.3.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 short-wave infrared Landsat 8 OLI band 6 (1.57–1.65 μm, SWIR) and thermal infrared Landsat 8 OLI band 10 (10.60–11.19 μm, TIR), as it is documented in Equation (10):

$$\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{a}\mathrm{I}={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}-\mathrm{T}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}+\mathrm{T}\mathrm{I}\mathrm{R}}}$$

(10)

NDBaI ranges between −1 and 1, with values above zero indicating an increasing proportion of bare land and values below zero typical for vegetated areas [67].

The TVDI (Temperature Vegetation Dryness Index) can assess soil moisture. It was determined based on the method proposed by Sandholt et al. [37]. 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 [37]. 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 is defined by the following equation (Equation (11)):

$$\mathrm{T}\mathrm{V}\mathrm{D}\mathrm{I}={\displaystyle \frac{\mathrm{L}\mathrm{S}\mathrm{T}-\mathrm{L}\mathrm{S}\mathrm{T}\mathrm{m}\mathrm{i}\mathrm{n}}{\mathrm{L}\mathrm{S}\mathrm{T}\mathrm{m}\mathrm{a}\mathrm{x}+\mathrm{L}\mathrm{S}\mathrm{T}\mathrm{m}\mathrm{i}\mathrm{n}}}$$

(11)

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 calculations, 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 equations:

LSTmin = (a + b × NDVI)

(12)

LSTmax = (c + d × NDVI)

(13)

Thus, for each pixel’s observed NDVI values, an LSTmin and LSTmax value was calculated and used for the TVDI calculation (Equation (11)). 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 calculations, regression outliers were eliminated using Cook’s Distance which summarizes how much a regression model changes when removing the ith observation. Data points from regression processing were excluded if Cooks’s Distance was higher than 4/n. After outlier correction, we obtained validated regression models with R² generally higher than 0.96.

2.3.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 a 1 m × 1 m horizontal resolution and a vertical accuracy of 0.1 m from campaigns in 2019 and 2020 [51].

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 (Version 3.28.4-Firenze). Also, hill shading attributes were calculated for each date of the different Landsat scenes, according to the corresponding azimuth and altitude of light [54], since the hill shading effect depends on the position of the sun, which changes over the year and day.

2.3.5. Data Downscaling

For all datasets, we generated a reference grid with a 30 m × 30 m horizontal resolution, which accurately overlaps the grid of the Landsat 8 bands (1 to 7). This grid was then used as a reference basis. For all spectral indices (NDVI, NDMI, NDBaI, TVDI) and for the retrieved LST (based on resampling during Level-1 product processing), one value per grid cell was given.

Due to the higher horizontal resolution of 1 m × 1 m of all the other vegetation and terrain variables used, their arithmetic means were calculated for each grid cell to accurately relate all variables to each other in further regression analyses.

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 (Equation (14)):

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

(14)

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 generally vary between 0 and $\pi$ as well as between 0 and −$\pi ,$ representing, in our case, east and west expositions of the tailing piles, respectively. However, since we assumed that, according to the azimuth of light, the north and south expositions of the tailing piles were more important for their LST pattern, then the absolute values (ABSs) 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 expositions of the tailing piles, respectively. However, we could no longer distinguish between the east and west expositions.

2.4. 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 described below.

2.4.1. Cluster Analyses

In order to differentiate between tailing pile LST-types, a k-mean classification was performed based on the five mean LST values of each summer date of each tailing pile. After clustering, we further characterized each cluster regarding pile vegetation status and topographical parameters.

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 errors in the Euclidean distance over all data points [68,69,70,71].

For the cluster analyses, we used the R (Version 4.2.2) package Factoextra, which automatically performs an initial principal component analysis (PCA) if there are more than two dimensions (variables). The corresponding plot shows the cluster data points according to the first two principal components explaining most of the variance. In our case, an initial PCA was calculated based on the LST values of the five different dates of the selected Landsat 8 scenes before clustering started.

The optimal number of clusters was identified using the Silhouette method as described in Kodinariya and Makwana [70]. Generally, this method compares within-cluster distances with between-cluster distances, where the greatest difference between clusters indicates the best number of clusters.

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. This means scores of more than 0.5 indicate high-quality cluster assignments.

For a final TP typification, the resulting clusters were characterized in terms of their mean total tailing pile area, spectral indices, and terrain attributes. Significant differences between these means of cluster characteristics were determined using ANOVA analyses. Since ANOVA requires assumptions of variance analyses such as normality and homogeneity, normality tests with Kolmogorov–Smirnov tests, Q-Q plots, histogram comparison, and Levene’s test were calculated according to Cilek et al. [72] before ANOVA testing was performed. Significant differences in cluster characteristics between the clusters were subjected to the Tukey test as a post hoc statistic with a significance level of 0.05.

2.4.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 m × 30 m). Due to the small sample size of 82 piles, we refrained from using the mean pile LST as a dependent variable. We also tested multiple linear regression (MLR), but they failed in prediction accuracy, with R² being lower than those of the random forest models. The RFR analyses generally belong to the more complex supervised learning algorithms [73,74] 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,74].

However, in our case, the LST values of all grid cells covering the tailing pile expansion were used as the dependent variable for RFR modeling, while separate RFR models were run to explain the LST distribution for the different dates of the selected Landsat summer scenes.

For the RFR models, various vegetation (NDVI, NDMI, and VH) and TP terrain characteristics (altitude, slope, curvature, aspect, and hill shading values) were used as independent and explaining variables.

For RFR modeling, as a first step, training and testing data for the different RF models were set to 70 and 30 percent, respectively, which is in accordance with numerous other studies [75,76,77]. To gain valid RF models, a random subset of training parameters was supervised iteratively in relation to the number of regression trees to reduce the error and attain greater accuracy without overfitting according to Belgiu and Dragut [76]. For each RFR model, the number of trees was optimized to obtain the most stable performance. The results showed that across all RFR models, 300 trees were sufficient for stable RFR performance since no further enhancement in model performance was attained beyond this. Based on this, the regression models were finally fitted using the training datasets in the third step, and the dependent variable (LST) was predicted using the trained models and the corresponding test datasets.

The predicted and the observed variables were used to calculate error metrics, such as the coefficient of determination R², the root mean squared error (RMSE), and the mean absolute error (MAE). These metrics were calculated separately for each date-dependent RFR model and were finally averaged for a better overview of the RFR performances.

To further identify the importance of the independent variables, variables are removed from RFR analyses, assessing how the prediction changes [78]. 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 removed randomly. A higher %IncMSE value represents higher variable importance [79]. 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 [80]. Thus, both indicators evaluate the variable’s importance in the RFR models, without distinguishing between positive or negative effects on the LST variance. Analysis was conducted in R (version 4.2.2) [62] using the Metrics package. The %IncMSE and the IncNodePurity were calculated separately for each date-dependent RFR model and were finally averaged for a better overview of the important variables. However, to describe the variability of the %IncMSE and of the IncNodePurity, their standard deviations across the different RFR models were also calculated.

RFR modeling was performed in R 4.2.2 by using the package “randomForest”.

3. Results

3.1. Tailing Pile Characteristics

As shown in Figure 1, most tailing piles are located within or near settlements as indicated by high Normalized Build-up Index (NDBI) values and, thus, are considered as urban greens. While some of the larger tailing piles are adjacent to periurban areas with lower NDBI values, the smaller ones are often embedded in 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, and cemeteries (

Table 2. Characterization of the four different classes of Urban Green Spaces (UGSs) in the Ruhr Area in terms of spatial expansion (% and total area, mean, min, max) and land surface temperatures (LST) with mean and standard deviations determined from thermal satellite images at five different summer dates.

			Tailing Piles (n = 82)	Parks (n = 2233)	Allotments (n = 1409)	Cemeteries (n = 640)
Spatial expansion	total area	%	0.51	1.46	0.79	0.69
	total area	km2	25.3	64.70	35.4	30.6
	Ø area	km2	0.323	0.029	0.025	0.048
	min area	km2	0.030	0.0001	0.002	0.0001
	max area	km2	22.61	2.14	0.19	1.04
	21 July 2013	°C	30.83 ± 2.96	29.79 ± 2.07	30.26 ± 1.47	28.98 ± 1.69
	6 June 2014		22.51 ± 1.92	23.38 ± 1.89	23.15 ± 1.14	22.77 ± 1.59
	30 August 2016		23.55 ± 2.78	23.15 ± 1.74	24.15 ± 0.05	22.38 ± 1.55
	14 June 2017		25.99 ± 3.67	25.73 ± 2.52	26.62 ± 1.85	24.88 ± 2.11
	9 July 2023		22.98 ± 2.23	23.02 ± 1.71	23.61 ± 1.33	22.53 ± 1.47

). Nevertheless, tailing piles comprise large areas with a mean size of 0.323 km², 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 2).

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

Thermal Footprints

Table 2 summarizes the mean LST values of all tailing piles and the other UGSs for the five different analyzed 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 is generally distinctly higher for the piles than for the other UGSs, indicating a higher heterogeneity within this UGS class.

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 4A). 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. When analyzing the LST pattern in relation to slope exposition, south- and south-east-facing slopes were the warmest and the northern to western slopes were the coolest (Figure 4B). However, the LST differences between the slope expositions are not significant, amounting to only 1.2 °C and, thus, only partly explaining the large LST variability shown in Figure 4A.

3.2. Thermal Typification of Tailing Piles

The Silhouette method of the k-mean classification showed that the optimal number of clusters was four (Figure 5A). Figure 5B and

Table 3. Cluster separation quality and characterization of the different thermal clusters according to terrain, vegetation, and soil attributes.

Mean Properties		Cluster	Cluster	Cluster	Cluster
		1	2	3	4
Cluster Size		7	15	40	19
Mean Cluster Silhoutte Scores 1		0.28	0.29	0.33	0.40
LST *	°C	28.17 ± 2.84	26.05 ± 1.43	24.31 ± 1.12	22.88 ± 0.71
Area	km2	0.64 a	0.40 a	0.27 a	0.19 a
Altitute	m	89.72 a	88.54 a	96.04 a	91.16 a
Slope	°	11.89 a	13.48 a	15.45 a	14.03 a
Curvature	°	0.02 a	0.05 a	0.02 a	0.03 a
NDVI *	-	0.229 a	0.340 a	0.385 b	0.411 c
NDMI *	-	0.070 a	0.147 b	0.191 c	0.221 d
VH	m	2.07 a	3.61 b	7.03 c	10.04 bc
NDBaI	-	−0.556 a	−0.536 c	−0.546 bc	−0.558 b
TVDI *	-	0.517 b	0.482 c	0.374 a	0.270 a

^abcd indicated significant difference of ANOVA analysis (TUKEY posthoc with a significance level of 0.05); ¹ based on an average silhoutte of 0.68 had been calculated; * cluster mean values of the summer sampling dates: 21 July 2013, 6 June 2014, 30 August 2016, 14 June 2017, 9 July 2023.

show the cluster assignments and their quality, respectively. Moreover, Figure 5B shows that the first two components of the PCA explain up to 90.5% of the data variance of the five LST summer Landsat scenes, guaranteeing valid cluster results that were mapped along these two dimensions. In this context, Cluster A comprises the tailing piles with the highest LST values, with a mean LST of 28.17 °C (Table 3) followed by Clusters B, C, and D with a descending order of LST means of 26.05, 24.31, and 22.88 °C, respectively. Although the clusters show 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 shows 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 allowed for further cluster characterization.

Mean terrain attributes such as the mean tailing pile altitude (DEM) and the mean slope and curvature showed no significant differentiation between the four clusters (ANOVA analyses, p-values > 0.05) (Table 3). The mean tailing pile area seemed to be related to the cluster LST characteristics since the mean area sizes declined with the mean cluster temperatures, but the mean areas 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 (Table 3). Soil moisture (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 negative values ranging from −0.558 to −0.536, the mean soil bareness index (NDBaI) indicated a high degree of vegetative cover and was quite similar among the four clusters, although differences were partly significant. This shows that vegetation parameters and soil moisture are closely related to the differentiation between the four pile LST types.

For a visualization of these results, Figure 6 shows one arbitrarily selected tailing pile from each cluster as aerial images (A1–D1), DEMs (A2–D2), and LST patterns (A3–D3). As the aerial images and the DEMs indicate, the tailing piles generally form hills higher than the surroundings 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 6(A-1)), belonging to the hottest cluster, is characterized by a high proportion of bare soil (Figure 6(A-3)), while the cool piles “Tetraeder” and “Lohberg Nord” (Figure 6C,D) are dominated by tree vegetation.

Figure 6. Representative tailing piles (Kohlenhuck, Hohewardt, Tetraeder, Lohberg Nord) for the four different thermal clusters (A–D) shown as orthophotos (first column), altitude (second column), and land surface temperature (third column).

Figure 6. Representative tailing piles (Kohlenhuck, Hohewardt, Tetraeder, Lohberg Nord) for the four different thermal clusters (A–D) shown as orthophotos (first column), altitude (second column), and land surface temperature (third column).

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 ran RFR models with the LST as dependent and various independent variables shown in Figure 7, 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 4. Model performance of the different random forest regression models for explaining the variability of mean LST values at the pixel level, separated for summer and winter months.

RFR Modells	Period	Model Performance
		R2	RMSE	MAE
RFR 1 (n = 32,007)	Summer 1
all TP data		0.85 ± 0.04	0.39 ± 0.04	0.29 ± 0.03
RFR 2 (n = 24,032)
TPs with DGMmax < 125 m		0.86 ± 0.03	0.38 ± 0.03	0.29 ± 0.03
RFR 3 (n = 7770)
TPs with DGMmax ≥ 125 m		0.83 ± 0.03	0.42 ± 0.03	0.31 ± 0.02
RFR 1 (n = 32,007)	Winter 2
all TP data		0.75	0.49	0.38
RFR 2 (n = 24,032)
TPs with DGMmax < 125 m		0.74	0.49	0.38
RFR 3 (n = 7770)
TPs with DGMmax ≥ 125 m		0.83	0.40	0.31

¹ mean values of the summer sampling dates: 21 July 2013, 6 June 2014, 30 August 2016, 14 June 2017, 9 July 2023; ² values of the winter sampling date: 28 February 2022.

). Generally, we obtained high RFR model performances with a mean R² of 0.85 and an RMSE of 0.39 (Table 4). To evaluate the importance of the independent variables across all satellite summer scenes, the %IncNodeMSE values are visualized as mean values in Figure 7A. In accordance with the cluster characterization, the LST variance of the piles was found to be mainly controlled by vegetation properties, which is 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, and 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 7B), resulting in a less accurate RFR model with an R² of only 0.75 (Table 4). However, in winter, with a low sun altitude, hillshade and aspect dominate the LST distribution of the piles, especially the aspect. Due to the dominance of annual herbs, deciduous trees, and shrubs, vegetation is inactive in winter and, thus, is of minor importance for LST variability during this season.

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

Figure 8. 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 in the aspect importance in the RFR models occurs.

Figure 8. 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 in 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 elevations above and below 125 m ASL. Both datasets were then subjected to the RFR analyses, with the results for both analyses shown in Table 4 and Figure 7. Concerning the Landsat summer scenes, the RFR model performance showed no considerable differences between the three datasets (all pile pixel data, pixel data of piles with DGMmax < 125 m ASL and with DGMmax > 125 m ASL), but with respect to the Landsat winter scenes model performance improved using pile data with DGMmax > 125 m. With an R² of 0.83, the LST variance of the pixel LST values of piles with DGMmax > 125 m 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 < 125 m ASL and with DGMmax > 125 m ASL, separately. For the Landsat summer scenes, the importance of the aspect and the VH increased when pixel data of piles with DGMmax > 125 m ASL were used, whereas simultaneously, the importance of the NDVI, NDMI, and NDBal decreased. 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 dataset separation concerning the maximal pile height resulted in more precise and more differentiated RFR modeling results that highlight the importance of the pile aspect and slightly modify the importance of vegetation properties. Nevertheless, vegetation characteristics were 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, 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.

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

Across all piles, mean TVDI values varied between 0.30 and 0.46, where the lowest values (moist soils) were found on the north and northwest slopes of the pile hills (Figure 9). Complementary to this, the highest values and, thus, soils with the lowest moisture were found in the south and southeast expositions of the piles. Interestingly, the spatial distribution of the mean TVDI seemed 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 Figure 9. Here, the highest vegetation heights are found in the northwest and west while the lowest occur in 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 similar to those of the conventional UGS parks, allotment gardens, and cemeteries across the Ruhr Area. Among these three classes, cemeteries generally have the lowest LST values, which is mainly attributed to the higher density and older age of trees on these sites [14], while allotments are warmest due to a higher presence of sealed surfaces from sheds and pathways [14,81,82]. The mean LST of the tailing piles is most similar to that of the urban parks in the study area, most likely due to a similar structure with open spaces interspersed by trees and forest stands. While the role of the conventional UGSs as cooling islands in UHI mitigation has been demonstrated many times [7,8,9,10,11,12,13,14], our study, for the first time, shows that rehabilitated tailing piles can be of similar significance. This is a more or less a 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, which was also in the Ruhr Area. Compared to flat urban parks, such piles have an additional cooling benefit during the nighttime when cool air can drain down the slopes into adjacent residential areas. This highlights an aspect we could not consider but that can 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 was also studied by Mathew [83] and Tan et al. [34], who generally explained this by a terrain-dependent variation in solar incident radiation.

4.2. Factors Controlling LST on Tailing Piles

The k-means analysis based on the mean pile LST values showed that the four separated clusters were mainly characterized by differences in vegetation parameters while terrain attributes of the piles showed no statistically significant differences between the LST clusters. Clearly, the vitality and density of the vegetation, as reflected in the NDVI, and even more its moisture content (NDMI) and the height of trees (VH), showed a pronounced and significant differentiation between the four clusters, with the highest values in the coolest pile cluster, thus underlining the generally acknowledged importance of evapotranspirative cooling from a dense and vital vegetation [18,19,20,33].

Similarly, the random forest analysis showed that 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. The RFR models reached a very high prediction accuracy of 83 to 86% rarely found in other similar studies [22,27,28]. This is most likely related to the complex tailing pile environment showing a high heterogeneity of grass, shrub, and tree vegetation, so that the NDVI and NDMI indices vary greatly and, thus, have a high impact on LST patterns. The most important factor in the RFR models was the NDMI index, which reflects the importance of plant available water as the limiting factor for transpirative cooling processes, as also reported in several other studies [33,34,35,84].

The influence of vegetation height (VH) on LST has only been analyzed in a few studies [85,86,87,88] because it requires data from expensive aerial laser scan campaigns. In the study of Chen et al. [88], the species composition of the tree canopy showed the highest impact on LST variability, but its vertical structure was also an important impact factor. Gage and Cooper [87] 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. In addition, tree height and aboveground biomass are generally closely related to the expansion of the rooting system [89], so that larger trees also have more access to soil water and, thus, can transpire even during longer dry summer periods when grasses, shrubs, and smaller trees may already suffer from water stress.

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,90,91]. 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 [90]. This opposing effect likely masks the slope’s potential influence 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 impeded on north-exposed hillsides only for these higher tailing piles. For the winter image, the aspect and the hillshade parameter modeled by QGIS became the most important controlling variables for LST because vegetation was bare or dormant, but also because the solar zenith angle dropped to about 15°, so that shading effects on slopes that were not exposed to the sun became more prominent, thus becoming an important factor for LST variability.

Studies like that of Bai et al. [41] and 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. In contrast to the mountainous areas from these studies, tailing piles are low hills with only few steep slopes, so that elevation effects on LST can be expected to be negligible and shading is of little importance during satellite overpassing around noontime on these moderate slopes. Furthermore, the man-made plant cover was varied mainly for aesthetic reasons and is, therefore, not related to any terrain attributes.

However, although the aspect parameter generally showed a low impact on 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 a 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 [92,93,94]. Thus, these studies frequently used the Topographic Wetness Index (TWI) to describe and predict vegetation height. For example, Mohamedou et al. [93] used the TWI 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, but it appears to be a promising factor in predicting vegetation heights. Unlike the TWI, the TVDI does not consider lateral water movement in soils, but this may not be relevant on the tailing piles since they were constructed from fragmented bedrock debris, so that 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 similar to that of other urban green spaces such as parks, cemeteries, and 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 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 and biochar [95,96], could be spread out on the surfaces or mixed into the topsoil to improve 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 the 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 and buildings at the hillfoot [45], 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 nighttime urban heat in the surrounding neighborhoods.

4.4. Limitations and Open Questions

In our study, we derived the LST values from band 10 of the Landsat 8 OLI/TIRS images using a mono-window algorithm, similar to other studies [57,58,59]. However, there are other common methods to assess the LST, for example, the split-window algorithm (SWA), which are known to retrieve LST values with another accuracy [60,61]. Based on the high diversity of algorithms used for LST calculation, we cannot guarantee to have used the one with the most accurate absolute LST prediction. Nevertheless, for our study, we mainly required relative LST differences, which can correctly be determined with the data from band 10 [57].

Many parts of this analysis are based on mean values from days with warm to hot air temperatures, low cloud cover, and no rainfall during the previous 2 weeks. This ensures 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 and plant phenology, possibly also management, affect the cooling potential differently in UGSs. Working with mean values allows for a better classification of the pile type 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 most tailings were rehabilitated before 2013 [97], the tree heights determined in 2019/20 may not represent the height distributions on the scenes from 2013 to 2016, or 2023. However, these were the best data available for the region and, once again, the model performance based on the mean values from the satellite scenes indicates that the found relationships are solid and meaningful.

Finally, the LSTs of the tailing piles may originate not only 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 [97,98]. 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 are currently unavailable at the required spatial resolution. In any case, such “hot spots” can only have minor effects on LST values since these are very well predicted from soil and vegetation parameters in our RFR models.

5. Conclusions

We identified the rehabilitated mine tailings in the metropolitan Ruhr Area as UGSs with a cooling function similar to that of urban parks or allotments, thus substantially contributing to the mitigation of the urban heat island. The analysis of tailing pile parameters derived from satellite, topographic, and LiDAR data showed that LST patterns of the piles during summer were largely controlled by the vegetation indices for vitality (NDVI) and for moisture (NDMI), as well as by vegetation height. The index for soil moisture (TVDI) was closely related to vegetation height and both parameters were lower on S- and SE-facing slopes; this suggests that soil moisture limits tree growth and, thus, the cooling potential on slopes that are subject to higher irradiation. To improve tree growth and the cooling potential on these slopes, the soil’s water-holding capacity could be enhanced through establishing thicker soil covers during tailing pile rehabilitation or through soil additives such as biochar and mulches that could be incorporated or distributed beneath the existing vegetation.