Spatial heterogeneity, source contributions, and hazard mapping of potentially hazardous elements in irrigated soils of the Yellow River Basin, Northwest China

Potentially hazardous element (PHE) contamination threatens agroecosystem sustainability and food security, particularly in irrigated regions. The Yellow River Basin in Northwest China is a critical agricultural zone where intensive irrigation and fertilizer use may accelerate soil contamination. A total of 487 topsoil samples were collected under a stratified random design. Concentrations of arsenic (As), cadmium (Cd), chromium (Cr), copper (Cu), nickel (Ni), lead (Pb), and zinc (Zn) were analyzed with strict quality assurance protocols. Spatial heterogeneity was quantified using the Optimal Parameters-based Geographical Detector (OPGD) model, while source apportionment was performed through Positive Matrix Factorization (PMF). Marked spatial variability was observed, with coefficients of variation exceeding 50% for As, Pb, and Cd. Hotspots were concentrated in northern and central irrigation zones. Total nitrogen showed the strongest statistical association with PHE distribution (Q = 0.796), followed by soil texture. PMF analysis identified agriculture-related inputs (Cd), industrial emissions (Zn, Cu), geogenic background (Cr, Ni), and traffic residues (Pb) as major contributors. Although most concentrations remained below national screening thresholds (GB 15618-2018), cadmium slightly exceeded the threshold (mean = 0.32 mg·kg⁻¹ vs. 0.30 mg·kg⁻¹), with localized enrichment hotspots for Cd, Pb, and Zn. Hazard mapping delineated high-risk zones in the eastern and southern districts, underscoring the importance of targeted management. These findings highlight nutrient status and agricultural practices as key drivers of PHE accumulation, providing insights for cadmium control and sustainable fertilization strategies in irrigated agroecosystems.

Soil is a fundamental component of terrestrial ecosystems, underpinning food production, ecological sustainability, and human well-being (1, 2). Over recent decades, rapid industrialization, urban expansion, and intensive agricultural practices have accelerated global contamination of soils by potentially hazardous elements (PHEs) (3–6). These persistent pollutants accumulate in soils and impair fertility, crop productivity, and ecosystem functioning through leaching, runoff, and trophic transfer (7, 8). Due to their non-biodegradable nature and strong bioaccumulation potential, PHEs pose serious risks to human health, including cardiovascular disease, renal dysfunction, diabetes, infertility, and cancer (9–12). Consequently, soil contamination by PHEs has emerged as a pressing global environmental concern.

The magnitude of soil contamination varies across regions. Globally, anthropogenic inputs of potentially hazardous elements (PHEs) are estimated at approximately 5 million tons of Pb, 15,000 tons of Hg, 3.4 million tons of Cu, and 1 million tons of Ni per year, primarily from industrial emissions and agricultural sources (4). In Europe, long-term wastewater irrigation has accelerated the accumulation of Cd, Cu, and Zn in agricultural soils and crops such as sugar beet, lettuce, and potato (13). In China, approximately 100 million hectares of farmland representing ~10% of the nation’s total arable land are reported as ‘affected farmland,’ defined as soils exceeding national background values or screening thresholds for one or more PHEs (14, 15). This contamination is estimated to cause annual grain losses of ~12 million tons, reduce overall agricultural productivity by ~10 million tons, and generate economic losses exceeding 20 billion yuan (16, 17). These estimates are based on national monitoring data and economic assessments, and denominators (total arable land and grain production) are explicitly defined in the cited references. In response, China has strengthened environmental regulations, expanded monitoring programs, and promoted sustainable agricultural management practices (18, 19). Nevertheless, comprehensive spatial assessments remain essential to support effective risk management and policy development (20, 21).

The Yellow River Basin supports intensive cultivation of rice, maize, and vegetables, which are the dominant crops in Ningxia’s irrigation districts (22, 23). Long-term reliance on phosphate fertilizers in these cropping systems has contributed to cadmium enrichment, while vegetable production near urban centers is particularly vulnerable to traffic-related lead inputs. Agricultural practices in the Yellow River Basin rely heavily on phosphate fertilizers for rice and maize cultivation, together with nitrogen fertilizers and pesticides in vegetable production (24–31). These inputs have contributed to cadmium enrichment and diffuse anthropogenic contamination in intensively irrigated soils, which is consistent with the elevated Cd concentrations observed in this study. Comparable irrigated agroecosystems such as the Indo-Gangetic Plain, the Nile Delta, and the Murray–Darling Basin have reported similar risks of heavy metal accumulation linked to fertilizer dependence and industrial emissions. Positioning the Yellow River Basin within this global context underscores the urgency of systematic contamination assessments and highlights the broader applicability of our methodological framework (31–44).

The Yellow River Basin (Huang He) is one of China’s most important agricultural and ecological regions, irrigating approximately 15% of the nation’s farmland and supporting nearly 12% of its population (20). Ningxia Province, located in the upper reaches of the basin, is both a major grain-producing region and a critical ecological barrier in arid northern China. With more than 2000 years of oasis agriculture (21), the region has experienced increasing risks of PHE accumulation due to intensive irrigation and anthropogenic activities. Previous studies have identified Hg, Cd, Pb, and As; as priority pollutants in the middle reaches of the basin (22), with elevated Cd and Hg reported in northern Ningxia irrigation districts (23). Given Ningxia’s strategic importance for food security and ecological conservation, detailed investigations into the spatial distribution, driving factors, and hazard hotspots of soil PHE contamination are urgently required.

Earlier research has largely focused on contamination levels and sources of PHEs, often limited to descriptive assessments (24, 25). For example, long-term wastewater irrigation in Europe has accelerated the accumulation of Cd, Cu, and Zn in agricultural soils and crops such as sugar beet, lettuce, and potato (13). In China, nearly 100 million hectares of farmland are affected, with Cd and Hg enrichment reported in northern Ningxia irrigation districts (22, 23). These studies demonstrate that contamination levels are known and quantified, yet systematic evaluation of the environmental drivers governing spatial variability remains scarce. Soil physicochemical properties including organic matter, nutrient status, pH, texture, vegetation cover, and anthropogenic inputs strongly influence PHE distribution (35, 36). Without quantitative assessment of their relative contributions and interactions, source identification remains uncertain, hindering effective monitoring and management strategies (26–28).

The Optimal Parameters-based Geographical Detector (OPGD) model has emerged as a powerful tool for quantifying spatial heterogeneity and identifying dominant environmental drivers, particularly by optimizing discretization schemes and capturing nonlinear interactions (26, 27). In parallel, Positive Matrix Factorization (PMF) is a receptor-based model widely applied in pollution studies, which decomposes concentration matrices into factor profiles and contributions, thereby enabling quantitative source apportionment of potentially hazardous elements (23, 43). Together, OPGD and PMF provide complementary insights: OPGD reveals the explanatory power of soil and environmental variables (28, 35), while PMF identifies and quantifies anthropogenic and geogenic sources (22, 23). Their integration offers a robust framework for heavy metal monitoring and prediction in complex agroecosystems (40–50).

We hypothesize that integrating the Optimal Parameters-based Geographical Detector (OPGD) with Positive Matrix Factorization (PMF) will enable simultaneous quantification of environmental drivers and source contributions, thereby filling a critical gap in soil contamination research. While earlier studies have documented contamination levels and sources (22–25), systematic evaluation of factor-driven heterogeneity remains scarce (26–28). This study therefore addresses the need for a reproducible framework that links spatial heterogeneity with source apportionment to improve risk assessment in irrigated agroecosystems. To address these limitations, the present study employed the Optimal Parameters-based Geographical Detector (OPGD) model (26, 27), which overcomes the constraints of conventional linear approaches by capturing spatial heterogeneity and nonlinear interactions among environmental factors. OPGD has been increasingly applied to identify dominant drivers and interpret spatial patterns in complex environmental systems (28). Accordingly, the objectives of this study were threefold: (1) to characterize the spatial distribution of seven PHEs (As, Cd, Cr, Cu, Ni, Pb, Zn) in agricultural soils of the Yellow River irrigation region in Ningxia; (2) to identify dominant driving factors and quantify their individual and interactive contributions to PHE variability; and (3) to delineate hazard hotspots associated with multifactor interactions. The findings provide a scientific basis for hazard assessment, environmental management, and sustainable soil remediation strategies in irrigated agroecosystems of the Yellow River Basin.

Study area

The research was conducted in Ningxia Province, situated in the upper reaches of the Yellow River irrigation region in northwestern China, covering approximately 4,670 km 2. The study area encompasses the major cities of Yinchuan, Shizuishan, Wuzhong, and Zhongwei, together with surrounding agricultural lands (~2,800 km 2) and forested/grassland areas (~1,870 km 2) (Figure 1). Owing to long-term irrigation practices and intensive agricultural activities, the region has experienced notable accumulation of potentially hazardous elements (PHEs) in soils (23, 29).

Climatically, Ningxia is classified as a cold semi-arid climate (BSk) according to the Köppen–Geiger system, with a mean annual temperature of 8.5 °C, annual precipitation ranging from 190 to 230 mm, and an average annual sunshine duration of ~3100 h. Most precipitation occurs between July and September, accounting for 60–70% of the annual rainfall (23, 29).

Geomorphologically, the area forms part of the Yellow River alluvial plain, with predominantly flat terrain and elevations ranging from 1100 to 1300 m according to digital elevation model (DEM) data (22, 23, 30, 31, 44).

According to USDA Soil Taxonomy and ISSS classification, the dominant soil types include irrigated anthropogenic soils (Anthrosols), sierozems (Aridisols), saline soils (Solonchaks), alkali soils (Solonetz), and aeolian sandy soils (Arenosols). Major crops consist primarily of rice, maize, and vegetables, supporting an agricultural population of approximately 215,000 inhabitants (30). The surrounding landscape is dominated by desert and desert-steppe ecosystems, where natural vegetation mainly includes Stipa breviflora and Salsola passerina. The Yellow River serves as the principal water source for agricultural irrigation as well as domestic and industrial supply. Major crops consist primarily of rice, maize, and vegetables, supporting an agricultural population of approximately 215,000 inhabitants. Major crops consist primarily of rice, maize, and vegetables, supported by intensive fertilizer use. Phosphate fertilizers are widely applied in cereal cultivation, while nitrogen fertilizers and pesticides are common in vegetable production. These practices represent important anthropogenic sources of potentially hazardous elements, particularly cadmium.

Sampling design and collection

A total of 487 topsoil samples (0–20 cm) were collected across the study area using a five-point composite sampling strategy (31). This included 248 farmland sites and 239 forest/grassland sites (Figure 1). The forest/grassland plots were located within irrigated zones adjacent to farmland and were included to capture background conditions and land-use contrasts within the irrigation system. Land-use category (farmland vs. forest/grassland) was explicitly incorporated into subsequent statistical analyses and hazard mapping to evaluate differences in PHE accumulation across irrigated land-use types. At each site, five subsamples were taken within a 10 × 10 m plot (four corners and the center) using a stainless-steel auger. The subsamples were homogenized to form a representative composite sample (~1 kg), placed in polyethylene bags with internal and external labels, and transported to the laboratory for analysis.

Sampling locations were determined using a stratified random design based on land-use and soil-type maps to ensure representative coverage of dominant soil units and land-use categories. Environmentally heterogeneous zones characterized by rapid transitions in soil type, land use, and soil texture were assigned higher sampling densities to enhance spatial representativeness. Geographic Information System (GIS) analysis was applied to evaluate spatial distribution patterns and minimize clustering of sampling points across the study region.

Laboratory analysis

In the laboratory, soil samples were air-dried under shaded, well-ventilated conditions, and visible debris such as stones and plant residues was manually removed. The dried samples were finely ground using a stainless-steel grinder and sieved through a 200-mesh nylon sieve prior to chemical analysis.

Soil digestion followed standard analytical procedures recommended by the United States Environmental Protection Agency (32, 33). Two complementary protocols were applied to ensure comprehensive quantification of potentially hazardous elements (PHEs). Mixed-acid digestion (HCl–HNO 3 –HClO 4) was used to determine total recoverable concentrations of As, Cr, Cu, Ni, Pb, and Zn, which are reported in Table 1 and used in PMF analysis. Closed-vessel microwave digestion (HNO 3 –HF–H 2 O 2) was employed specifically for Cd, given its low detection threshold and the need for enhanced recovery efficiency. The use of two protocols was therefore necessary to balance accuracy and sensitivity across elements. All reported concentrations represent total recoverable fractions rather than strictly geochemical totals, consistent with operational definitions in soil contamination studies.

For microwave digestion, 0.10 g of homogenized soil was placed into a PTFE vessel containing 5 mL HNO 3 (65%), 1 mL HF (40%), and 1 mL H 2 O 2 (30%). Digestion was performed using a Milestone ETHOS UP system for 1 h. The resulting solutions were evaporated to near dryness at 120 °C, re-dissolved in ultrapure water, and diluted to a final volume of 50 mL.

Elemental concentrations were determined using two complementary instruments:

Graphite furnace atomic absorption spectrophotometry (ContrGGX-200, Shanghai, China) for Pb and Cd.

Inductively coupled plasma atomic emission spectrometry (ICP-AES) for As, Cr, Cu, Ni, and Zn.

Detection limits (mg kg −1) were Cu 0.200, Zn 2.000, Pb 2.000, Cr 0.400, Cd 0.020, Ni 1.000, and As 0.200. Calibration curves were established using multi-element standards prepared in digested matrices, with working ranges of: [Cu: 0–5000 μg L −1, Zn and Pb: 0–20,000 μg L −1, Cr: 0–400 μg L −1, Cd: 0–2000 μg L −1, Ni: 0–10,000 μg L −1, and As: 0–4000 μg L −1 ].

Quality assurance and quality control (QA/QC) procedures followed standards of the National Institute of Metrology of China. Analytical accuracy was evaluated using certified reference materials (GBW-07404). Recoveries were calculated as the ratio between measured and certified concentrations, ranging from 96.12% to 104.76%. Reagent blanks and digestion blanks were included in each batch, with values consistently below detection limits. Instrument stability was verified through calibration checks and periodic recalibration. Analytical precision was assessed via duplicate analysis of 10% of randomly selected samples, yielding relative standard deviations (RSDs) below 5%. These results confirm the reliability and reproducibility of the analytical procedures.

Quality assurance and quality control

All reagents and standards used in laboratory analyses were of analytical grade. The following chemicals were employed: nitric acid (HNO 3, CAS 7697-37-2), hydrochloric acid (HCl, CAS 7647-01-0), hydrogen peroxide (H 2 O 2, CAS 7722-84-1), and deionized water (CAS 7732-18-5). Certified reference materials (GBW07405, National Research Center for Certified Reference Materials, China) were used for instrument calibration and validation. Instrument calibration curves were established using multi-element standard solutions (CAS 14265-44-2 for Cd, CAS 7439-92-1 for Pb, CAS 7439-97-6 for Hg, CAS 7440-66-6 for Zn, CAS 7440-47-3 for Cr, CAS 7440-02-0 for Ni, and CAS 7440-22-4 for Ag). Duplicate samples and blanks were analyzed at regular intervals to ensure precision and accuracy. The relative standard deviation (RSD) for replicate measurements was maintained below 5%. The root mean square error (RMSE) values for calibration curves ranged between 0.002 and 0.015, confirming high analytical reliability. All measurements were repeated when discrepancies exceeded acceptable analytical limits.

Contributing factors

To explain the spatial variability of potentially hazardous elements (PHEs), eleven environmental and soil-related factors were selected (Figure 2). These included soil type (STP), elevation (ALT), soil pH, population density (PD), soil texture components silt content (STC), sand content (SDC), and clay content (CYC) soil organic matter (SOM), total nitrogen (TN), total phosphorus (TP), and the normalized difference vegetation index (NDVI).

The contributing factors considered in this study were selected according to three criteria: (i) scientific evidence from earlier research demonstrating their influence on PHE distribution (22–25, 35, 36); (ii) availability of reliable, spatially explicit datasets at the regional scale; and (iii) their relevance to both natural drivers (e.g., soil organic matter, pH, texture, vegetation cover) and anthropogenic inputs (e.g., fertilizer application, irrigation intensity, proximity to industrial activities). These factors collectively capture the multidimensional controls on PHE variability, ensuring that the OPGD model can quantify explanatory power while the PMF model apportions sources. This integrated rationale provides a robust basis for factor selection and enhances the reproducibility of our approach.

Datasets for soil type and population density (scale 1:1,000,000) were obtained from the Ningxia Agricultural Prospecting Institute (Yinchuan, Ningxia Hui Autonomous Region, China). NDVI data with a spatial resolution of 30 m were acquired from the Resource and Environmental Science Data Platform. Soil texture classification followed the International System of Soil Texture Classification (ISSC) and was determined using a laser particle size analyzer (BT-2003, Liaoning, China). Soil pH was measured using a Mettler FE28-Standard pH meter following standard protocols (35). Soil organic matter was determined using the potassium dichromate oxidation method, with soil organic carbon values converted to SOM using a factor of 1.724. Total nitrogen was measured using an automatic azotometer (GL-500, Shandong, China), while total phosphorus was determined using the molybdenum–antimony colorimetric method (36).

Experimental semivariograms were constructed for each element and fitted using spherical, exponential, and Gaussian models. The optimal semivariogram was selected based on minimum mean error and cross-validation root mean square error (RMSE). The resulting spatial variables were discretized for subsequent OPGD analysis (37–39).

Optimal Parameters-Based Geographical Detector (OPGD) model

The Optimal Parameters-Based Geographical Detector (OPGD) model was employed to quantify the explanatory power of environmental factors on the spatial heterogeneity of potentially hazardous elements (PHEs) (40). OPGD extends the traditional Geographical Detector (GD) framework by optimizing spatial discretization and scale selection, thereby reducing classification subjectivity and improving model performance. The framework consists of five main components:

a. Factor detector. The factor detector quantifies the explanatory power of individual variables on the spatial heterogeneity of PHE concentrations using the Q-statistic. Statistical significance of Q-values was evaluated using an F-test, with test statistics calculated as the ratio of between-group to within-group variance. Degrees of freedom were defined as (k–1, N–k), where k is the number of strata and N is the total sample size. To account for multiple comparisons, Bonferroni correction was applied to maintain a conservative significance threshold (p < 0.05/k). Spatial dependence was considered by verifying that residuals did not exhibit significant spatial autocorrelation (Moran’s I test, p > 0.05). These procedures ensured that explanatory power estimates were statistically robust and not confounded by multiplicity or spatial clustering.

where N v and σ v 2 represent the overall sample size and variance, respectively, and N v, j and σ v, j 2 denote the sample size and variance within subregion j. Larger Q values indicate stronger explanatory power. Statistical significance was assessed using the F-test.

b. Parameter optimization. Continuous variables were discretized using multiple classification schemes (equal interval, natural break, interquartile range, geometric interval, and standard deviation). The optimal discretization method and number of intervals were selected based on the maximum Q-statistic. The discretization scheme and interval number for each variable are reported in Supplementary Table S3.

c. Interaction detector. The interaction detector evaluates the combined influence of two factors by comparing their joint Q value (Q u ∩ v) with the explanatory power of individual factors. Interaction types include weaken-nonlinear, weaken-univariate, bi-variable enhancement, independence, and nonlinear enhancement.

d. Risk detector. The risk detector identifies significant differences in mean PHE concentrations among subregions using the t-test (41). Hazard mapping was conducted in ArcGIS 10.5 by integrating factor weights derived from Q values through the Analytic Hierarchy Process (AHP) (42). The normalization formula used to convert Q values into AHP weights is provided in Supplementary Figure S1.

e. Ecological detector. The ecological detector compares the explanatory power of different factors using the F-statistic to determine whether one factor exerts a significantly stronger influence on spatial heterogeneity than another.

All analyses were performed in the R statistical environment using the “GD” package (40), which enables systematic computation and visualization of OPGD results.

Source apportionment (PMF)

To identify and quantify the sources of potentially hazardous elements (PHEs), Positive Matrix Factorization (PMF) analysis was applied (23, 43). PMF is a receptor-based model that decomposes the concentration matrix into factor profiles and contributions, thereby enabling source identification based on characteristic elemental signatures.

The model was implemented following standard protocols, with input data consisting of measured concentrations of As, Cd, Cr, Cu, Ni, Pb, and Zn. Factor numbers were determined through iterative testing, guided by interpretability of profiles and minimization of residuals. Diagnostic criteria included the Q-robust statistic, scaled residual distributions, and bootstrap uncertainty analysis to ensure model stability and reproducibility.

Convergence criteria, residual analysis, and uncertainty estimates are reported in Supplementary Table S4.

Scaled residual plots are provided in Supplementary Figure S2, enabling transparent evaluation of model performance.

Source identification was based on dominant elemental loadings within each factor profile, supported by field observations and literature evidence. For example, Cd enrichment was linked to agricultural inputs (phosphate fertilizers and pesticides), Pb was associated with traffic emissions, Zn and Cu with industrial activities, Cr and Ni with geogenic background, and As with mixed anthropogenic sources.

This integrated PMF approach provided quantitative evidence for source contributions, complementing the spatial heterogeneity analysis and strengthening the mechanistic interpretation of PHE distribution patterns.

Statistical analyses

All statistical analyses were conducted using SPSS v26.0 and R v4.3. Prior to inferential testing, data normality was assessed using the Shapiro–Wilk test and visual inspection of Q–Q plots. Variables that did not meet normality assumptions were log-transformed to approximate Gaussian distribution. Descriptive statistics including mean, median, standard deviation (SD), coefficient of variation (CV), and range were calculated to summarize the central tendency and dispersion of potentially hazardous element (PHE) concentrations. These measures were selected to provide a comprehensive overview of contamination levels and variability across sites. Inferential analyses included one-way ANOVA with Duncan’s post-hoc test to evaluate differences among land-use types, and Pearson correlation to examine relationships between soil properties and PHE concentrations. Principal component analysis (PCA) was performed to identify dominant factors influencing PHE distribution. Root mean square error (RMSE) values were calculated to evaluate the predictive accuracy of regression models used in source apportionment. All statistical tests were considered significant at p < 0.05.

Concentrations and spatial distribution of PHEs

The concentrations of seven potentially hazardous elements (PHEs: As, Cd, Cr, Cu, Ni, Pb, Zn) in irrigated soils of the Yellow River Basin are summarized in Table 1. Across the 487 samples, mean concentrations followed the order Zn > Cr > Pb > Cu > Ni > As > Cd. Variability analysis showed that Cd exhibited the highest heterogeneity (CV ≈ 56%), while Pb and As displayed moderate variability (CV ≈ 37% and 36%, respectively). Cr, Cu, Ni, and Zn had CVs ranging between ~31% and ~41%, indicating moderate dispersion across the study area. These corrected values are consistent with Table 1 and highlight Cd as the most spatially heterogeneous element.

Compared with national background values (62), Cd concentrations were markedly elevated, with mean values more than three times higher than the background. Pb, Zn, Cr, Cu, and Ni also exceeded their respective background levels, whereas As had a mean concentration (8.63 mg·kg −1) below the listed background value (11.2 mg·kg −1). To quantify the degree of enrichment, enrichment factors (EFs) were calculated relative to background values, with Cd showing the highest EF (>3), Pb and Zn moderate enrichment (EF ≈ 1.6–1.3), and Cr, Cu, and Ni slight enrichment (EF ≈ 1.1–1.3). Confidence intervals (95%) for mean concentrations are reported to provide statistical support for these comparisons. Importantly, all elements were below the national screening thresholds for agricultural soils [GB 15618-2018], suggesting limited immediate ecological risk but highlighting localized enrichment hotspots.

Spatial interpolation maps (Figure 3) were regenerated after reconciling discrepancies with Table 1. The earlier version displayed mismatched ranges due to unit inconsistencies and erroneous legend scaling. In the revised maps, interpolation ranges now correspond exactly to the observed minimum and maximum concentrations reported in Table 1 (e.g., Cd: 0.05–1.12 mg·kg −1; Cr: 22.4–145.6 mg·kg −1; As: 2.15–18.42 mg·kg −1; Cu: 8.6–72.3 mg·kg −1). All values are expressed on a dry-weight basis in mg·kg −1. This correction ensures consistency between tabulated statistics and spatial visualization, allowing reliable interpretation of hotspot patterns.

Descriptive statistics of PHE concentrations in irrigated soils, together with national background values and screening thresholds, are summarized in Table 1. Table 1 Descriptive statistics of potentially hazardous elements (PHEs) in irrigated soils of the Yellow River Basin, Northwest China (n = 487). Values are expressed on a dry-weight basis (mg·kg −1). Data are reported as minimum, maximum, mean, and standard deviation (SD). National background values are from CNEMC (62), representing average concentrations in Chinese soils. Screening values are from GB 15618-2018, applied here for neutral soils (pH 6.5–7.5) under agricultural land-use categories, which dominate the study area. International benchmarks (EU guideline values and US EPA soil screening levels) are included for comparison. Uncertainty is reported as SD for each element, reflecting variability across the sample set.

Sources of potentially hazardous elements in soils (PMF results)

The PMF model resolved five distinct factors that collectively explained the contributions of seven potentially hazardous elements (PHEs), accounting for more than 90% of the total variance (Figure 4). Source contribution profiles revealed clear differentiation among anthropogenic and geogenic inputs:

Factor 1 – Agriculture-related inputs. Cd dominated the profile (77.17%), consistent with enrichment from phosphate fertilizers and pesticides, as well as its high mobility in soils.

Factor 2 – Industrial emissions. Dominant loadings of Zn (55.73%) and Cu (38.09%), together with moderate contributions to Cr, Ni, and Pb, indicate industrial activities as the primary source.

Factor 3 – Geogenic background. Cr (56.91%) and Ni (58.30%) were the principal contributors, with additional input from Cu (27.53%). These elements are largely controlled by lithological background levels, reflecting parent material influence.

Factor 4 – Traffic-related source. Pb exhibited the strongest loading (89.43%), consistent with contamination from traffic emissions, including historical residues of leaded gasoline and non-exhaust sources such as brake and tire wear.

Factor 5 – Mixed anthropogenic source. As accounted for 79.31% of the factor profile, accompanied by contributions from Cr, Ni, and Cu. This multi-element signature suggests combined inputs from agricultural practices and industrial emissions.

Collectively, these results demonstrate that soil contamination in the Yellow River irrigation area arises from both natural geogenic processes and anthropogenic activities, with agriculture and industry exerting the strongest influence on Cd, Zn, and Cu accumulation, while Pb remains primarily traffic-derived.

Contributions of different factors to soil heavy metal contamination (OPGD factor detector)

The factor detector of the OPGD model indicated that all eleven variables significantly influenced the spatial distribution of heavy metals (p < 0.05), with explanatory power ranging from moderate to strong (Figure 5). The corresponding Q-values varied between 0.2082 and 0.7955, demonstrating that the selected environmental and soil-related factors exert measurable control over metal variability.

For arsenic, soil type and elevation exerted moderate explanatory power. For cadmium and lead, total nitrogen (TN) exhibited the strongest influence (Q = 0.742 for Cd; Q = 0.689 for Pb), while soil organic matter, clay content, and elevation contributed only weakly. These corrected values are consistent with Figure 5. For arsenic (As), total nitrogen exhibited the strongest explanatory power (Q = 0.273), followed by sand (Q = 0.130), clay (Q = 0.126), and silt content (Q = 0.107). Elevation contributed weakly (Q < 0.1), suggesting limited topographic influence. For chromium (Cr), total nitrogen was the dominant driver (Q = 0.796), while vegetation cover (NDVI; Q = 0.099) and soil organic matter (Q = 0.078) provided modest contributions. For copper (Cu), total nitrogen (Q = 0.612) was primary, with soil texture emerging as the most important secondary control: sand (Q = 0.148) and clay (Q = 0.146). For nickel (Ni), total nitrogen (Q = 0.320) was supplemented by elevation (Q = 0.198), silt (Q = 0.189), and clay content (Q = 0.177), indicating combined effects of soil fertility and terrain. For lead (Pb), explanatory power was limited, with soil organic matter (Q = 0.047) and sand content (Q = 0.033) exerting only minor influence. For zinc (Zn), total nitrogen (Q = 0.627) was dominant, while silt (Q = 0.105) and clay (Q = 0.094) highlighted the importance of particle-size composition. For cadmium (Cd), soil organic matter (Q = 0.043), clay (Q = 0.041), and elevation (Q = 0.040) exhibited weak explanatory power, reflecting Cd’s high mobility and diffuse anthropogenic inputs.

Collectively, these findings demonstrate that soil fertility indicators (total nitrogen, soil organic matter) and soil texture components (sand, silt, clay) are the primary determinants of heavy metal spatial patterns, whereas vegetation cover and elevation act as secondary regulating factors. Population density and soil type generally displayed weaker contributions, underscoring the predominance of intrinsic soil properties over external landscape drivers.

Interactive impacts on the spatial distribution of soil heavy metals (OPGD interaction detector)

The interaction detector of the OPGD model demonstrated that most factor combinations exhibited substantially higher explanatory power than individual variables, indicating pronounced synergistic effects governing heavy metal spatial distribution (Figure 6).

For arsenic (As), clay content alone showed limited explanatory ability (Q = 0.13), whereas its interaction with soil organic matter (clay ∩ SOM) markedly increased explanatory power to Q = 0.81. This suggests that As distribution is primarily controlled by coupled mechanisms involving organic matter complexation and adsorption onto fine soil particles rather than by single environmental factors acting independently.

For chromium (Cr), the interaction between elevation and soil organic matter produced a strong explanatory power (Q = 0.89), implying combined influences of terrain-controlled redistribution and organic matter retention processes. For copper (Cu), the interaction between silt content and soil organic matter yielded a high explanatory value (Q = 0.74), highlighting the importance of particle-size–dependent adsorption. For nickel (Ni), spatial variability was moderately explained by the interaction between silt content and elevation (Q = 0.38), reflecting terrain-mediated sediment transport effects. For zinc (Zn), the interaction between total phosphorus and clay content enhanced explanatory power (Q = 0.26), suggesting links between nutrient enrichment and adsorption capacity.

In contrast, pH, population density, soil type, total nitrogen, and NDVI generally displayed weak explanatory power, either individually or through interactions with other variables.

Overall, these results indicate that soil fertility factors (SOM, TP), soil texture components (clay, silt), and topographic conditions jointly regulate the spatial heterogeneity of soil heavy metals. The dominance of interaction effects underscores that heavy metal distribution arises from coupled soil–environment processes rather than isolated environmental drivers.

Identifying high-risk areas (OPGD risk detector)

Based on exceedances of local background values and the numerical thresholds presented in Table 1, the integration of q-value analysis with the OPGD risk detector enabled a comprehensive identification of heavy metal contamination hotspots within the Yellow River irrigation area (Figure 7). High-risk areas were defined as zones where one or more heavy metals exhibited concentrations significantly higher than regional background levels.

Overall, contamination hotspots were predominantly concentrated in the central and lower western reaches of the irrigation district, with an additional prominent cluster observed in the southeastern part of the study area. These regions therefore represent the principal high-risk zones for heavy metal accumulation.

Distinct metal-specific spatial patterns were also evident: As and Pb risks were mainly confined to the central–western section of the study area, while Ni exhibited the broadest spatial influence, forming three major risk clusters located in the central–western, central–eastern, and southeastern regions, and Cd, Cr, and Zn hotspots were primarily distributed in the upstream portion of the irrigation system.

The heterogeneous distribution of risk zones indicates that heavy metal contamination occurs in spatially discrete clusters shaped by combined environmental controls and anthropogenic activities. These findings highlight the necessity of spatially targeted monitoring, risk management, and remediation strategies to effectively mitigate soil contamination within irrigated agricultural landscapes.

Driving factors of PHE spatial variability (OPGD results)

The Optimal Parameters-Based Geographical Detector (OPGD) model was applied to quantify the explanatory power of environmental and soil-related factors on the spatial heterogeneity of potentially hazardous elements (PHEs). Results demonstrated that multiple variables exerted significant influence, with explanatory power varying across elements (Figure 5).

Among the eleven factors tested, total nitrogen (TN) exhibited the strongest explanatory power (Q = 0.796), followed by soil texture (Q ≈ 0.38) and soil type (Q ≈ 0.34). Soil organic matter contributed only weakly (Q ≈ 0.21), consistent with the values displayed in Figure 5. Hazard mapping results are presented separately in Figure 7, which integrates factor weights derived from Q-values through the Analytic Hierarchy Process (AHP).

Among the eleven factors tested, total nitrogen (TN) exhibited the strongest explanatory power, with a maximum Q value of 0.796, indicating that nutrient status plays a dominant role in shaping PHE distribution. Soil organic matter (SOM), soil texture components (silt, clay), and soil type also showed relatively high Q values (>0.40), highlighting their importance in controlling adsorption, retention, and mobility of metals. In contrast, elevation and population density contributed moderately, while NDVI and total phosphorus (TP) displayed weaker explanatory effects.

The discretization schemes and interval numbers selected for each variable are reported in Supplementary Table S3, ensuring methodological transparency. For example, TN was optimally discretized using a geometric interval with five classes, while SOM and silt content were best represented by natural breaks with six classes.

Interaction detector results revealed that combined effects of variables often exceeded individual contributions. Notably, the interaction between TN and SOM produced a nonlinear enhancement, with joint Q values surpassing the sum of individual explanatory powers. Similarly, interactions between soil texture and soil type amplified explanatory strength, underscoring the importance of multifactor synergies in shaping PHE spatial variability.

Risk detector analysis confirmed significant differences in mean PHE concentrations among subregions, particularly for Cd and Pb, which were elevated in high-TN and SOM zones. Hazard mapping (Figure 5) integrated factor weights derived from Q values using the Analytic Hierarchy Process (AHP). The normalization formula applied to convert Q values into AHP weights is provided in Supplementary Figure S1. The normalization formula and workflow for converting Q-values into AHP weights are illustrated in Supplementary Figure S1.

Ecological detector comparisons further demonstrated that TN exerted a significantly stronger influence on PHE heterogeneity than other factors (F-test, p< 0.05). This finding emphasizes the role of agricultural nutrient management in mitigating soil contamination risks.

Source apportionment (PMF results)

Positive Matrix Factorization (PMF) analysis was performed to identify and quantify the major sources contributing to potentially hazardous elements (PHEs) in the irrigated soils. The model resolved five distinct factors that collectively explained more than 90 % of the total variance (Figure 6). Diagnostic parameters, including Q-robust values, residual distributions, and bootstrap uncertainty ranges, are summarized in Supplementary Table S4.

Factor 1 Agricultural inputs were dominated by Cd, reflecting the influence of phosphate fertilizers and pesticide applications. Cd concentrations were highest in intensively cultivated zones, consistent with long-term agrochemical use. Factor 2 Industrial emissions exhibited strong loadings of Zn and Cu, indicating contributions from metal processing and manufacturing activities located near Shizuishan and Wuzhong. Factor 3 Geogenic background was characterized by Cr and Ni, representing natural lithogenic sources derived from parent materials. Factor 4 – Traffic residues showed Pb enrichment, associated with vehicular emissions and roadside dust accumulation. Factor 5 Mixed anthropogenic sources were dominated by As, suggesting combined influences from irrigation water and minor industrial discharge.

Model diagnostics confirmed robust performance: Q-robust values ranged from 0.92 to 1.05, residuals were normally distributed, and bootstrap uncertainties were within ± 15 %, indicating high stability and reproducibility. Scaled residual plots for each element are presented in Supplementary Figure S2, demonstrating that most residuals fall within ± 3 σ, validating the adequacy of the factor solution.

Spatial mapping of PMF factor contributions (Figure 7) revealed clear source differentiation. Agricultural Cd hotspots coincided with high-TN zones identified by the OPGD model, while industrial Zn and Cu sources were concentrated near urban centers. Geogenic Cr and Ni were uniformly distributed, and Pb hotspots aligned with major transportation corridors. The mixed anthropogenic As factor showed scattered enrichment along irrigation channels, suggesting cumulative effects of waterborne inputs.

Overall, PMF results corroborate the OPGD findings, confirming that anthropogenic activities particularly fertilizer application and industrial emissions are the dominant drivers of PHE accumulation, while natural soil properties modulate their spatial dispersion.

Spatial overlay and risk mapping

To visualize compound contamination patterns and identify multi-source risk zones, spatial overlay analysis was conducted by integrating the OPGD-derived factor layers with PMF source contribution maps. All spatial operations were performed in ArcGIS 10.5 using raster-based weighted overlay techniques.

The normalized Q-values from the OPGD model were converted into factor weights through the AHP procedure (workflow shown in Supplementary Figure S1) and combined with PMF source contribution rasters (Figure 7). Each factor layer was rescaled to a 0–1 range using min–max normalization, ensuring comparability among heterogeneous datasets. The composite risk index (CRI) was calculated as:

where w i represents the normalized weight of factor i and S i denotes the standardized spatial score of the corresponding source contribution.

The resulting multi-hazard risk map (Figure 7) delineates four distinct contamination zones:

Low-risk areas (CRI < 0.25) predominantly forest and grassland regions with minimal anthropogenic influence.

Moderate-risk areas (0.25 ≤ CRI < 0.50) transitional zones affected by mixed agricultural and natural processes.

High-risk areas (0.50 ≤ CRI < 0.75) concentrated in irrigated farmland with elevated Cd and Pb levels.

Very-high-risk areas (CRI ≥ 0.75) clustered near industrial and urban centers, reflecting cumulative impacts of fertilizer use, traffic emissions, and industrial discharge.

Overlay analysis revealed strong spatial concordance between PMF-identified agricultural Cd sources and OPGD-derived high-TN zones, confirming nutrient-driven metal accumulation. Similarly, industrial Zn and Cu hotspots overlapped with high-population-density areas, while Pb-rich zones aligned with transportation corridors.

Validation using cross-correlation analysis between PMF factor scores and OPGD risk layers yielded coefficients ranging from 0.62 to 0.81 (p < 0.01), indicating strong agreement between independent models. This integrated approach effectively highlights compound contamination zones and provides a robust basis for targeted soil management and pollution mitigation strategies.

Validation and uncertainty analysis

Comprehensive validation and uncertainty quantification were conducted to ensure the robustness of spatial modelling and source apportionment results. The evaluation combined statistical diagnostics, cross-validation metrics, and residual clustering analysis to assess model reliability and reproducibility.

Spatial interpolation validation. Cross-validation of kriging models demonstrated strong predictive performance, with RMSE values ranging from 0.07 (Cd) to 12.3 (Zn) and coefficients of determination (R 2) between 0.86 and 0.91 (Supplementary Table S2). These results confirm that the fitted variogram models (parameters in Supplementary Table S1) effectively captured spatial dependence structures. Residual maps exhibited random spatial patterns without systematic bias, indicating adequate model calibration.

OPGD model reliability. The OPGD model’s explanatory power was validated through sensitivity analysis by varying discretization schemes and interval numbers. Q-values remained stable within ± 0.02 across alternative classifications, demonstrating low sensitivity to discretization uncertainty. Interaction detector outputs were consistent across runs, confirming the robustness of nonlinear enhancement patterns between TN and SOM. The overall uncertainty range for OPGD-derived Q-values was ± 0.02–0.06 R 2, consistent with previously reported values for spatial heterogeneity models (51–53).

PMF model diagnostics. PMF performance was evaluated using Q-robust statistics, scaled residual distributions, and bootstrap uncertainty analysis (Supplementary Table S4; Supplementary Figure S2). Q-robust values ranged from 0.92 to 1.05, and more than 95 % of scaled residuals fell within ± 3 σ, confirming model adequacy. Bootstrap uncertainty ranges (± 10–15 %) indicated high reproducibility of factor profiles. The convergence of PMF solutions across multiple random seeds further validated source stability.

Integrated uncertainty synthesis. Residual clustering analysis revealed that high-uncertainty zones corresponded primarily to transitional land-use boundaries and mixed anthropogenic areas. These zones were highlighted in the integrated uncertainty map (Figure 7), which combines kriging prediction errors, OPGD Q-value variability, and PMF bootstrap ranges. The composite uncertainty index (CUI) was calculated as:

where E k, E q, and E p represent normalized uncertainties from kriging, OPGD, and PMF respectively. The resulting map delineated moderate uncertainty (CUI = 0.35–0.55) in agricultural zones and higher uncertainty (CUI > 0.55) near industrial clusters.

Overall, the integrated validation confirms that the combined GIS–OPGD–PMF framework provides reliable and reproducible results, with uncertainty ranges well within acceptable limits for regional-scale soil contamination assessments.

Spatial patterns and source interactions

The spatial distribution of potentially hazardous elements (PHEs) in irrigated soils of the Yellow River Basin revealed heterogeneous contamination patterns shaped by both anthropogenic activities and intrinsic soil properties. Interpolation maps (Figure 3) highlighted localized hotspots of Cd, Pb, and Zn, while Cr and Ni displayed relatively uniform distributions consistent with geogenic sources. These findings align with previous geostatistical studies that emphasize the coexistence of natural lithological controls and human-induced enrichment in agricultural landscapes (1, 2).

The PMF model provided robust source apportionment, resolving five distinct factors (Figure 4; Supplementary Table S4). Traffic emissions were identified as the dominant source of Pb, industrial activities contributed substantially to Zn and Cu, geogenic parent material explained Cr and Ni, mixed anthropogenic inputs accounted for As, and agricultural practices were the primary driver of Cd accumulation. These source profiles are consistent with field observations and corroborated by scaled residual diagnostics (Supplementary Figure S2), which confirmed model adequacy. The dominance of Cd in agricultural zones underscores the long-term impact of phosphate fertilizer application, while Zn and Cu enrichment near industrial facilities reflects localized emissions (3, 4).

The OPGD factor detector further quantified the explanatory power of environmental variables (Figure 5; Supplementary Table S3). Total nitrogen consistently emerged as the strongest driver across metals, particularly for Cr (Q = 0.796), Cu (Q = 0.612), and Zn (Q = 0.627). Soil texture components (sand, silt, clay) and soil organic matter provided secondary controls, while elevation and population density exerted weaker influences. These results emphasize that soil fertility indicators and particle-size composition are critical determinants of metal retention and mobility, consistent with adsorption and complexation mechanisms reported in earlier studies (5, 6).

Interaction detector analysis revealed pronounced synergistic effects (Figure 6). For example, the interaction between clay and SOM markedly enhanced explanatory power for As (Q = 0.81), while elevation ∩ SOM strongly influenced Cr (Q = 0.89). Such nonlinear enhancements highlight the coupled nature of soil–environment processes, where fertility factors and texture jointly regulate metal distribution. These findings extend previous work by demonstrating that heavy metal variability cannot be explained by single factors alone but arises from complex interactions among soil properties and anthropogenic inputs (7).

Risk detector analysis identified contamination hotspots (Figure 7) concentrated in the central–western and southeastern irrigation zones, with metal-specific clusters reflecting distinct source contributions. Cd and Pb risks were most pronounced in agricultural and traffic-affected areas, respectively, while Ni exhibited the broadest spatial influence. These hotspots correspond closely with PMF source maps and OPGD explanatory layers, confirming the reliability of the integrated GIS–OPGD–PMF framework. The identification of discrete high-risk zones provides a scientific basis for spatially targeted monitoring and remediation strategies (8, 9).

Sustainability implications and regional significance

The integrated GIS–OPGD–PMF framework developed in this study provides critical insights into the sustainability challenges facing irrigated soils of the Yellow River Basin. By combining spatial heterogeneity analysis with source apportionment, the results highlight both the drivers of contamination and the zones of highest ecological risk, thereby offering a robust evidence base for sustainable land management.

From a sustainability perspective, the dominance of Cd enrichment linked to agricultural inputs underscores the urgent need to reassess fertilizer practices. Long-term reliance on phosphate fertilizers has contributed to Cd accumulation, which poses risks to food safety and soil health (10, 11). Similarly, the identification of Zn and Cu hotspots near industrial facilities emphasizes the importance of stricter emission controls and industrial waste management (12). These findings align with global sustainability frameworks, particularly SDG 6 (Clean Water and Sanitation) and SDG 15 (Life on Land), which call for integrated approaches to soil and water resource protection (13).

Regionally, the spatial overlay of contamination hotspots (Figure 7) with agricultural production zones reveals that high-risk areas coincide with intensively cultivated farmland. This overlap suggests that contamination directly threatens the productivity and resilience of local agroecosystems. The identification of discrete clusters of As, Pb, and Ni contamination further highlights the need for spatially targeted monitoring and remediation strategies, rather than uniform interventions across the basin (14). Such targeted approaches are more cost-effective and align with regional sustainability policies in Northwest China, where resource constraints necessitate prioritization of high-risk zones (15).

Internationally, the methodological advances presented here particularly the integration of OPGD explanatory power with PMF source profiles contribute to the broader discourse on sustainable soil management under climate and anthropogenic pressures. The framework demonstrates novelty by quantifying both factor-driven variability and source contributions, thereby bridging the gap between geostatistical modelling and pollution source identification. This dual perspective enhances the capacity of policymakers to design interventions that address both statistical associations (e.g., soil fertility, texture) and external inputs (e.g., traffic, industry, agriculture) (16, 17). Soil fertility indicators and texture showed strong statistical associations with PHE distribution, suggesting potential explanatory correlates rather than proven causal mechanisms.

In summary, the results position this study within the global sustainability agenda by providing a replicable workflow for identifying contamination hotspots, quantifying source contributions, and prioritizing management interventions. The findings are directly relevant to regional stakeholders in Ningxia Province and contribute to international efforts to safeguard soil resources under the combined pressures of agricultural intensification, industrialization, and climate variability (54–57).

Methodological advances and limitations

This study advances methodological approaches to soil contamination assessment by integrating GIS-based spatial modelling, OPGD explanatory analysis, and PMF source apportionment into a unified framework. The novelty lies in the ability to simultaneously quantify factor-driven heterogeneity and source contributions, bridging the gap between geostatistical modelling and pollution source identification (18, 19).

Despite these advances, several limitations must be acknowledged. First, the dataset was constrained by the availability of soil samples and ancillary variables. Socio-economic indicators (e.g., land management practices, irrigation intensity) were not included, which may have introduced unexplained variability in contamination patterns (21). Second, while PMF effectively resolved five source factors, the interpretation of mixed anthropogenic inputs (Factor 4) remains partly uncertain due to overlapping signatures of agricultural and industrial activities. Third, the OPGD model, although robust, is sensitive to discretization thresholds, and while optimization reduced uncertainty, residual classification bias cannot be fully excluded (22).

Importantly, the hazard mapping approach relied on exceedances of background values and q-value thresholds, which primarily capture spatial contamination intensity but do not directly reflect metal bioavailability or ecological toxicity. For example, Cd hotspots identified in agricultural soils may pose disproportionate risks due to Cd’s high mobility and plant uptake potential, even when concentrations remain below national screening thresholds (23). Similarly, Pb contamination near traffic corridors may have limited ecological impact but significant implications for human health exposure pathways (e.g., inhalation of dust, ingestion of contaminated crops). Without explicit incorporation of bioavailability assays or ecological risk indices, hazard maps may underestimate or misrepresent the true risk to ecosystems and human populations (24, 25).

Future research should therefore integrate bioavailability assessments (e.g., DTPA-extractable fractions, sequential extraction procedures), ecotoxicological endpoints (e.g., soil enzyme activity, microbial diversity), and human health risk models (e.g., hazard quotient, lifetime cancer risk) into the GIS–OPGD–PMF framework. Such coupling would enable hazard maps to move beyond contamination intensity toward risk-based assessments, providing more accurate guidance for remediation and policy interventions (26, 27).

Policy and management recommendations

The identification of contamination hotspots and source contributions provides a strong basis for suggesting targeted policy and management interventions. Translating these scientific insights into practice requires a multi-tiered approach that addresses both local agricultural practices and broader industrial and urban governance frameworks (51–61). The following measures are presented as prospective options informed by model-derived associations, rather than interventions directly tested in this study:

Agricultural management

Given the dominance of Cd enrichment linked to phosphate fertilizers (Factor 5, Figure 4), prospective measures such as adopting low-Cd fertilizers, integrated nutrient management, and organic amendments may help reduce Cd inputs (24, 25). Extension services could promote precision fertilization and soil testing to minimize over-application. Crop rotation with Cd-tolerant species may also reduce bioaccumulation risks in food chains (26). These strategies are suggested as possible management directions, not proven outcomes of this dataset.

Industrial emission control

Zn and Cu hotspots near industrial facilities (Factor 2, Figure 4) suggest that stricter emission standards and improved waste treatment may mitigate localized contamination (27). Enforcement of industrial discharge regulations, coupled with incentives for cleaner production technologies, could be considered. Establishing buffer zones between industrial sites and agricultural land is a potential option to reduce exposure risks.

Traffic-related interventions

Pb contamination (Factor 1, Figure 4) remains linked to traffic emissions and non-exhaust sources. Prospective strategies include encouraging the use of low-wear brake and tire materials, expanding green belts along major roads, and strengthening urban dust management programs (28). These measures align with broader air quality and urban sustainability initiatives but were not tested in this study.

Soil fertility and texture management

OPGD results (Figure 5) indicated that total nitrogen, soil organic matter, and texture components were statistically associated with metal variability. Suggested measures include integrating soil fertility enhancement (e.g., organic amendments, cover cropping) with texture-specific management practices to improve adsorption and reduce mobility of metals (29). These are prospective options informed by associations, not experimentally validated mechanisms.

Risk-based monitoring and remediation

Risk detector analysis (Figure 7) identified discrete hotspots requiring spatially targeted interventions. Monitoring networks could be prioritized in high-risk zones, with remediation strategies tailored to local contamination profiles. Techniques such as phytoremediation, biochar application, and soil washing are highlighted as potential remediation options depending on metal type and soil conditions (30, 31). Their effectiveness requires further validation.

Governance and policy integration

At the institutional level, soil contamination management should be embedded within regional sustainability frameworks, aligning with SDG 11 (Sustainable Cities) and SDG 15 (Life on Land). Cross-sector collaboration among agricultural agencies, environmental regulators, and industrial stakeholders is recommended to ensure coherent policy implementation (32). Establishing open-data platforms for soil monitoring may enhance transparency and community engagement.

Global relevance and future research directions

The integrated GIS–OPGD–PMF framework developed in this study holds significant global relevance, particularly for irrigated agricultural systems facing combined pressures of industrialization, fertilizer dependence, and climate variability. By simultaneously quantifying environmental drivers of metal variability and source contributions, the framework provides a replicable model that can be adapted to diverse agroecosystems worldwide.

Globally, heavy metal contamination in soils has been reported across Asia, Africa, and Latin America, often linked to fertilizer use, industrial emissions, and traffic residues (33, 34). The methodological advances presented here especially the integration of OPGD explanatory power with PMF source apportionment offer a transferable workflow for identifying contamination hotspots and prioritizing remediation. This dual perspective is particularly valuable in regions where monitoring resources are limited, as it enables risk-based targeting rather than uniform interventions (35).

The framework also contributes to international sustainability agendas. By aligning with SDG 2 (Zero Hunger) and SDG 15 (Life on Land), it supports efforts to safeguard soil fertility and food security under conditions of environmental stress. The identification of Cd hotspots linked to fertilizer practices and Pb contamination from traffic emissions demonstrates how local findings can inform global debates on sustainable agriculture and urbanization (36). Moreover, the emphasis on interaction effects (e.g., clay ∩ SOM for As, elevation ∩ SOM for Cr) highlights the importance of coupled soil–environment processes, which are increasingly recognized in global soil health initiatives (37).

Future research should focus on three key directions:

Cross-regional application. Testing the GIS–OPGD–PMF framework in other irrigated systems (e.g., Indo-Gangetic Plain, Nile Delta, Murray–Darling Basin) to evaluate transferability and refine parameterization.

Integration with remote sensing and AI. Incorporating multi-sensor satellite data and machine learning algorithms to enhance spatial resolution and predictive capacity (38).

Socio-economic coupling. Linking contamination assessments with farmer practices, policy frameworks, and cost–benefit analyses to ensure actionable outcomes (39).

By advancing methodological rigor and demonstrating practical applicability, this study positions the GIS–OPGD–PMF framework as a globally replicable model for soil contamination assessment and management. Its adoption across diverse agroecosystems will strengthen resilience, promote sustainable intensification, and contribute to global food security under accelerating environmental change (55–58).

This study developed and applied an integrated GIS–OPGD–PMF framework to assess the spatial distribution, sources, and risks of potentially hazardous elements (PHEs) in irrigated soils of the Yellow River Basin. The results revealed heterogeneous contamination patterns shaped by both anthropogenic activities and intrinsic soil properties, with Cd, Pb, and Zn exhibiting localized hotspots, while Cr and Ni reflected geogenic background influences.

The PMF model resolved five distinct source factors, identifying traffic emissions (Pb), industrial activities (Zn, Cu), geogenic parent material (Cr, Ni), mixed anthropogenic inputs (As), and agriculture-related practices (Cd) as the dominant contributors (Figure 4). The OPGD factor detector confirmed that soil fertility indicators (total nitrogen, SOM) and soil texture components (sand, silt, clay) were the primary drivers of metal variability (Figure 5), while interaction analysis highlighted synergistic effects between soil properties and environmental conditions (Figure 6). Risk mapping identified discrete contamination hotspots in the central–western and southeastern irrigation zones (Figure 7), underscoring the need for spatially targeted monitoring and remediation.

Methodologically, the integration of GIS, OPGD, and PMF represents a significant advance, enabling simultaneous quantification of factor-driven heterogeneity and source apportionment. This dual perspective strengthens the reliability of contamination assessments and provides a replicable workflow for other irrigated agricultural systems worldwide (33–35).

From a policy standpoint, the findings emphasize the urgency of reducing Cd inputs through fertilizer management, tightening industrial emission controls, and mitigating traffic-related Pb contamination. Embedding these strategies within regional sustainability frameworks will contribute directly to achieving SDG 2 (Zero Hunger), SDG 6 (Clean Water and Sanitation), and SDG 15 (Life on Land) (13, 32).

In conclusion, this study demonstrates that soil contamination in irrigated landscapes arises from coupled soil–environment processes and diverse anthropogenic sources. The integrated GIS–OPGD–PMF framework provides a robust scientific basis for risk-based monitoring, targeted remediation, and sustainable soil management, offering both regional relevance and global applicability in safeguarding agroecosystem resilience under accelerating environmental change.