Legacy of non-ferrous metallurgy driving heavy metal-induced geochemical transformation of urban soils

Urban soils in long-established mining and metallurgical cities represent persistent reservoirs of potentially toxic elements (PTEs), posing long-term ecological and human health risks. Ust-Kamenogorsk (Oskemen), one of Kazakhstan’s major non-ferrous metallurgical centers, has experienced more than seven decades of intensive ore processing and smelting. Although previous studies reported elevated concentrations of Cd, Pb, Zn, and Cu, a comprehensive assessment integrating contamination status, ecological risk, spatial variability, and environmental controls has remained limited. Surface soils from 29 sampling sites representing industrial, urban, and background functional zones were analyzed for Zn, Cu, Cd, and Pb. Contamination was evaluated using the contamination factor (CF), geoaccumulation index (Igeo), pollution index (PI), Nemerow pollution index (PN), ecological risk factor (Er), and potential ecological risk index (RI). Spatial patterns were mapped using inverse distance weighting (IDW), while Spearman correlation analysis, principal component analysis (PCA), and non-metric multidimensional scaling (NMDS) were used to evaluate relationships among heavy metals, soil properties, and contamination patterns. Heavy metal concentrations ranged from 5.2-1675.0 mg kg - ¹ (Zn), 0.3-57.8 mg kg - ¹ (Cu), 0.8-39.0 mg kg - ¹ (Cd), and 0.7-59.5 mg kg - ¹ (Pb). Industrial soils exhibited pronounced polymetallic contamination, with Zn dominating overall enrichment and Cd contributing most to ecological risk. The maximum RI reached 2188.9, while the mean RI (674.9) indicated a generally very high ecological risk. Spatial analyses identified contamination hotspots surrounding major metallurgical facilities. Significant correlations among Zn, Cd, and Pb (ρ = 0.68-0.90, p < 0.001), together with PCA and NMDS, were consistent with a strong anthropogenic influence, whereas pH, organic carbon, carbonate content, and soil texture significantly influenced metal accumulation. The integrated assessment demonstrates that combining pollution indices, ecological risk metrics, spatial modeling, and multivariate statistics provides a robust framework for identifying contamination hotspots and environmental drivers of legacy pollution. These findings provide a scientific basis for environmental monitoring, risk-based remediation, and sustainable urban land management in long-established metallurgical cities.

Soils constitute one of the most persistent environmental reservoirs of potentially toxic elements (PTEs), integrating contaminant inputs over extended periods and reflecting the cumulative impacts of anthropogenic activities (1, 2). Unlike degradable organic contaminants, heavy metals persist in soils and occur in chemical fractions with markedly different mobility and bioavailability (3). Their retention, release, and redistribution are controlled by soil properties, particularly pH, organic matter, carbonate minerals, and particle-size composition (4–6). Their long-term accumulation degrades soil quality, disrupts microbial communities and biogeochemical processes, reduces ecosystem functioning, and increases the risk of human exposure through direct contact, inhalation of contaminated dust, and transfer into agricultural products (7–9). Consequently, urban soils have become a critical focus of environmental monitoring because they record both historical and ongoing industrial emissions while serving as an important interface between environmental quality and public health.

Mining and metallurgical cities commonly experience simultaneous inputs from ore processing, smelting, industrial production, traffic, and atmospheric deposition, resulting in complex and persistent contamination patterns (1, 10). Non-ferrous metallurgy, ore processing, smelting, fossil-fuel combustion, traffic emissions, and atmospheric deposition collectively generate complex contamination patterns that persist long after emission rates decline (11). As a result, urban soils frequently exhibit pronounced spatial heterogeneity controlled by emission intensity, prevailing wind direction, land use, topography, and soil characteristics. Understanding these spatial patterns is essential for identifying contamination hotspots, evaluating ecological risks, and supporting evidence-based environmental management in industrial regions.

Kazakhstan hosts several major mining and metallurgical centers where long-term exploitation of polymetallic deposits has resulted in widespread soil contamination. Previous investigations in Ust-Kamenogorsk (Oskemen) demonstrated elevated concentrations of Cd, Pb, Zn, and Cu and provided valuable information on metal fractionation, mobility, and bioavailability (12). Nevertheless, important knowledge gaps remain. Most available studies focused on individual aspects of contamination and were based on relatively limited spatial coverage, whereas integrated assessments combining contamination indices, ecological risk evaluation, spatial distribution, soil physicochemical controls, and multivariate statistical analyses remain scarce. Furthermore, updated assessments are needed because contamination patterns in industrial cities evolve under continuing anthropogenic pressures, land-use changes, and environmental management practices.

Ust-Kamenogorsk represents one of the largest non-ferrous metallurgical centers in Central Asia and has experienced more than seven decades of intensive mining and metallurgical activity (12, 13). The coexistence of zinc and lead smelting, titanium-magnesium production, uranium-related industries, transport infrastructure, and densely populated residential districts has created a complex legacy of polymetallic contamination. These characteristics make the city an important natural laboratory for investigating the interactions between industrial emissions, soil properties, and heavy-metal accumulation, while providing insights applicable to other long-term industrial regions with similar environmental histories (14).

Previous investigations conducted in Oskemen have reported elevated concentrations of Zn, Cd, Pb, and Cu and have provided valuable insights into metal fractionation, mobility, and bioavailability in contaminated soils (12). Nevertheless, despite the considerable body of research devoted to environmental contamination in the city, several important knowledge gaps remain. Most previous studies have focused on individual heavy metals, specific contamination indicators, or relatively limited areas of the urban environment. Comprehensive investigations integrating spatial distribution analysis, quantitative pollution assessment, ecological risk evaluation, soil physicochemical controls, and multivariate statistical analyses across different urban functional zones remain comparatively scarce. Furthermore, the spatial variability of heavy metals in relation to land use, proximity to industrial emission sources, and the combined influence of anthropogenic activities and soil properties has not yet been comprehensively evaluated. Updated assessments are also required because contamination patterns continue to evolve under changing industrial activities, urban expansion, and environmental management practices.

To address these knowledge gaps, the present study applied an integrated framework to evaluate heavy-metal contamination in urban soils of Ust-Kamenogorsk. Geochemical analyses were combined with multiple pollution indices, ecological risk assessment, spatial distribution mapping, correlation analysis, and multivariate statistical techniques to quantify Zn, Cu, Cd, and Pb concentrations, characterize contamination patterns, assess ecological risks, and identify the principal environmental factors governing heavy-metal accumulation and redistribution. By integrating these complementary approaches, this study provides a comprehensive assessment of heavy-metal contamination across contrasting urban functional zones and advances current understanding of the spatial behaviour of heavy metals in a long-established metallurgical city. The findings provide a robust scientific basis for environmental monitoring, ecological risk assessment, and the development of sustainable soil-management strategies for industrial urban environments subjected to prolonged anthropogenic pressure.

Study area

The study was conducted in Ust-Kamenogorsk (Oskemen), East Kazakhstan Region, Republic of Kazakhstan (Figure 1). Oskemen (49°57′00″ N, 82°37′00″ E) has been one of the major centers of non-ferrous metallurgy in Central Asia since the mid-twentieth century, with large-scale production of Pb, Zn, Cu, Ti, and rare metals. More than six decades of metallurgical activities have resulted in persistent atmospheric deposition of potentially toxic elements and the formation of technogenic geochemical anomalies within urban ecosystems.

The investigated area included urban and industrial functional zones located within the city and a local reference site situated approximately 20 km southwest of the urban center. Soil sampling sites were distributed across various geomorphological positions, including river terraces, gentle slopes, interfluves, and foothill plains. The inclusion of these geomorphological units was motivated by their potentially different capacities for sediment accumulation and secondary redistribution of contaminated particles and by the need to capture spatial heterogeneity in contaminant distribution patterns.

The dominant soil types according to WRB (2022) were fluvic kastanozems and haplic kastanozems in urban and background areas, whereas technosols predominated within industrial zones. The predominance of technosols reflects long-term anthropogenic transformation of the soil cover caused by industrial construction, waste deposition, and continuous emissions from metallurgical enterprises. Particularly intensive industrial activities are concentrated within the northeastern and central parts of the city, where metallurgical facilities are located in close proximity to residential districts and major transportation corridors. Detailed environmental characteristics of the sampling sites are presented in Supplementary Table 1.

Oskemen hosts several major industrial enterprises, including the Kazzinc metallurgical complex, the Ust-Kamenogorsk Titanium and Magnesium Plant (UKTMP), and the Ulba Metallurgical Plant (UMP). Kazzinc is one of the world’s leading producers of zinc, lead, copper, and precious metals and represents the principal source of atmospheric emissions of potentially toxic elements in the city. The UKTMP specializes in the production of titanium sponge and magnesium products, whereas the UMP is a major producer of uranium fuel components, beryllium products, and tantalum materials. Long-term operation of these industrial facilities, together with associated waste disposal sites and transportation infrastructure, has resulted in significant accumulation of heavy metals and other pollutants in urban environmental compartments.

The climate of the region is sharply continental, characterized by cold winters, warm summers, and pronounced seasonal temperature fluctuations. The mean annual precipitation is approximately 400–450 mm, most of which occurs during the warm season. The relatively low annual precipitation together with strong seasonal contrasts may favor long-term accumulation of atmospheric fallout in soils, whereas intense precipitation events during the warm season may promote episodic redistribution of contaminated particles through surface runoff.

Several industrial sampling sites were located within relatively low geomorphological positions adjacent to the Irtysh River valley, which may facilitate accumulation of fine technogenic sediments and potentially contribute to elevated concentrations of heavy metals. Topographic analysis based on the Copernicus Digital Elevation Model (DEM) GLO-30 dataset demonstrated considerable relief heterogeneity, with elevations ranging approximately from 320 to 750 m above sea level (Figure 1). Such topographic variability may strongly influence hydrological processes, sediment transport, and localized accumulation of contaminated particles. River terraces and lower geomorphological positions may function as depositional environments for fine technogenic sediments enriched with heavy metals, whereas elevated and exposed areas are more susceptible to erosion and secondary contaminant redistribution.

The coexistence of intensive metallurgical emissions, heterogeneous relief conditions, and contrasting functional zones makes Oskemen a particularly suitable case study for disentangling the relative contributions of industrial activities and environmental factors to heavy metal accumulation in urban soils. In addition, prevailing northwestern and southeastern winds may further influence the directional transport and deposition of atmospheric pollutants, potentially contributing to the formation of spatial contamination gradients across the city. Site S29, located approximately 20 km southwest of Oskemen and outside the principal industrial influence zone, was selected as a local reference site for comparison with the urban and industrial sampling sites. The site was chosen based on its distance from the major industrial facilities and the absence of direct local industrial sources. Because only one reference site was included in the sampling design, S29 was used as a local comparative reference rather than as a representative estimate of the regional natural geochemical background.

Soil sampling and physicochemical characterization

Soil sampling was conducted during the summer season (June-August) in urban and industrial zones of Ust-Kamenogorsk (Oskemen), Kazakhstan. A total of 29 sampling sites were investigated, including urban (n = 12), industrial (n = 16), and local reference site (n = 1) locations (Supplementary Table 1). Sampling sites were selected based on proximity to major metallurgical enterprises, prevailing wind directions, geomorphological conditions, and potential pollutant accumulation zones identified during field reconnaissance and GIS analysis.

Composite soil samples were collected from the upper 0–30 cm horizon using the envelope (diagonal composite) method. Five subsamples were collected at each site and combined into a representative composite sample. At each sampling site, five field subsamples were combined to form one representative composite soil sample. Each composite sample was subsequently analyzed in three analytical replicates (n = 3) for heavy metal concentrations. This depth interval was selected because it represents the soil layer most strongly influenced by long-term atmospheric deposition and anthropogenic disturbances.

Geographic coordinates were recorded using a handheld GPS (GPSMAP 66S, Garmin) receiver. Mesotopography, functional zones, and soil classification according to WRB (2022) were determined based on field observations and physicochemical soil properties. Following transportation to the laboratory, soil samples were air-dried, homogenized, and sieved through a 2-mm mesh prior to analysis.

Soil pH was measured potentiometrically using a laboratory pH meter (SevenCompact S220, Mettler Toledo, Greifensee, Switzerland) in a 1:2.5 soil-to-1 M KCl suspension. Total organic carbon (TOC) content was determined by the dichromate oxidation method followed by titrimetric quantification (15). The analysis was based on oxidation of organic carbon using potassium dichromate (K 2 Cr 2 O 7) in concentrated sulfuric acid (H 2 SO 4), with the excess dichromate determined by titration with ferrous sulfate solution using standard laboratory titration equipment (automatic burette system, Brand GmbH,Wertheim, Germany) (16). Total nitrogen (N) was measured using the Kjeldahl digestion method (17). Available phosphorus was determined with a spectrophotometer by extracting 0.5 M NaHCO 3 (pH 8.5) (18). Exchangeable potassium (K), calcium (Ca), and magnesium (Mg) were extracted with ammonium acetate (AAc) (1 M, pH 7.00).

Sample digestion and heavy metal determination

Heavy metal concentrations in soils were determined following a modified nitric acid-hydrogen peroxide digestion procedure commonly applied in environmental geochemical studies. Approximately 0.5 g of air-dried and homogenized soil sample was transferred into a digestion vessel and treated with 5 mL of concentrated nitric acid (65%, analytical grade). The mixture was pre-digested at room temperature for 30 min and subsequently heated on a digestion block at 120 °C until near dryness. After cooling, 2 mL of hydrogen peroxide (30%) was added to ensure complete oxidation of organic matter. The resulting digest was diluted to a final volume of 25 mL with ultrapure water produced by a Milli-Q purification system (Merck Millipore, Darmstadt, Germany), filtered through quantitative filter paper, and subsequently analyzed.

Concentrations of Pb, Zn, Cd, and Cu were determined by flame atomic absorption spectrometry (FAAS) using an Agilent 240FS atomic absorption spectrometer (Agilent Technologies, Santa Clara, CA, USA) equipped with single-element hollow cathode lamps and an air-acetylene flame system. Analytical measurements were performed at wavelengths of 283.3 nm for Pb, 213.9 nm for Zn, 228.8 nm for Cd, and 324.8 nm for Cu. The gas flow rate was maintained at 50 dm³ h - ¹ with an aspiration rate of 5 cm³ min - ¹.

The investigated metals were selected because Zn, Pb, Cd, and Cu represent the principal pollutants associated with long-term non-ferrous metallurgical activities in Oskemen and have repeatedly been identified as priority contaminants in previous environmental investigations of the region. Their accumulation in urban soils is primarily related to atmospheric emissions from the Kazzinc metallurgical complex, the Ust-Kamenogorsk Titanium and Magnesium Plant, the Ulba Metallurgical Plant, as well as associated transportation and historical waste disposal activities.

Assessment of contamination levels was performed using maximum permissible concentrations (MPCs) commonly adopted in Kazakhstan and neighboring post-Soviet countries. The applied MPC values (mg·kg - ¹) were 32.0 for Pb, 23.0 for Zn, 2.0 for Cd, and 3.0 for Cu (19). These threshold values were subsequently used for calculating contamination factors, integrated pollution indices, and ecological risk parameters.

Method validation and quality assurance

Quantitative determination of heavy metals was performed using an external calibration approach. Calibration curves for Pb, Zn, Cd, and Cu were constructed using certified atomic absorption stock solutions (1.0 g·L - ¹; Agilent Technologies, Santa Clara, CA, USA) through successive dilution with ultrapure water. Working calibration standards covered the concentration range of 1.0-10.0 mg·L - ¹, which adequately encompassed the concentrations observed in soil extracts from urban and industrial environments of Oskemen.

Method validation included assessment of calibration linearity, analytical precision, recovery, and method sensitivity. Each calibration curve was established using four concentration levels, with each level analyzed in five independent replicates. Excellent linear relationships between absorbance and analyte concentration were obtained for all investigated metals, with coefficients of determination (R²) exceeding 0.9995.

Quality assurance and quality control (QA/QC) procedures included procedural blanks, duplicate analyses, calibration verification standards, and recovery tests using certified reference standard solutions. The measured concentrations showed good agreement with certified values at the 95% confidence level, confirming the accuracy of the analytical procedure.

Mean recovery values ranged from 99.1% to 100.3% for Pb, Zn, Cd, and Cu, while analytical precision evaluated through duplicate analyses showed relative standard deviation (RSD) values below 5% for all investigated elements, indicating satisfactory repeatability and reliability of the method.

Limits of detection (LOD) and limits of quantification (LOQ) were determined using repeated measurements of low-concentration standards (1, 2, 5, and 10 mg·L - ¹). The obtained LOD and LOQ values ranged from 0.04 to 0.08 mg·L - ¹ and from 0.12 to 0.24 mg·L - ¹, respectively. These sensitivity levels were considered sufficient for reliable determination of trace metal concentrations in soils affected by long-term metallurgical emissions and elevated concentrations of potentially toxic elements in Oskemen.

Detailed analytical characteristics, including calibration equations, recovery values, limits of detection, and limits of quantification, are presented in Supplementary Tables 2 and 3.

Pollution indices and ecological risk assessment

To comprehensively evaluate soil contamination and ecological risk associated with heavy metal accumulation, several widely used geochemical and ecological indices were calculated.

The contamination factor (CF) was used to evaluate the degree of contamination for individual metals and was calculated as:

where C n is the measured concentration of the metal in soil and B n represents its geochemical local reference concentration.

The Nemerow Pollution Index was calculated using the following equation (20):

where C F m e a n is the arithmetic mean of the contamination factors for all analyzed heavy metals and was calculated as:

where n is the total number of analyzed heavy metals.

The maximum contamination factor was determined as:

Where C F 1, C F 2, …, C F n represent the contamination factors of each analyzed heavy metal, and C F m a x corresponds to the highest contamination factor at each sampling site.

The PN values were classified as follows: PN ≤ 0.7 indicates clean soil; 0.7 < PN ≤ 1.0 indicates the warning level of pollution; 1 < PN ≤ 2 indicates slight pollution; 2 < PN ≤ 3 indicates moderate pollution; and PN > 3 indicates heavy pollution.

The geoaccumulation index (Igeo) was calculated according to the equation proposed by Müller (21):

where the constant factor 1.5 compensates for natural lithogenic variability and minor background fluctuations (Supplementary Table 4).

The ecological risk factor (Er) was determined as:

where T r denotes the toxic-response coefficient specific to each heavy metal. The toxic-response coefficients used in ecological risk calculations were 30 for Cd, 5 for Cu and Pb, and 1 for Zn according to Hakanson (Supplementary Table 5) (22).

The pollution index (PI) was calculated as the arithmetic mean of the ratios between the measured concentrations of individual heavy metals and their corresponding maximum permissible concentrations (MPCs):

Potential ecological risk index (RI) values were calculated as the cumulative sum of ecological risk factors:

Although the potential ecological risk index (RI) proposed by Hakanson was originally developed for assessing contamination in aquatic sediments, both RI and the ecological risk factor (Er) have been widely adopted for evaluating heavy metal pollution in contaminated soils, including mining- and smelting-affected environments. Their broad application is attributed to the integration of contamination intensity with metal-specific toxic-response coefficients, providing a standardized framework for comparing the potential ecological risks posed by different heavy metals across diverse environmental settings.

Because reliable regional geochemical background values are unavailable for the Oskemen area and truly pristine soils unaffected by long-term mining and metallurgical activities are difficult to identify, the local reference concentrations used to calculate the contamination factor (CF), geoaccumulation index (Igeo), Nemerow Pollution Index (PN), and Potential Ecological Risk Index (RI) were derived from the local reference site (S29). This site is located outside the principal zone of industrial influence and exhibited generally low heavy metal concentrations and pollution index values, making it a suitable local reference point for evaluating anthropogenic enrichment. The concentrations measured at S29 were therefore used as site-specific comparative reference values and should not be interpreted as representative of the regional natural geochemical background.

To achieve a comprehensive assessment of soil contamination, the Nemerow Pollution Index (PN) was interpreted together with the contamination factor (CF), geoaccumulation index (Igeo), Pollution Index (PI) and Potential Ecological Risk Index (RI). While CF and Igeo quantify the enrichment and accumulation of individual heavy metals relative to local reference concentrations, PI evaluates the average level of multi-element contamination relative to the corresponding MPC values, RI estimates the integrated ecological risk by incorporating metal-specific toxic-response coefficients, and PN simultaneously considers both the average contamination level and the maximum contamination intensity. The complementary application of these indices provides an integrated assessment of contamination severity, potential ecological risk, and the spatial characteristics of heavy metal pollution relative to the local reference conditions.

Statistical analysis

Statistical analyses and graphical visualizations were performed in the R statistical environment (version 4.5.1; R Core Team), whereas spatial interpolation and cartographic analyses were conducted using QGIS software (version 3.40.7).

Descriptive statistics, including minimum, maximum, mean, median, standard deviation (SD), and interquartile range (IQR), were calculated to characterize the variability of soil physicochemical properties and heavy metal concentrations among sampling sites and functional zones. Data distributions were visualized using boxplots and violin plots, while soil texture composition was illustrated using USDA ternary diagrams.

Relationships among soil physicochemical properties, heavy metals, and pollution indices were evaluated using Spearman rank correlation analysis. Correlation matrices were constructed using Z-score standardized data and visualized as correlation heatmaps. Statistical significance was adjusted using the false discovery rate (FDR) correction procedure to minimize type I errors associated with multiple comparisons.

Principal component analysis (PCA) was performed to identify the major environmental factors controlling heavy metal distribution and to evaluate potential contamination sources in urban soils. Prior to PCA, quantitative variables were standardized (z-score transformation) to eliminate scale effects among soil properties and metal concentrations.

Spatial distribution patterns of heavy metal contamination were visualized using inverse distance weighting (IDW) interpolation implemented in QGIS. The concentration-to-MPC ratios (Zn/MPC, Cu/MPC, Cd/MPC, and Pb/MPC) were used as interpolation variables. A distance coefficient (power) of 2 and an output grid resolution 30 × 30 m were applied. Spatial analyses were performed in the WGS 84/UTM zone 44N coordinate reference system (EPSG:32644), and the resulting raster surfaces were clipped to the study-area polygon. The predictive performance of the IDW interpolation was evaluated using leave-one-out cross-validation (LOOCV) based on the 28 sampling sites (S1–S28) used for spatial interpolation. The local reference site (S29), which was not included in the original IDW interpolation, was also excluded from the cross-validation. Each of the 28 sites was sequentially omitted, and its value was predicted from the remaining 27 sites using the same distance coefficient (p = 2). Prediction errors were summarized using mean error (ME), mean absolute error (MAE), and root mean square error (RMSE). Statistical significance was accepted at p < 0.05.

Differences between urban and industrial functional zones were assessed using the Mann-Whitney U test, a non-parametric alternative to the independent-samples t -test, because several variables deviated from normality and the sample sizes among functional zones were unequal. To reduce the probability of type I errors associated with multiple comparisons, raw p -values were adjusted using the Benjamini-Hochberg false discovery rate (FDR) procedure. The magnitude of zonal differences was quantified using the rank-biserial correlation coefficient (rrb), which provides an effect-size estimate for Mann-Whitney tests. Following conventional interpretation, absolute rrb values approaching 1 indicate very strong separation between groups, whereas values close to 0 indicate negligible differences. Statistical significance was considered at FDR-adjusted p < 0.05. Violin plots with embedded boxplots and individual sampling points were used to visualize variable distributions and support the interpretation of statistical comparisons between functional zones.

Soil physicochemical properties and texture characteristics

The investigated soils exhibited considerable physicochemical heterogeneity (Supplementary Table 6). Soil pH ranged from 5.73 to 7.48, with a mean value of 6.81, indicating predominantly neutral to slightly alkaline conditions across the study area. Carbonate contents were generally low, varying from 0.61 to 4.52%, whereas TOC showed substantial variability, ranging from 0.42 to 9.88% with a coefficient of variation (CV) of 66.2%. Available phosphorus exhibited the highest variability among the investigated soil properties (CV = 82.7%), followed by exchangeable magnesium and carbonate contents (Table 1).

Descriptive comparison of functional zones showed some variation in soil physicochemical properties (Supplementary Table 7). Industrial soils exhibited higher median available P 2 O 5 concentrations than urban soils (62.0 vs. 46.0 mg·kg - ¹), whereas median exchangeable K 2 O concentrations were identical between the two zones (480 mg·kg - ¹). In contrast, median Ca² + concentrations were lower in industrial soils (106.65 vs. 133.35 mmol(+)/kg). Median pH and CaCO 3 values were similar between the two functional zones.

Particle-size composition also demonstrated pronounced spatial variability. Sand content ranged from 14.1 to 83.0%, whereas silt varied between 13.9 and 71.5%, indicating substantial heterogeneity in sedimentary environments and parent materials (Supplementary Table 8). According to USDA classification, the investigated soils were predominantly represented by sandy loam and silt loam textures, with several sites classified as loam or loamy sand (Supplementary Figure 1). Industrial zones were generally enriched in coarse fractions, exhibiting a median sand content of 65.6%, whereas urban soils contained considerably higher proportions of silt (median = 48.3%). Clay contents remained uniformly low across all functional zones, with median values below 2.5%.

Spatial variability and enrichment of heavy metals in urban soils

Concentrations of heavy metals in Oskemen soils exhibited pronounced spatial heterogeneity and strong enrichment (Supplementary Table 9). Among the investigated elements, Zn was the dominant contaminant, with concentrations ranging from 5.2 to 1675.0 mg·kg - ¹ and an average concentration of 419.6 mg·kg - ¹, approximately 23 times higher than the concentration recorded at the local reference site (S29). Cadmium represented the second most enriched element, varying between 0.8 and 39.0 mg·kg - ¹ with a mean concentration of 10.4 mg·kg - ¹, approximately 13 times higher than the concentration recorded at S29. Copper and lead exhibited mean concentrations of 9.0 and 14.5 mg·kg - ¹, respectively, although their maximum values reached 57.8 mg·kg - ¹ and 59.5 mg·kg - ¹ at individual sites. Coefficients of variation exceeded 110% for all investigated metals, demonstrating substantial spatial variability in their concentrations.

Distinct differences were observed between functional zones. Urban soils generally contained relatively low concentrations of heavy metals, with median concentrations of 24.1, 1.15, 1.55, and 1.65 mg·kg - ¹ for Zn, Cu, Cd, and Pb, respectively. In contrast, industrial soils exhibited substantially elevated concentrations, with median values increasing to 704.5 mg·kg - ¹ for Zn, 11.5 mg·kg - ¹ for Cu, 15.8 mg·kg - ¹ for Cd, and 22.1 mg·kg - ¹ for Pb. These values indicate approximately 29-, 10-, 10-, and 13-fold enrichment of Zn, Cu, Cd, and Pb, respectively, in industrial areas compared with urban environments.

The highest concentrations of Zn, Cd, and Pb were predominantly recorded at sites S16-S25, located within the central industrial sector of the city and in close proximity to major metallurgical enterprises and transportation corridors. Site S19 exhibited the maximum concentrations of Zn (1675.0 mg·kg - ¹) and Pb (59.5 mg·kg - ¹), whereas the highest Cd and Cu concentrations were observed at sites S20 (39.0 mg·kg - ¹) and S25 (57.8 mg·kg - ¹), respectively. Conversely, the local reference site (S29), located outside the principal industrial influence zone, exhibited low concentrations of all investigated metals, although Zn concentrations at several urban sites (S3, S5, S6, and S9) were lower than that recorded at S29.

Comparison with maximum permissible concentrations (MPCs) further emphasized the severity of contamination. Zinc concentrations exceeded the MPC (23 mg·kg - ¹) at numerous industrial sites by up to approximately 73 times, while Cd exceeded its threshold value (2 mg·kg - ¹) by nearly 20 times. Copper concentrations reached up to 19 times the MPC, whereas Pb exceedances were comparatively lower, reaching approximately two-fold above the regulatory threshold at the most contaminated sites (Supplementary Figure 2).

The Mann-Whitney U test revealed that most soil physicochemical properties did not differ significantly between urban and industrial functional zones (Figure 2A, Supplementary Table 7). Exchangeable Na + was the only soil property that remained significant after Benjamini-Hochberg correction (FDR-adjusted p = 0.002) and was significantly higher in urban soils. In contrast, industrial soils contained significantly higher concentrations of Zn, Cu, Cd, and Pb than urban soils, while pollution indices (PI, PN, and RI) were also markedly elevated. All differences remained significant after FDR correction (FDR-adjusted p < 0.001) (Figure 2B). The corresponding effect size analysis further confirmed these patterns, with the strongest zonal separation observed for Zn, PN, and PI, each exhibiting a rank-biserial correlation coefficient (rrb) of 1.00, whereas RI also showed a very large effect size (rrb = 0.948) (Figure 2C).

Integrated assessment of soil pollution and ecological risk

The average pollution index (PI) varied from 0.26 to 25.08, while the Nemerow pollution index (PN) ranged between 0.37 and 54.27 (Figure 3A, Supplementary Table 10). Most urban and peripheral sites exhibited PN values below 1, corresponding to clean or warning levels of contamination, whereas several sites located within the industrial belt substantially exceeded the threshold for heavy pollution (PN > 3). The highest PN values were recorded at sites S19 (54.27), S20 (51.77), S17 (40.77), and S25 (38.18) (Figure 3A).

Contamination factor (CF) values varied substantially among individual metals (Figure 3B). Maximum CF values reached 93.06 for Zn, 48.75 for Cd, 192.67 for Cu, and 74.38 for Pb. The highest Cu contamination was observed at site S25.

The geoaccumulation index further confirmed the severity of anthropogenic impacts (Supplementary Figures 3, 4A). Several industrial sites were classified as heavily to extremely contaminated with respect to Zn and Cd, with maximum Igeo values reaching 5.96 and 5.02, respectively. Copper generally fell within moderate to extremely contaminated categories and reached the highest Igeo value among all investigated elements (Igeo = 7.01 at site S25). In contrast, Pb exhibited a more heterogeneous pattern, with contamination largely confined to selected industrial locations (Supplementary Figures 3, 4A).

Ecological risk indices showed even greater spatial contrasts than pollution indices and identified cadmium as the dominant environmental threat in Oskemen soils (Supplementary Figures 4B,C). The integrated potential ecological risk index (RI) ranged from 41.0 to 2188.9, with approximately one-third of the sampling sites exceeding the threshold for very high ecological risk (RI > 600). The highest RI values were recorded at sites S25 (2188.9), S20 (2163.1), S19 (1966.2), and S17 (1627.2), forming a distinct ecological risk belt within the central industrial zone (Supplementary Figure 4B).

Among individual metals, Cd overwhelmingly dominated ecological risk. Ecological risk factors (Er) for Cd ranged from 30.0 to 1462.5, substantially exceeding those of Zn and Pb and often surpassing the threshold for very high ecological risk (Er > 320) (Figure 3C, Supplementary Figure 5). Copper represented the second most important contributor to ecological risk, with Er values reaching 963.3 at site S25, whereas Zn exhibited relatively low ecological risk despite extremely high concentrations, with Er values remaining below 100 at all sites. Lead generally contributed only moderately to total ecological risk and was largely restricted to areas immediately surrounding industrial emission sources (Supplementary Figures 4C, 5).

Environmental drivers and multivariate patterns of heavy metal accumulation

Strong positive correlations were observed among Zn, Cu, Cd, and Pb (ρ = 0.51-0.90, p < 0.01). Particularly strong associations between Zn and Cd (ρ = 0.90), Zn and Pb (ρ = 0.85), and Cu and Cd (ρ = 0.87) (Figure 4A).

Heavy metals also exhibited significant relationships with several soil properties. Positive correlations between Zn, Cd, and available phosphorus (ρ = 0.38-0.48) and exchangeable potassium (ρ = 0.46-0.85). In contrast, pH was negatively correlated with Zn (ρ = −0.46), Cd (ρ = −0.23), and Pb (ρ = −0.29), indicating that higher total metal concentrations tended to occur at sites with relatively lower pH values. Exchangeable Ca² + similarly showed negative relationships with Zn and Pb.

Particle-size composition exerted a particularly strong influence on metal distribution. Sand content was negatively correlated with Zn (ρ = −0.46), Cu (ρ = −0.58), Cd (ρ = −0.51), and Pb (ρ = −0.57), whereas silt exhibited significant positive relationships with these metals (ρ = 0.42-0.61). The almost perfect negative relationship between sand and silt fractions (ρ = −0.98).

Principal component analysis further clarified the factors controlling contamination patterns. The first two principal components explained 58.9% of the total variance, with PC1 and PC2 accounting for 33.4% and 25.5%, respectively (Figure 4B). PC1 was strongly associated with Zn, Cd, Cu, Pb, and sand content and represented a distinct anthropogenic contamination gradient separating industrial and urban sites. Sampling locations characterized by high ecological risk (S13-S26) clustered on the negative side of PC1 and were closely associated with elevated heavy metal concentrations and pollution indices. Conversely, low-risk urban sites (S1-S12 and S29) were primarily associated with higher pH values, finer soil fractions, and elevated concentrations of exchangeable base cations.

PC2 reflected the influence of soil physicochemical properties, particularly TOC, nutrient status, and exchangeable cations. Sites S2, S4, S5, and S12 were positively associated with silt content, Ca² +, K +, and nutrient enrichment, whereas sites characterized by elevated sand content and lower nutrient availability occupied the lower portion of the ordination space.

Non-metric multidimensional scaling (NMDS) further supported these findings (Figure 4C). The ordination demonstrated excellent representation of the dataset (stress = 0.028), while PERMANOVA analysis confirmed highly significant differences among ecological risk groups (R² = 0.841, p = 0.001). Pollution indices (RI, PI, and PN) were strongly aligned with the industrial cluster, whereas pH, CaCO 3, and silt content were associated with low-risk urban soils.

Spatial distribution patterns and identification of contamination hotspots

Spatial interpolation revealed pronounced heterogeneity in the distribution of heavy metals across Oskemen and identified several distinct contamination hotspots (Figure 5). Zn, Cu, and Cd showed elevated concentrations predominantly within the central industrial sector, whereas lower concentrations were generally observed in peripheral urban areas.

Zinc displayed the most extensive contamination footprint among all investigated elements (Figure 5A). Elevated Zn concentrations occupied a large portion of the central and southwestern sectors of the study area, with contamination levels exceeding 20–60 MPCs and locally surpassing 60 MPCs in several industrial hotspots. The most pronounced Zn anomalies were concentrated around sites S16-S20 and S24-S26.

Copper displayed a more localized distribution pattern (Figure 5B). Although elevated Cu concentrations were also concentrated within the industrial sector, contamination hotspots were spatially less extensive and primarily associated with sites S20, S24, and S25. Maximum contamination levels exceeded 15 MPCs, but the area affected by severe Cu contamination was considerably smaller than that observed for Zn and Cd.

Cadmium exhibited a spatial distribution similar to that of Zn (Figure 5C). Extremely high Cd contamination (>15 MPCs) was concentrated within the central industrial zone, particularly near sites S16-S20 and S25.

Lead exhibited the least extensive contamination among the investigated metals, although elevated Pb concentrations were still evident across parts of the central industrial sector (Figure 5D). The highest Pb/MPC ratios occurred near sites S16, S17, S19, S24, and S25, whereas most urban and peripheral areas remained below or only slightly above the MPC threshold.

The interpolation maps showed a clear spatial contrast between the central industrial sector and the northeastern part of Oskemen, where comparatively low concentrations of Zn, Cu, Cd, and Pb were observed. Overall, the spatial analysis showed broader distributions of elevated Zn and Cd concentrations, whereas Cu and particularly Pb exhibited more localized contamination patterns.

Spatial analysis demonstrates that Oskemen is characterized by a well-defined industrial contamination gradient, with Zn and Cd representing the most widespread pollutants, whereas Cu and particularly Pb exhibit more localized distribution patterns. Leave-one-out cross-validation indicated variable predictive performance of the IDW interpolation, with R² values ranging from 0.187 to 0.467. Detailed validation statistics, including mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) for each metal, are provided in Supplementary Table 11.

The spatial distribution of heavy metals in Oskemen is consistent with a substantial long-term legacy of non-ferrous metallurgical activity rather than contemporary emissions alone. For more than six decades, ore mining and metallurgical processing have released substantial quantities of Zn, Cd, Cu and Pb into the urban environment, resulting in persistent technogenic geochemical anomalies that remain evident despite improvements in industrial emission control. The agreement among chemical analyses, pollution indices, multivariate statistics, and spatial interpolation indicates that these anomalies are spatially coherent and is consistent with a substantial long-term industrial legacy.Similar contamination patterns have been described in historical metallurgical centres worldwide, where urban soils continue to function as reservoirs of heavy metals long after reductions in industrial emissions because metals are efficiently retained within the upper soil horizons and only slowly removed through natural processes (23–25).

The contamination profile observed in Oskemen was dominated by Zn, followed by Cd, Cu and Pb, although each element exhibited distinct spatial behaviour reflecting differences in industrial sources and environmental transport. Zinc showed the highest concentrations and the broadest spatial distribution, reaching 1675 mg kg - ¹, approximately twenty-three times the local reference concentration. Such enrichment is characteristic of regions affected by polymetallic ore processing because Zn is one of the principal components of lead-zinc ores and is emitted in large quantities during roasting, smelting and refining. Its relatively fine particle size facilitates atmospheric transport over considerable distances before deposition, explaining the extensive contamination belt identified around the industrial zone. Comparable Zn enrichment has been documented in long-established smelting districts of China, Poland and Spain, where atmospheric deposition remains the dominant pathway of soil contamination even after improvements in emission-control technologies (26–28).

Although cadmium occurred at considerably lower absolute concentrations than Zn, it represents the most environmentally significant contaminant identified in the present study. Cadmium enrichment exceeded thirteen times the local reference concentration and showed an exceptionally strong correlation with Zn, suggesting that both elements may share common or closely related sources.This relationship is expected because Cd commonly substitutes for Zn within sphalerite and is released simultaneously during metallurgical processing of sulphide ores (29). The almost identical spatial distribution of Zn and Cd observed in the interpolation maps provides additional evidence that these metals originate predominantly from the same emission sources rather than from independent anthropogenic activities. Such coupled Zn-Cd enrichment has repeatedly been recognised as a characteristic geochemical signature of non-ferrous metallurgy and has been reported in numerous mining and smelting regions throughout Europe and Asia (30, 31). Consequently, the contamination pattern observed in Oskemen is consistent with global observations from historical polymetallic mining districts.

Copper and lead displayed markedly different spatial behaviour. In contrast to Zn, elevated Cu concentrations were confined to relatively small areas surrounding specific industrial facilities, while Pb exhibited the most localized distribution among the investigated metals. These differences probably reflect variations in the physical properties of emitted particles. Coarser Cu- and Pb-bearing particles generally settle rapidly after emission, limiting their atmospheric transport, whereas finer Zn- and Cd-rich aerosols remain suspended for longer periods and are dispersed over larger areas (5, 32). The exceptionally high contamination factor recorded for Cu at individual sampling sites therefore most likely reflects localized industrial emissions rather than regional atmospheric transport. Likewise, the restriction of Pb anomalies to the vicinity of industrial enterprises supports previous observations that lead contamination decreases rapidly with increasing distance from emission sources because of its comparatively low atmospheric mobility (33, 34).

The comparison between functional zones provides further evidence of a strong anthropogenic influence on heavy metal accumulation in Oskemen. Median concentrations of Zn, Cu, Cd, and Pb in industrial soils exceeded those measured in urban soils by approximately one order of magnitude, while the local reference site (S29) generally exhibited low metal concentrations relative to the industrial sites. Such pronounced differences are unlikely to reflect natural variability in parent material or soil development alone and are consistent with a substantial contribution from long-term metallurgical activities. Similar industrial-urban gradients have been reported in numerous historical smelting cities, where contamination intensity decreases progressively away from production facilities (35–37).

Although industrial emissions clearly determine the magnitude of contamination, the environmental behaviour of heavy metals after deposition is strongly influenced by soil physicochemical properties. The soils investigated in the present study were predominantly neutral to slightly alkaline, conditions generally considered favourable for reducing metal mobility through adsorption onto clay minerals, Fe-Mn oxides and soil organic matter (38, 39). Nevertheless, despite these relatively favourable geochemical conditions, industrial sites remained heavily contaminated. This finding demonstrates that soil buffering capacity modifies the mobility of heavy metals but cannot compensate for continuous atmospheric inputs accumulated over several decades. In other words, physicochemical properties regulate contaminant retention rather than their primary source.

The correlation analysis provides further insight into the mechanisms controlling heavy metal distribution within the urban environment. The strong negative relationships between sand content and all investigated metals, together with positive correlations involving the silt fraction, indicate that fine-textured soils serve as efficient sinks for airborne industrial particulates. This observation agrees with well-established principles of soil geochemistry, whereby finer particles possess greater specific surface area, higher cation-exchange capacity and a larger number of reactive sorption sites capable of retaining trace metals (4, 40). Once deposited, metal-rich particulates become preferentially associated with silt- and clay-sized fractions, reducing their susceptibility to rapid removal while simultaneously promoting long-term accumulation in surface soils. This mechanism also explains why the most contaminated areas identified by spatial interpolation coincide with landscapes characterised by relatively fine-textured soils.

Moderate positive relationships between Zn, Cd, available phosphorus and exchangeable potassium probably reflect the co-deposition of metallurgical dust enriched in both trace metals and alkaline mineral particles rather than direct geochemical interactions among these variables. Similar associations have been observed in industrial regions where atmospheric emissions modify not only heavy metal concentrations but also selected soil fertility indicators through continuous deposition of mineral aerosols (4, 5, 41). Conversely, the negative correlations between pH and several metals indicate that sites with relatively lower pH tended to contain higher total metal concentrations. This relationship should not be interpreted as direct evidence that lower pH simultaneously promotes metal accumulation and mobility. In general, decreasing pH can increase the solubility and potential mobility of many metals by reducing adsorption to soil surfaces and increasing competition for sorption sites (42). However, the present study measured total metal concentrations rather than mobile or bioavailable fractions; therefore, the observed correlations primarily describe spatial co-variation between soil pH and total metal concentrations. Given the relatively narrow pH range across the study area, pH is likely to act as a secondary factor modifying metal retention and redistribution after deposition rather than as the primary cause of the contamination gradient.

The integrated application of contamination indices provided a more comprehensive understanding of soil quality than individual metal concentrations alone. Although CF, Igeo, PN and RI quantify different aspects of contamination, these approaches consistently identified the industrial sector of Oskemen as the area of greatest environmental concern. Their agreement provides complementary evidence of the pronounced contamination gradient between industrial and urban soils, while the absolute magnitude of reference-based indices should be interpreted in relation to the local reference concentrations measured at S29. Similar conclusions have been reported by Kowalska et al. (43), who emphasized that combining several pollution indices substantially improves the interpretation of contaminated soils because each index captures different aspects of environmental degradation. The consistency between PN values and the interpolation maps obtained in the present study demonstrates that integrated pollution indices accurately reflect the spatial distribution of heavy metals across Oskemen (12, 44).

The contamination factor and geoaccumulation index provided complementary information on element-specific enrichment. Zinc exhibited the greatest overall enrichment across the study area, whereas copper reached the highest contamination factor at individual industrial sites. Maximum Igeo values classified Cu, Zn and Cd as heavily to extremely contaminated at specific locations. These findings illustrate that average concentrations alone may underestimate localized contamination and support the complementary use of multiple assessment indices.

Perhaps the most important outcome of the ecological risk assessment is the dominant contribution of cadmium. Although Zn was by far the most abundant contaminant, Cd accounted for the overwhelming proportion of total ecological risk throughout the industrial area. This pattern is consistent with the Hakanson framework (22), in which ecological risk depends not only on contamination intensity but also on the toxic-response coefficient assigned to each element. Because cadmium has substantially higher toxic-response coefficient than Zn, Cu and Pb, comparatively lower Cd concentration can generate disproportionately high ecological risk. Its environmental relevance is further supported by its relatively high mobility and potential biological availability under a range of soil conditions (12). Similar dominance of Cd in ecological risk assessments has been reported in other mining and metallurgical regions (30, 45, 46).

The distinct clustering of industrial sites (S13-S26) in the PCA ordination demonstrates that contamination is systematic rather than random, with soils adjacent to metallurgical facilities exhibiting a characteristic geochemical signature defined by elevated heavy metal concentrations and pollution indices (42, 47, 48). Similar ordination patterns have been reported in industrial regions where anthropogenic variables dominate the primary environmental gradient (30, 45). The NMDS analysis provided complementary evidence. Together, the ordination results support the pronounced geochemical differentiation between industrial and urban soils and are consistent with previous studies identifying non-ferrous metallurgy as a major contributor to heavy metal accumulation in Oskemen (12).

Spatial interpolation complemented the multivariate analyses by illustrating how the contamination gradient was expressed across the urban landscape. Zn and Cd showed broader spatial anomalies, whereas Cu and Pb exhibited more localized distributions (31, 49).

Local topography may additionally influence contaminant retention and redistribution; however, these processes were not quantified directly in the present study. Moreover, LOOCV indicated variable predictive performance of the IDW interpolation (R² = 0.187-0.467). The interpolated surfaces should therefore be interpreted primarily as representations of broad spatial contamination patterns rather than as precise predictions at unsampled locations. Future studies with higher sampling density, particularly within spatially heterogeneous industrial areas, would improve the robustness of spatial predictions.

The persistence of elevated metal concentrations within industrial sectors indicates that contaminated urban soils remain long-term reservoirs of heavy metals and potential sources of secondary redistribution through dust resuspension, erosion, and surface runoff (24, 50). Environmental management should therefore consider both current emissions and contaminants already accumulated in surface soils. Importantly, remediation priorities should not be based solely on total metal concentrations: although Zn dominated in concentration and spatial extent, Cd contributed disproportionately to ecological risk. Risk-based monitoring and targeted management of Cd-enriched hotspots may therefore provide a more appropriate basis for remediation planning, while stabilization approaches using mineral amendments, organic sorbents, or phytostabilization may represent suitable long-term management options (9, 51, 52).

Although this study provides a comprehensive assessment of heavy metal contamination, several limitations should be acknowledged. An important limitation is the use of a single local reference site (S29) for deriving the reference concentrations used in the contamination and ecological risk indices. Although S29 is located outside the principal industrial influence zone and provides a consistent local basis for comparison, a single site cannot fully capture the natural spatial variability of soil geochemistry across the Oskemen region. Consequently, the absolute values and classification levels of reference-based indices, particularly CF, Igeo, Er, and consequently RI, may be sensitive to the concentrations measured at S29. The resulting indices should therefore be interpreted primarily as measures of contamination relative to this local reference rather than as estimates based on a regional natural geochemical background. Future studies incorporating multiple reference sites across comparable soil types and geomorphological settings would provide a more robust estimate of regional background variability. The investigation was based on total metal concentrations and ecological risk indices, whereas metal speciation and bioavailability were not determined. Because the mobility and ecological effects of trace metals depend strongly on their chemical forms, future studies should combine total concentration measurements with sequential extraction techniques to better characterize contaminant stability and potential mobility. In addition, incorporating human health risk assessment and source apportionment approaches would provide a more comprehensive evaluation of the environmental consequences of historical metallurgical activity.

Long-term non-ferrous metallurgical activity has contributed to pronounced contamination gradients, localized geochemical anomalies, and areas of elevated ecological risk in the urban soils of Oskemen. The integrated assessment provides strong evidence of anthropogenic influence, with Zn representing the most spatially extensive contaminant and Cd constituting the principal contributor to ecological risk because of its high toxicity and toxic-response coefficient.

Spatial and multivariate analyses consistently showed a pronounced industrial–urban contamination gradient, while soil physicochemical properties and texture may influence the retention and redistribution of deposited metals.

The findings have direct implications for environmental management. Monitoring and remediation efforts should prioritize highly contaminated industrial areas rather than applying uniform management strategies across the entire urban area. In particular, ecological risk assessments should account not only total metal concentrations but also for differences in metal toxicity, as cadmium represented the dominant contributor to ecological risk despite the greater abundance of zinc.

The main limitations of this study include the use of a single local reference site, the assessment of total rather than speciated or bioavailable metal concentrations, and uncertainty in spatial predictions at unsampled locations. Future investigations incorporating multiple reference sites, higher spatial sampling density, metal fractionation, source-apportionment approaches, and human health risk assessment would provide a more comprehensive understanding of contamination and strengthen the scientific basis for risk-oriented environmental management.