The study area covers approximately 35 km² and is delineated by the Fenhe Reservoir hydrological station upstream and the Zhaishang station downstream, encompassing 173 sub-watersheds across Gujiao City and Loufan County in Shanxi Province (Fig. 1a). The area is located in the central Loess Plateau, a major coal resource zone in China, and features highly dissected terrain with densely incised gullies and steep slopes. Elevations range from 954 to 1,804 m, with most slopes exceeding 15°. The region experiences a temperate continental climate, with a mean annual temperature of 9.5°C and average annual precipitation of 450 mm, primarily occurring between June and September. Surface soils are mainly cinnamon and loessial types derived from loess parent material, exhibiting weak structural stability and high erosion susceptibility (Shi et al., 2004).

The Gujiao section lies within the Xishan coalfield, a typical resource-dependent industrial mining area with a long history of coal extraction and a high density of mining operations (Fig. 1b). Coal-based solid wastes (e.g., gangue and fly ash) have been openly stockpiled between gullies and slopes over prolonged periods, causing substantial landscape disturbance, altered hydrological processes, and increased ecological fragmentation. Some waste dumps exceed 20 ha in area and reach heights of 30–40 m, with accumulated leachate at their bases threatening downstream cropland, forestland, and water bodies with potential contamination (Wang et al., 2025). Processes such as spontaneous combustion, weathering dust, and heavy metal migration contribute to multi-media contamination of air, water, and soil. This area provides a representative setting for assessing watershed-scale ecological and environmental risks.

Four types of raw data were used in this study: remote sensing imagery, basic geographic information, socio-economic data, and digital vector maps (Table 1). Landsat OLI imagery was obtained from the United States Geological Survey (USGS) Earth Explorer platform (https://earthexplorer.usgs.gov/) to generate various remote sensing indices. Terrain data for land surface relief were extracted from the ALOS PALSAR Digital Elevation Model (DEM), provided by the Japan Aerospace Exploration Agency (JAXA). Meteorological data, including annual average temperature and precipitation, were acquired from the Taiyuan meteorological station through the China Meteorological Data Service Center. Hydrological data were collected from the Zhaishang and Fenhe Reservoir hydrological stations. Population data were obtained from the WorldPop (https://www.worldpop.org), while GDP and other socio-economic data were sourced from the Statistical Yearbook of Shanxi Province and the official statistical bulletins of Gujiao City. The boundaries of coal-based solid waste dumps were delineated manually through visual interpretation of historical Google Earth imagery and refined using field survey data. Water quality and soil heavy metal data were derived from publicly available environmental monitoring reports issued by the Shanxi Provincial Department of Ecology and Environment. In areas with missing values, data from adjacent regions reported in published literature were used for spatial interpolation. Vector data for administrative boundaries and river networks were acquired from the National Geographic Information Monitoring Cloud Platform (http://www.dsac.cn) and OpenStreetMap (https://www.openstreetmap.org).

All spatial datasets were projected to the WGS 1984 UTM Zone 49N coordinate reference system, rasterized, and resampled to a spatial resolution of 30 m × 30 m. This resolution served as the minimum evaluation unit for watershed-scale environmental risk modeling. To ensure consistency in spatial and attribute dimensions across heterogeneous datasets, data preprocessing including cleaning, format conversion, projection transformation, and normalization was carried out using ArcGIS, ENVI, and Python.

To analyze the spatiotemporal evolution of coal-based solid waste dumps in the study area, we employed a multi-temporal remote sensing approach integrating Landsat satellite imagery and historical Google Earth data from 2010 to 2025 at approximately three-year intervals. These datasets facilitated long-term visual interpretation and temporal comparison of dump boundaries. Considering the irregular geometry, heterogeneous surface texture, and spectral similarity of dumps to surrounding land covers, a visual interpretation method was adopted. Interpretation criteria included surface color, topography, proximity to mining infrastructure, and morphological features indicative of waste accumulation (Kupidura et al., 2023; Thiruchittampalam et al., 2023). We manually digitized the dump boundaries using ArcGIS 10.5 to ensure spatial accuracy and temporal consistency. To enhance delineation precision, we derived the normalized difference built-up index (NDBI) and bare soil index (BSI) from Landsat TM and OLI images, which improved the visual contrast of dump features (Arif et al., 2024; Salas et al., 2023). Finally, we overlaid dump polygons from different years to analyze expansion patterns, growth directions, and spatial clustering characteristics.

To assess the environmental impacts of coal-based solid waste dumps on watershed systems, this study employs the drivers, pressures, states, impacts, and responses (DPSIR) framework to develop an integrated environmental risk assessment indicator system. The DPSIR was developed based on the concept of pressure–state–response (Hu et al., 2019) and was further developed to support systematic analysis of the complex interactions between anthropogenic activities and environmental change (Tang et al., 2024; Kazuva et al., 2018; Yussif et al., 2023). This five-dimensional analytical structure enabled a comprehensive assessment of disturbance pathways and the feedback mechanisms by which human-induced pressures affect watershed environmental systems. Based on the geomorphological features and pollution evolution patterns of the typical gully-type watershed in the Xishan mining area, a three-tier assessment framework was developed, comprising a target layer, a criterion layer, and an indicator layer.

This study developed an 18-indicator environmental risk assessment system based on the DPSIR framework, covering driving forces, pressures, states, impacts, and responses (Fig. 2). Driving forces include precipitation, temperature, and terrain relief, reflecting the natural background influencing pollutant dispersion. Pressure indicators encompass population density, economic level, waste dump density, and proximity to water bodies, representing direct human disturbances. State indicators describe ecosystem conditions through landscape diversity, vegetation health, and coverage. Impact indicators capture ecological vitality, water quality, and soil contamination consequences. Response indicators assess social governance and ecological restoration via healthcare capacity and land reclamation area. Each indicator was classified as either benefit-oriented (+) or cost-oriented (−) according to its environmental risk attributes (Boori et al., 2022).

Within the DPSIR framework, natural background conditions such as annual precipitation, temperature, and surface relief are considered driving forces that shape the baseline ecological processes in the watershed. External pressures, including population agglomeration, economic development, and the high-density distribution of coal-based waste dumps, result in significant environmental stress. These pressures lead to changes in ecosystem conditions, such as reduced biological abundance, vegetation degradation, and diminished ecological vitality. Such alterations further trigger compound environmental risks, including soil heavy metal contamination and surface water deterioration. Ultimately, these impacts necessitate government and institutional responses through ecological restoration, land reclamation, and public resource reallocation. The indicator selection and classification were informed by a comprehensive literature review, previous studies, data from government and non-governmental organizations, and expert interviews within the study region, thereby providing a robust theoretical and methodological foundation for subsequent weight assignment and spatial risk stratification (Kaur et al., 2020; Sun et al., 2015). Table 2 illustrates the stepwise indicator analysis process for watershed environmental risk under the DPSIR model.

To integrate multiple indicators from diverse sources with different units and scales, it is essential to first standardize all variables onto a common scale. Subsequently, the relative weights of these indicators are determined by using a combination of multi-criteria decision-making methods, enabling a comprehensive evaluation of multidimensional ecological risks.

To identify the dominant environmental risk factors driving spatial heterogeneity in watershed-level eco-environmental risk, we applied two complementary methods, Grey Relational Analysis (GRA) and Manhattan distance analysis (MDA). These approaches were chosen for their robustness in handling multidimensional datasets and their suitability for evaluating complex environmental systems under data limitations and uncertainties. GRA was used to quantify the strength of association between each indicator and the overall ERI. As a non-parametric technique rooted in grey system theory, GRA is effective in identifying key variables in systems with incomplete information (Zhang et al., 2024). The analysis involved constructing a reference sequence from the ERI values and comparing it with sequences derived from each normalized indicator. The grey relational coefficient was calculated for each pair at every spatial unit, and the resulting grey relational grade represented the overall correlation strength. Indicators with higher relational grades were interpreted as having greater influence on spatial variation in environmental risk. In parallel, MDA evaluated the spatial dissimilarity between each risk indicator and the composite risk surface (Nur et al., 2024). Manhattan distance was calculated by summing the absolute differences between indicator values and the ERI across all grid cells. Indicators with lower total distance values showed spatial distributions more similar to the ERI, indicating greater explanatory power in shaping the composite risk pattern. By integrating the results of GRA and Manhattan distance analysis, indicators that consistently ranked among the top in both methods were identified as the dominant contributors to watershed environmental risk. This combined approach enhanced the robustness of factor identification and provided a solid foundation for policy prioritization and targeted intervention planning.

Based on Google imagery from 2010, 2013, 2016, and 2019, together with Sentinel imagery from 2022 and 2025, visual interpretation was used to delineate the spatial distribution of coal gangue, fly ash, and other coal-based solid waste dumps and their associated reclamation areas in the study area (Fig. 4a). The dump area, reclamation area, and total impacted area for each period are summarized in Table 5. From 2010 to 2025, the dump area decreased from 3.7936 km² to 1.9038 km², while the reclamation area increased steadily from 0.4303 km² to 2.3457 km². The total impacted area reached a maximum of 4.7186 km² in 2013, dropped to a minimum of 3.7084 km² in 2019, and then rose slightly to 4.2495 km² in 2025 (Fig. 5).

Temporally, the dump area remained at a relatively high level from 2010 to 2013, decreased between 2013 and 2019 by approximately 18.1%, and then stabilized between 2019 and 2025. Over the same periods, reclamation area increased by approximately 55.4% from 2013 to 2019 and continued to grow thereafter. Spatially, dumps exhibited a discontinuous distribution pattern, occurring as discrete clusters in valleys and on gentle slopes near coal preparation plants, coal chemical facilities, and pithead power plants. Reclamation areas were generally located along gentle slopes and ridge belts surrounding the dumps, with some overlapping runoff convergence zones (Fig. 4b).

The relative stability of coal-based solid waste dump areas during 2019–2025 prompted the use of 2020 data, along with other indicators from the same year, to construct the EERI for environmental assessment (Fig. 6). The spatial distribution of the EERI showed marked heterogeneity, with high- and low-risk zones distributed unevenly across the watershed.

High-risk areas were primarily concentrated in sub-watersheds 150, 160, 25, and 26. Sub-watersheds 150 and 160 contained a high density of coal-based solid waste dumps, were located near major drainage channels, and experienced relatively high rainfall–runoff volumes. Sub-watersheds 25 and 26 encompassed the urban area of Gujiao City, characterized by dense building coverage and high population density. These locations exhibited significant spatial overlap of multiple environmental pressures, including dense dump distribution, intensive human activity, and low vegetation cover, jointly resulting in elevated risk levels. In contrast, upstream and peripheral zones had notably lower EERI values due to minimal surface disturbance, low population density, and limited human activity. Reclaimed areas generally showed moderate risk levels, indicating partial ecological improvement but also exposing deficiencies in runoff control and vegetation restoration. Overall, the EERI spatial pattern reflected the combined influence of land-use intensity, topographic features, and socioeconomic factors, providing spatially explicit guidance for targeted risk mitigation and ecological restoration planning.

To identify the key environmental risk factors driving spatial heterogeneity in watershed-scale eco-environmental risk, GRA and MDA were applied. Both methods are recognized for their robustness in handling multidimensional datasets and their suitability for assessing complex environmental systems with limited or uncertain information. In this study, nine indicators from the driving-force and pressure dimensions were selected for analysis (Fig. 7).

GRA was used to quantify the degree of association between each individual indicator and the overall EERI. The analysis involved constructing a reference sequence based on EERI values and comparing it with sequences derived from each standardized indicator, followed by the calculation of grey relational coefficients and association rankings. The results showed that dump density, runoff volume, and hydrogeological conditions had the highest association rankings, indicating strong correlations with EERI. Meanwhile, Manhattan distance analysis measured the spatial dissimilarity between each standardized indicator and the EERI surface by calculating the sum of absolute differences across all grid cells. Indicators such as population density, GDP, dump density, and runoff volume exhibited the smallest total Manhattan distances, suggesting close alignment of their spatial patterns with the EERI surface. Integrating the ranking results from both methods revealed that dump density and runoff volume consistently ranked among the top contributors in both analyses, and thus were identified as the primary drivers of spatial variability in watershed-scale environmental risk.

Building on the observed temporal and spatial patterns, further analysis was conducted to explore the policy, industrial, and hydrological factors underlying these changes. The temporal trends observed in dump contraction and reclamation expansion from 2010 to 2025 reflect the combined influence of policy interventions, industrial restructuring, and restoration initiatives. The stable high dump footprint in 2010–2013 coincided with the commissioning of new coal preparation, coal chemical, and pithead power facilities in Gujiao City and the Xishan mining area, which increased solid waste generation and storage demand. The 2014 revision of the Coal Gangue Comprehensive Utilization Management Measures set clear requirements for the resource utilization of coal gangue and fly ash, and the 2015–2016 provincial and local campaigns targeting “scattered and polluting” enterprises removed open-air dumps and illegal disposal sites. These measures drove a marked reduction in dump areas and rapid growth in reclamation between 2013 and 2019. From 2019 to 2025, the implementation of the 14th Five-Year Plan Guidance on Comprehensive Utilization of Bulk Solid Waste and Shanxi’s pilot utilization projects promoted value-added uses such as ultrafine fly ash production and high-value coal gangue applications. Although these projects introduced limited new storage demand, strict reclamation requirements kept the total impacted area stable.

The spatial patterns observed are also strongly shaped by hydrological structure and topographic setting. Many dumps are located in valley bottoms or gentle slopes within sub-catchment convergence zones, where micro-topographic depressions channel runoff into shared gully networks. This connectivity accelerates pollutant mobilization and downstream transport during high-flow events, increasing cumulative impacts at confluence points (Chen et al., 2021). In several reclaimed sites, the absence of integrated drainage guidance or interception measures limits their ability to disrupt pollutant pathways, reducing the effectiveness of ecological restoration. These findings underscore the interconnected nature of environmental risks in coal-mining watersheds, where industrial layout, topography, and hydrology interact with policy measures to shape both spatial patterns and risk intensity (Bao et al., 2021). Effective management should integrate dump siting and reclamation design with catchment runoff characteristics, ensuring that restoration incorporates hydrologically informed infrastructure to achieve both waste footprint reduction and watershed ecological security. The integration of high-resolution remote sensing and visual interpretation has proven effective for monitoring these dynamics, providing a spatially explicit evidence base for targeted land-use planning and risk mitigation (Garrett et al., 2000; Zhou et al., 2023; Giri et al., 2021).

The spatial heterogeneity of EERI results from the interplay of natural factors and anthropogenic pressures, interpreted through the Drivers-Pressures-State-Impact-Response (DPSIR) framework. In the driver dimension, precipitation regimes and topographic relief regulate hydrological connectivity and pollutant transport, with steep valleys and convergence zones facilitating runoff generation and directing contaminants along terrain-controlled flow paths, thereby increasing risk in lower reaches. In the pressure dimension, risk intensifies where anthropogenic stressors overlap. High dump density, elevated population concentration, and proximity to major drainage corridors amplify pollutant loading and accelerate the transport of leachate and surface runoff. In highly responsive terrain units, these pressures lead to rapid convergence of contaminants along flow paths, forming spatial clusters of elevated EERI that align with identified hotspots. The state dimension mediates these pressures: sub-watersheds with intact vegetation cover, heterogeneous and stable landscape patterns, and healthy vegetation indices demonstrate a stronger capacity to buffer and attenuate pollutants, resulting in consistently lower EERI values. In contrast, degraded vegetation and fragmented landscapes diminish resilience. The impact dimension is evident in high-risk zones, where declines in NRSEI, deterioration of downstream water quality, and persistent heavy metal accumulation indicate multi-media and long-lasting ecological risks. The response dimension appears in areas with targeted reclamation and infrastructure improvements, where risk levels are generally moderate to low, demonstrating that well-designed interventions can effectively reduce pressures, mitigate impacts, and enhance watershed resilience (Fernández-Caliani et al., 2021).

The quantitative assessment of risk drivers using Grey Relational Analysis and Manhattan distance analysis further reinforces these interpretations. Both methods identified dump density and runoff volume as the most influential contributors to EERI variability, underscoring the dual role of pollutant sources and hydrological transport capacity in shaping spatial risk patterns. The prominence of hydrogeological conditions in GRA rankings suggests that subsurface characteristics also play a role in mediating contaminant infiltration and migration, while the high ranking of socio-economic indicators such as population density and GDP in Manhattan distance analysis reflects the compounding effects of human activity intensity. These results confirm that spatial risk hotspots emerge where pollutant generation potential and efficient transport pathways converge, and where socio-economic drivers amplify pressure on environmental systems (Teng et al., 2021).

The integration of spatial pattern analysis and quantitative driver identification highlights the need for multi-pronged management strategies. Risk mitigation should prioritize reducing dump density through consolidation or reclamation, interrupting hydrological connectivity in high-runoff zones via engineered drainage controls, and strengthening vegetation restoration in erosion-prone slopes to disrupt contaminant transport. Additionally, land-use planning should address the spatial co-location of high-intensity human activities and sensitive hydrological corridors. The use of high-resolution remote sensing and spatially explicit modelling in this study demonstrates the value of combining spatial analysis with quantitative attribution techniques to guide targeted, evidence-based interventions in watershed-scale environmental risk management.

The management recommendations informed by the identification of dominant environmental risk drivers emphasize the need for differentiated and targeted ecological governance strategies in coal-based solid waste dump watersheds. Priority should be given to mitigating pollutant sources through accelerated consolidation and reclamation of waste dumps, alongside the implementation of engineered controls such as impermeable liners and leachate collection systems to limit contaminant release. Hydrological interventions are crucial to disrupt pollutant transport pathways, including drainage control, slope stabilization, and the establishment of vegetative buffer zones, while natural treatment systems like constructed wetlands and retention basins can further enhance water quality by filtering runoff (Zhang et al., 2021).

In addressing socio-economic pressures, land-use planning must strategically regulate high-intensity human activities near sensitive hydrological corridors to reduce pollutant flushing associated with urbanization and industrial development. Ecosystem restoration efforts focused on vegetation rehabilitation, particularly using native and deep-rooted species on erosion-prone slopes and riparian zones, are vital to strengthening the landscape’s buffering capacity and improving pollutant attenuation (Bashir et al., 2024). The success of these management measures relies on coordinated actions among environmental agencies, industry stakeholders, and local governments, supported by continuous monitoring through high-resolution remote sensing and spatial modelling. This integrated, adaptive management approach is essential for sustaining watershed resilience and effectively mitigating environmental risks over the long term.