Mapping Soil Contamination in Arid Regions: A GIS and Multivariate Analysis Approach

Abstract

Heavy metal soil contamination is a global environmental issue that poses serious threats to human health, agricultural advancement, and ecosystem systems. Thirty-five soil samples from various parts of Jazan, Southwest Saudi Arabia, were collected. To create spatial pattern maps for nine potentially toxic elements (PTEs) (As, Co, Cr, Cu, Fe, Ni, Pb, V, and Zn), Ordinary Kriging (OK) was utilized. The variability of the soil metal concentration was estimated using multivariate analysis, including principal component analysis (PCA) and cluster analysis. In addition, the levels of soil contamination in the research area were assessed using contaminations indices, namely, Enrichment Factor (EF), Contamination Factor (CF), and geoaccumulation index (Igeo), and modified contamination degree (mC_d). Normalized Difference Vegetation Index (NDVI) and land use/land cover (LULC) were assessed to evaluate vegetation density and identify different forms of land cover and land use. The results showed that the Gaussian model fitted As well, whereas the spherical model fitted Co, Cr, Cu, Ni, and Zn. An exponential model was fitted to Fe and V. Pb also suited the Stable model. In each of the selected metals, the root mean square standardized error (RMSSE) values were close to one, and the mean standardized error (MSE) values were almost zero for each fitted model. Moreover, the findings showed that there was a tendency for the concentration of heavy metals in the research area to rise from west to east. The cluster analysis divided the data in this investigation into two clusters. Significant alterations in Co, Cr, Cu, Fe, Ni, V, and Zn were revealed by the acquired data. However, the total As and Pb concentrations in the two clusters did not differ significantly. The mC_d value of the research region often fell into one of three classes, with areas of 148.20 km² (nil to very low degree of contamination), 26.16 km² (low degree of contamination), and 0.495 km² (moderate degree of contamination). The findings indicated that the minerals connected to the Arabian Shield’s basement rocks are the main source of these PTEs. It is crucial to monitor PTEs contamination because the research region is highly cultivated, as shown by the NDVI and LULC status. Given the potential for future pollution due to human activity, PTEsPTEs decision-makers may use the findings of the spatial distribution maps of pollutants and their concentrations as a basis for future monitoring of PTEs concentrations in the study area.

1. Introduction

Heavy metal contamination of soil is regarded as one of the major concerns worldwide, due to its direct and indirect effects on human health [1]. Agricultural soil immediately impacts human health and plays a critical role in maintaining food safety [2]. Significant increases in industrial activity, rapid urbanization, and population growth have occurred throughout the past ten years [2]. Heavy metals (PTEs) have spread into all environmental compartments as a result of the large volumes of solid and liquid waste that have been produced. As a consequence, there has been a significant decline in soil and water quality, endangering human health as well as marine life [3,4]. A global environmental concern, soil heavy metal pollution, has garnered significant public attention largely due to growing worries about the safety of agricultural goods [5]. Five million sites of soil pollution exist worldwide, spanning 500 million hectares of land, where various heavy metals or metalloids have contaminated the soils. The cumulative annual economic impact of soil contamination from heavy metals is projected to be more than $10 billion worldwide [6].

Globally, Saudi Arabia’s progress during the past three decades, particularly in the agriculture sector, has been astounding. A significant achievement in a long time is the conversion of vast swaths of desert into agricultural land, even though Saudi Arabia is a nation with one of the lowest annual rainfall rates in the world, receiving an average of just 101.6 mm. Nevertheless, there are several environmental issues associated with this development [7]. PTEs are a class of chemical elements that can be extremely harmful to the environment and human health, even in relatively low concentrations. These elements include, among others, As, Pb, Cd, Hg, and Cr [8,9]. Particularly, even in trace amounts, Cd, Hg, Pb, and Cr have detrimental effects that cause both acute and long-term toxicities in the body [10]. PTEs occur naturally as constituents in the chemical composition of rocks and are thus released through wind and water weathering as well as volcanic eruptions [11]. Conversely, anthropogenic sources of PTEs are especially important and include airborne sources, coal ash, mining of metals, waste disposal, pesticides, inorganic fertilizers, lead-based paints, petrochemicals, and leaded petrol [12,13,14]. Through the inhibition of enzyme activity, competition with necessary cations, and generation of oxidative stress, an excess of PTEs generates abiotic stress for soil biota [15]. As a result, the plant life cycle—from seed germination to maturity stage—is negatively impacted, which ultimately lowers agricultural production and quality [16]. Food, feed, and fodder crops, known as hyperaccumulators, can absorb metals and are extremely resilient to metal stress and excess heavy metals, which then transfer to their aerial sections. A key problem is figuring out which environmental elements influence vegetation dynamics and how important they are about one another. Based on the differential absorption of the red and near-infrared spectral bands [17,18], the Normalized Difference Vegetation Index (NDVI) measures the greenness of vegetation about its vigor and extent [19,20]. The soil health, geomorphology, CO₂ levels, nitrogen deposition, and climatic restrictions are all closely related to NDVI [21,22]. As a result, there are serious health dangers when metals that have accumulated in soils are concentrated in animal and human organs through the food chain [23]. Furthermore, the PTEs have the potential to contaminate groundwater aquifers by leaking into them [24,25].Therefore, to be able to create an appropriate remediation plan and lessen adverse effects, a precise assessment of soil pollution based on PTEs is essential [16]. The first stages in effectively managing soil pollution are identifying the sources of contamination and comprehending the geographical distribution of heavy metals. Geographic information systems (GISs) thus aid in the spatial distribution mapping of soil attributes [26,27]. An approach that makes it possible to analyze spatial data and subsequently forecast the location of unsampled data is geostatistical analysis. Two widely used geostatistical analysis methods are Kriging and Inverse Distance Weighting (IDW) [28,29]. Numerous techniques, such as the index methodology, quotient method, fuzzy comprehensive assessment, geoaccumulation index, potential ecological risk index, and pollution load index, are used to evaluate soil ecological risk [30,31]. Assessing soil contamination requires an understanding of the relationships between heavy metals, particularly when there are a lot of data. Certain elements frequently have a greater influence on pollution than others [32]. To model the regressions between the various variables, principal component analysis (PCA) can help reduce a large number of variables to a small number of comprehensive principle components [33,34]. Hierarchical cluster analysis (HCA) looks at the separations across samples, grouping the most comparable points into a single cluster. HCA is a method for unsupervised classification that iteratively combines the closest two clusters. The results of PCA may aid in the comprehension of soil variable trends and patterns, find appropriate ways to clean up contaminated soil and create bold plans to attain sustainable development and raise the standard of food production [35]. This study’s primary goal is to answer the following research questions: (1) Where do the majority of the heavy metals (As, Co, Cr, Cu, Fe, Ni, Pb, V, and Zn) in some areas of Jazan, Southwestern Saudi Arabia, come from? (2) Does the soil exhibit any gradient in the levels of heavy metal pollution? (3) What causes the metals to be concentrated in the soil? (4) Is there a suitable way to lessen or overcome the intensity of pollution? In order to improve future management, the spatial distribution of these metals in the research region can help identify their origins and evaluate the mobility of the contaminants. The information gathered from this study may be useful in giving environmental organizations crucial information to update their programs and policies for guaranteed food safety and quality, as well as raising community awareness of the dangers of environmental contamination.

2. Materials and Methods

2.1. Description of the Study Area

The study area is situated between latitudes 16°40′0′′ and 16°50′30′′ N and longitudes 42°50′0′′ and 43°1′30′′ E in the subtropical zone with an area of 174.85 km² (Figure 1) and has an average monthly temperature ranging between 25.8 °C in January and 33.4 °C in July. The average relative humidity ranges from 55% and 72.5% [36]. The yearly rainfall average is roughly 206 mm. The three main periods that comprise the geologic setting of Jazan are the late Proterozoic rocks of the Arabian Shield, the igneous and sedimentary rocks that were deposited in the Red Sea basin between the middle Tertiary and the present, and the remains of sedimentary rocks that were deposited on the Shield between the Cambrian and early Tertiary periods [37].The late Proterozoic schists that were rich in quartz and highly foliated to gneissic granite make up the rocks of the Arabian Shield. The Cambrian to lower Tertiary rock deposition was regulated by changes in the locations of the nearby marine basins’ shorelines [38]. The groundwater level is less than two meters deep, and the Jazan ground surface is composed of deposits of silty clay and soft sand. Additionally, Jazan has been affected by earthquakes throughout history and is situated in an earthquake-prone region of the southern Red Sea earthquake-active zone [39].

The area produces a wide range of fruits, vegetables, and food grains, making it one of the nation’s top agricultural regions. The city also boasts a large power plant, a seaport, an oil port, a desalination plant, and enormous car factories in addition to the copper, aluminum, and steel industries [40,41]. The Digital Elevation Model (DEM) emerged from the SRTM (the NASA Shuttle Radar Topographic Mission) and has a 30 m spatial resolution (https://earthexplorer.usgs.gov, accessed on 1 March 2024).

2.2. Calculation of the Normalized Difference Vegetation Index (NDVI)

Sentinel-2 imagery was used in Google Earth Engine (GEE) to calculate the NDVI, which provided a high-spatial-resolution evaluation of vegetation density and health. Sentinel-2 is a satellite sensor that was ideal for NDVI analysis because it can record multispectral data [42]. Equation (1) presents this dimensionless index, which shows the density and presence of vegetation and may distinguish between several agricultural traits within a field, including plant height and productivity [18,43]:

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

(1)

where NIR is the near-infrared reflectance and R is the red reflectance.

The NDVI revealed a significant correlation between soil characteristics and vegetation [44]. The median of the NDVI was calculated from January 2024 to March 2024.

2.3. Identification of Land Use and Land Cover (LU/LC) Within the Study Area

The ability to identify and recognize different forms of land cover and land use is a critical function of remote sensing. In the past few years, machine learning classification techniques for crop classification have been evolving. Google Earth Engine is a cloud-based platform which offers users the capacity to analyze different satellite datasets using a range of advanced classification techniques [45,46,47]. Using Sentinel-2 (March 2024) (10 m) high-resolution optical data and the Random Forest (RF) approach, crops were classified [48,49]. The model’s accuracy was evaluated using 5000 ground truth data points, of which 70% were reserved for model calibration, allowing for the fine-tuning of features and parameters. By ensuring that the model can consistently generalize its findings beyond the calibration dataset, the calibration−validation split enhances the model’s credibility. The remaining 30% was reserved for the validation of the model. For agricultural planning, management, and observation, this technology aids in the more accurate and dependable identification of crop trends [50].

2.4. Field Work and Laboratory Analysis

To assess the level of PTEs contamination in the soil, 35 surface soil samples were randomly collected from the study area in 2024 (Figure 1). A firm plastic hand trowel was used to gather the samples at a depth of less than 10 cm. A representative sample was created at each station by combining four subsamples into a composite sample, which was then sealed in plastic bags and kept in an ice box. The samples were then sieved through a 2 mm mesh, and large rocks and organic material were removed from the samples in the lab. They were then allowed to air-dry before being manually broken up with an agate mortar [12]. The ALS Geochemistry Laboratory in Jeddah, Saudi Arabia, employed inductively coupled plasma atomic emission spectrometry (ICP-AES) to analyze As, Co, Cr, Cu, Fe, Ni, Pb, V, and Zn. In a graphite heating block, 0.50 g of each sample was incubated with HNO₃-HCl aqua regia for 45 min. The resulting solution was cooled and then diluted to 12.5 mL using deionized water, mixed and analyzed. The linearity, limits of quantification (LOQs), and limits of detection (LODs) of the ICP-AES method were evaluated (0.99; Table S1) [4]. To ensure the accuracy of the analysis, three samples were examined twice. The ALS Geochemistry Laboratory employed a typical analytical batch that included two certified reference materials from Western Australia (CRM11:EMOG-17_24112236 and CRM2:GBM321-8_24112236) to ensure data precision before release and a reagent blank to measure the background. The recovery percentages ranged from 85% to 100.51% (Tables S2–S4).

2.5. Evaluation of Heavy Metal Contamination

Examples of single indices are EF, CF, and Igeo [51]. They can separate pollution from lithologic sources and evaluate the ecological and environmental consequences of contamination [52]. In an effort to restore this resource, geochemical evaluation of soils using suitable indicators and pollution indices has drawn a lot of interest recently [53].

2.5.1. The Geoaccumulation Index (Igeo)

The geoaccumulation index is a widely used metric to assess the level of heavy metal contamination in sediments. Muller [54] introduced the geoaccumulation index, which may be computed using the following Formula (2) and compares current concentrations with background values to assess metal pollution in sediments:

$$I_{geo}={\log }_2\left[{\displaystyle \frac{Cn}{1.5Bn}}\right]$$

(2)

where Bn represents the trace element’s geochemical background concentration (middle crust) and Cn represents the trace element concentration as measured in the soil [55]. To reduce the impact of possible variations in background values that can be caused by geological variation in sediments, the constant 1.5 was added to Equation (2). Since soil is a layer of the Earth’s crust and its chemical composition is linked to that of the crust, the focus here is on the relationship between the concentration of elements in the crust and the concentration obtained from current analysis [56]. In Table S5, the Igeo categorization is displayed. The average Earth’s upper crust value, according to Wedepohl [57], serves as the background value.

2.5.2. Contamination Factor (CF)

In the current study, the contamination status of the soil is ascertained by the contamination factor and degree of contamination. The factor of contamination was determined using the following Equation (3):

$$CF=Cm/Cn$$

(3)

where $Cm$ is the measured concentration of a metal in soil samples and $Cn$ is the value of the metal in background values. In Table S6, the $CF$ categorization is displayed [58].

2.5.3. Enrichment Factor (EF)

According to Franco-Uria [59], EF is regarded as an important tool for assessing the extent of environmental pollutants. Geochemical normalization of the heavy metals data to a conservative element, like Fe, was used to detect aberrant metal concentrations [60,61]. Calculation of the EF was conducted according to the following Equation (4):

$$EF=\left(M÷Fe\right)sample/\left(M÷FE\right)background$$

(4)

The ratio of metal and iron concentrations in a sample is called the (M/Fe) sample, while the ratio of metal to iron concentrations in a backdrop is called the (M/Fe) crust, where M is the metal concentration. While EF values of >2 imply that the sources are more likely to be anthropogenic, EF values of <2 show that the metal comes exclusively from crustal materials or natural processes [62,63,64]. Iron has also been included in this work as a cautious tracer to distinguish between anthropogenic and natural components. The average value of the Earth’s upper crust, as reported by Wedepohl [57], serves as the background value. In Table S7, the EF categorization is displayed according to Sutherland [65].

2.5.4. The Modified Degree of Contamination (mC_d)

The modified degree of contamination (mC_d) denotes the integrated metal pollution status for each location, taking into account the contamination factors for all nine elements examined. The modified degree of contamination (mC_d), which was developed by Abraham and Parker and Cheng et al. [66,67], was used to evaluate the total contamination of several heavy metals in each soil sample. To calculate mC_d, the following Equation (5) was utilized:

$${mC}_d=\sum {\displaystyle \frac{CF}{n}}$$

(5)

where $CF$ is the contamination factor for each element and $n$ is the number of elements under study.

The terms used to categorize the mC_d are shown in Table S8.

2.5.5. Statistical Analysis

The statistical analysis of the investigated heavy metals and soil parameters was performed using SPSS version 25, and Z-scores were used to normalize the variables [68]. Principal component analysis (PCA) was used to separate the dataset into PC variables to remove multicollinearity between the original variables. Linear relationships between the analyzed variables were displayed using the Pearson correlation coefficient. The Kaiser–Meyer–Olkin (KMO) method was used to evaluate the samples’ appropriateness for PCA. If the KMO values were higher than 0.5, the data were suitable for PCA [69] (Table S5). The Bartlett’s test was used to further validate the data’s fitness for PCA, and the findings showed that p < 0.05 [70].

2.6. Spatial Data Analysis of Studied Heavy Metals

To determine the different patterns of specific soil heavy metals (As, Co, Cr, Cu, Fe, Ni, Pb, V, and Zn), the spatial distribution maps of soil heavy metals were created using the intersection Kriging method. ArcGIS Spatial Analyst Pro 3.4 incorporates the Kriging interpolation algorithm into its geostatistical analyses [71]. To find the best model matched for the chosen heavy metals, the current study employed semi-variogram models (Gaussian, Exponential, Stable, and Spherical). According to Johnston et al. [72], the mean standardized error (MSE) and root mean square standardized error (RMSSE) have been used to assess the correctness of the various models. Higher model accuracy is indicated by MSE values that are closer to zero and RMSSE values that are closest to one, and vice versa [73].

3. Results and Discussion

3.1. Vegetation Status (NDVI) and Land Use/Land Cover (LULC) for the Study Area

The research area has relatively high NDVI values, ranging from −0.08 to 0.78 (Figure 2a). According to Linage [74], bare soil is defined as having an NDVI value between 0 and 0.2, while cultivated land is defined as having a value between 0.2 and 0.7. Water has negative values associated with it. The southeast portion of the study area usually has high NDVI values (Figure 2a). The land use within the research region is depicted in Figure 2b. Trees, field crops, grasses, shrubs, urban areas, and bare soils were the six categories noted. Field crops make up about 71.30% of the total study area, while bare soils make up about 3%. Our findings are consistent with those of Najmi et al. [75], who noted that the Jazan region has a wealth of vegetation and that growing a range of fruits, vegetables, and food grains is a popular activity.

3.2. Land Surface Parameters and Heavy Metal Concentrations in the Study Area

The analysis of the Digital Elevation Model (DEM) revealed that the ground surface elevation varies from approximately 59 m above sea level (asl) to 103 m (asl) (Figure 3a). The study research slope ranged from flat (0–0.2%) to steep (30–60%) [76,77] (Figure 3b). The concentrations of nine PTEs in 35 surface soil samples are summarized in

Table 1. Quantitative analysis of selected PTEs within the study area.

	As	Co	Cr	Cu	Fe	Ni	Pb	V	Zn
Measuring unit	mg kg−1
N	35
Minimum	2	5	14	8	13,300	11	2	31	20
Maximum	6	22	67	51	48,600	54	11	109	106
Mean	2.94	12.31	41.17	24.11	32,508	30.17	4.97	73.00	50.40
STD.	0.93	4.51	11.39	10.84	8241	10.37	1.85	18.97	20.03

Note: STD.—standard deviation.

. Fe is the most prevalent PTEs (average of 32508 mg kg⁻¹), followed by V (73 mg kg⁻¹), Zn (50.40 mg kg⁻¹), Cr (41.17 mg kg⁻¹), Ni (30.17 mg kg⁻¹), Cu (24.11 mg kg⁻¹), Co (12.31 mg kg⁻¹), Pb (4.97 mg kg⁻¹), and As (2.94 mg kg⁻¹) (Table 1).

All concentrations of PTEs increased from west to east and south of the study area (Figure 4a–i). According to the findings, the spherical model fitted Co, Cr, Cu, Ni, and Zn, whereas the Gaussian model fitted as well. For Fe and V, an exponential model was fitted. Furthermore, Pb fitted the Stable model. The MSE values were nearly zero in every metal that was chosen, and the RMSSE values were nearly one (

Table 2. Semi-variogram parameters for modeling spatial data.

Variables	Measuring Unit	Models	ME	MSE	RMSEE
As	mg Kg−1	Gaussian	0.0002	−0.004	1.07
Co		Spherical	−0.033	−0.002	0.98
Cr		Spherical	−0.004	0.000	0.99
Cu		Spherical	0.10	0.02	0.97
Fe		Exponential	−12.60	0.000	0.97
Ni		Spherical	−0.04	−0.003	1.00
Pb		Stable	0.03	0.009	1.01
V		Exponential	−0.35	−0.017	0.94
Zn		Spherical	0.29	0.019	1.00

). All PTEs content varied greatly (STD. > 1), according to the data except As, which presented a high degree of uniformity in the research area (STD. < 1.00) (Table 1) [78].

Table 3. Comparison of the PTEs concentrations in the current study with background values, official quality thresholds, and PTEs from other areas.

	As	Co	Cr	Cu	Fe	Ni	Pb	V	Zn
Present study	2.94	12.31	41.17	24.11	32,508	30.17	4.97	73.00	50.40
Jazan coastal area, Red Sea, Saudi Arabia [83]	-	4.10	33	31.60	-	20	2.30	-	28.50
Al-Khobar, Arabian Gulf, Saudi Arabia [82]	1.60	4.80	51	183	-	75	5.40	-	52.70
Arabian Gulf, Saudi Arabia [81]	31	4.00	64	297	-	77		-	48.30
Background values [57]	2.00	11.60	35	14.30	30,890	18.60	17	53	52
Average natural concentration of heavy metals in rocks [80]	5.5–12	1.3–10	15–70	9.9–39	140,00–28,000	1.8–18	2.6–27	20–93	37–68
Recommended concentration based on the Department of Environmental Affairs (DEA) [79]	5.8	300	6.5	16	-	91	20	150	240

compares the average levels of current metal concentrations in the study area with those of other coastal regions in Saudi Arabia and Bradi’s, Wedepohl’s, and the Department of Environmental Affairs (DEA) recommended concentrations [57,79,80,81,82]. The average As concentration was lower than other values shown in Table 3 [79,80,81] and greater than those from the Arabian Gulf, Al-Khobar [82], and Wedepohl’s [57] recommended concentrations. Co concentration was higher than all other values included in Table 3, except the concentrations suggested by the DEA [79]. However, aside from Wedepohl’s [57] suggested concentrations, the average of Cu was lower than the other values listed in Table 3. The average level of Fe exceeded the other levels listed in the Table 3. The Ni value was lower than the values obtained from Alharbi and El-Sorogy [82] and El-Sorogy et al. [81] and higher than the values obtained from Kahal et al. [83] and Wedepohl’s [57] and Bradi’s [80] recommended concentrations. The V concentration in the current study was higher than recommended values [57] and lower than the other values listed in Table 3. The average Zn concentration was lower than the values listed in Table 3, except the values found by El-Sorogy et al. [81] in the Arabian Gulf, Saudi Arabia, [83] in the coastal region of Jazan, Red Sea, Saudi Arabia. Some PTEs, including Fe, Ni, and Zn, are necessary for vital nutritional processes but only require small amounts. However, an excessive amount of these PTEs can lead to serious health issues like diabetes, neurological and renal conditions, and cardiovascular diseases [84].

3.3. Single Contamination Indices

The degree of contamination can be measured using the contamination factor (CF) [58]. Except Pb and Zn, which had low contamination (less than 1), the results of the CF study showed that all of the PTEs present in the soil under investigation had moderate CF (1.05–1.69;

Table 4. Statistical analysis of individual PTE contamination indices in the study area.

Contamination Indices					CF
PTEs	As	Co	Cr	Cu	Fe	Ni	Pb	V	Zn
Minimum	1.00	0.43	0.40	0.56	0.43	0.59	0.12	0.58	0.38
Maximum	3.00	1.90	1.91	3.57	1.57	2.90	0.65	2.06	2.04
Mean	1.47	1.06	1.18	1.69	1.05	1.62	0.29	1.38	0.97
STD.	0.47	0.39	0.33	0.76	0.27	0.56	0.11	0.36	0.39
					EF
Minimum	0.82	0.73	0.91	1.04	1.00	1.14	0.14	1.13	0.59
Maximum	2.32	1.21	1.39	2.27	1.00	2.00	0.52	1.77	1.30
Mean	1.41	0.98	1.11	1.54	1.00	1.51	0.28	1.31	0.90
STD.	0.30	0.13	0.11	0.34	0.00	0.21	0.09	0.14	0.17
					Igeo
Minimum	−0.58	−1.80	−1.91	−1.42	−1.80	−1.34	−3.67	−1.36	−1.96
Maximum	1.00	0.34	0.35	1.25	0.07	0.95	−1.21	0.46	0.44
Mean	−0.09	−0.60	−0.41	0.03	−0.56	0.02	−2.45	−0.17	−0.74
STD.	0.44	0.54	0.44	0.66	0.39	0.52	0.53	0.40	0.57

Note: STD.—standard deviation.

). The origin of heavy metals can be ascertained with the help of the enrichment factor (EF) [85]. The PTEs’ average EF values, displayed in descending order, are as follows: Cu (1.54) > Ni (1.51) > As (1.41) > V (1.31) > Cr (1.11) > Fe (1) > Co (0.98) > Zn (0.9) > Pb (0.28) (Table 4). This implies that the research area had depletion to minimum enrichment, while background concentrations of Co, Pb, and Zn were present [86]. The study area soil samples appeared to be uncontaminated based on the average Igeo for the examined PTEs (Table 4).

3.4. Multivariate Analysis of PTEs in the Investigated Area

All heavy metals presented a strong positive association, suggesting that these elements share a common source [87]. Strong positive relationships between iron and other elements suggested that it originated naturally, specifically from the chemical weathering of the basement rocks (made up of igneous and metamorphic rocks and found beneath sedimentary layers) (

Table 5. Correlations relations between PTE_S in the study area.

Correlations
	As	Co	Cr	Cu	Fe	Ni	V	Pb	Zn
As	1
Co	0.831 **	1
Cr	0.810 **	0.934 **	1
Cu	0.828 **	0.977 **	0.917 **	1
Fe	0.780 **	0.971 **	0.930 **	0.951 **	1
Ni	0.824 **	0.974 **	0.972 **	0.956 **	0.943 **	1
V	0.656 **	0.897 **	0.836 **	0.846 **	0.922 **	0.852 **	1
Pb	0.574 **	0.591 **	0.621 **	0.681 **	0.547 **	0.625 **	0.343 *	1
Zn	0.826 **	0.928 **	0.915 **	0.960 **	0.893 **	0.939 **	0.731 **	0.773 **	1

** Correlation is significant at the 0.01 level; * correlation is significant at the 0.05 level.

), which located in the east of the research region, in the adjacent Arabian Shield mountains [88,89].The significance level is less than 0.001, and the results of the Kaiser−Meyer−Olkin (KMO) and Bartlett’s tests of sphericity were 0.88 and 613.96, respectively (Table S9). These numbers show that PCA is appropriate for this investigation [90]. Only one of the components with a cumulative variability of 82.72% had eigenvalues larger than 1, as shown in the data in

Table 6. PCA of PTEs in the study area.

PTEs	PC1
As	0.862
Co	0.985
Cr	0.963
Cu	0.984
Fe	0.966
Ni	0.982
Pb	0.682
V	0.865
Zn	0.963
Eigenvalue	7.44
Variability (%)	82.72
Cumulative (%)	82.72

, and the other components were eliminated based on a Kaiser test [91]. The principal component analysis yielded one principal component (Table 6), which showed loadings for each PTE and thus validated the correlation matrix of the geogenic source for the PTEs under investigation [2]. The data in this study were separated into two primary groups (clusters) by the HCA. The dendrogram in Figure 5 illustrates the differences between the two clusters, which include all of the sites that were further subdivided into other clusters. There were 18 observations in the first cluster and 17 in the second cluster, with varying ranges, averages, and STDs for every variable. Significant concentrations in Co, Cr, Cu, Fe, Ni, V, and Zn were found in the collected analysis. In contrast, there was no discernible difference in the total As and Pb concentrations between the two clusters (

Table 7. Mean concentrations of PTEs in two clusters (C1 andC2).

Mean of PTEs	N	As	Co	Cr	Cu	Fe	Ni	Pb	V	Zn
C1	18	3.50 a	15.94 a	49.61 a	32.78 a	39,150 a	38.22 a	5.67 a	87.28 a	65.67 a
C2	17	2.35 a	8.47 b	32.24 b	14.94 b	25,476 b	21.65 b	4.24 a	57.88 b	34.24 b

Note: significant differences are indicated by the means of variables with different letters.

). Cluster 2 can be explained by the fact that PTEs vary less in these locations. However, cluster 1 can be ascribed to the significant PTE variances in the east of the study area. This is because PTEs and other metallic minerals are abundant in the Arabian Shield’s basement rocks in Saudi Arabia. These include chromite (FeCr₂O₄) in peridotite units, vanadiferous magnetite (Fe,V)₃O₄ connected to mafic to ultramafic rocks, sphalerite (ZnS) in hydrothermal vein systems, and chalcopyrite (CuFeS₂) connected to intrusive bodies and volcanic rocks. Pentlandite ((Fe,Ni)₉S₈) is also found in mafic–ultramafic complexes, galena (PbS) is linked to hydrothermal veins, arsenopyrite (FeAsS) is found with sulphide mineralization in hydrothermal systems, and both magnetite (Fe₃O₄) and haematite (Fe₂O₃) are found in sedimentary and volcanic rocks, which were eroded after being weathered by rainfall [92,93].

3.5. Percentages of Single Contamination Indices Across the Clusters

The Igeo results indicated that all samples were uncontaminated. Combinations of cluster and CF analysis can be used to assess the degree of contamination by heavy metals in soil [28]. Both clusters (C1 and C2) presented that the CF fell within the range of 1 ≤ CF < 3, meaning that the samples under study exhibited a moderate degree of As contamination. (Figure 6a,b). The CF of Co indicated that 5.55% of samples were low-contaminated and 94.44% were moderately contaminated in cluster 1 (C1), while in cluster2 (C2) they were 100% moderately contaminated (Figure 6a,b). The results showed that C2 had two different degrees of Cr contamination, with about 41% for moderate contamination and 59% for low contamination. In contrast, the CF of Cr in cluster 1 was 100% moderate contamination. Figure 6a,b show that C1 recorded moderate Cu (CF) levels in 88.88% of soil samples and considerable Cu (CF) values of 11.11%. In contrast, 47.06% of soil samples in C2 showed moderate contamination, while 52.94% had low contamination. All of the soil samples in C1 had moderate levels of Fe contamination, according to the data, whereas C2 had low levels in 94.11% and moderate levels in 5.89% of the samples. According to Hökanson [58], for Ni and V, C1 demonstrated that CF belongs to the moderate class. Furthermore, 23.52% of samples had low contamination levels and 76.48% of samples had moderate contamination levels for V in C2, while 29.41% and 70.59% of soil samples had low and moderate contamination levels of Ni in C2, respectively. Every soil sample in clusters 1 and 2 had low levels of Pb contamination (Figure 6a,b). The results demonstrated that the CF of Zn in both clusters displayed two distinct levels of Zn contamination. About 5.55% of the samples had low contamination, and 94.96% of the samples had moderate contamination in C1. In C2, 94.11% of the samples had low contamination, and 5.89% of the samples had moderate contamination.

The EF of As indicated that C1 had low enrichment and moderate enrichment, with percentages of 94.11 and 5.89%, respectively. Furthermore, C2 showed a relative proportion of 88.89% for low enrichment and 11.11% for natural origin. In C1, the EF of Co were 38.88% natural origin and 61.12% low enrichment, whereas in C2 they were 88.23 and 11.77%, respectively. While the EF values of Cr in C2 were 11.76% natural origin and 88.24% low enrichment, the values in C1 were 22.22% natural origin and 77.78% low enrichment (Figure 6c,d). The findings showed that all samples in C2 were classified as low enrichment, while 88.88% of the cu soil-concentration samples in C1 had low enrichment and 11.12% had moderate enrichment. All of the C1 samples belonged to the low enrichment class for Ni, whereas 94.11% of the C2 samples belonged to this class with the remaining samples showing moderate enrichment. The findings showed that while the EFs of V in clusters 1 and 2 were 100% low enrichment, the EFs of pb in clusters 1 and 2 were 100% natural origin (Figure 6c,d). The findings demonstrated that the EF of Zn in C1 showed that 44.46% of soil samples had low enrichment and 55.55 percent had natural origin. On the other hand, C2 showed comparatively low contamination; 94.11% of its soil samples were of natural origin, while 5.89% belonged to the low enrichment class. These results showed that categorizing the study region into two zones—each with a different concentration and pattern of heavy metals—was made possible by a unique outcomes of integrating PCA and HCA [35].

Modified Degree of Contamination (mC_d)

According to the mC_d results calculated in each location, the overall study area can be divided in three sub-regions. The most extensive one, spanning 148.20 km², is characterized by nil to a very low degree of contamination and there are two smaller areas, 26.16 and 0.495 km², respectively, with low and moderate degrees of contamination. According to the NDVI and LULC data, the research area is highly cultivated, so it is crucial to monitor PTEs contamination there since the natural weathering process may be the source of PTEs. Future contamination may rise as a result of human activity, particularly in the 0.5 km² of the study region that has moderate contamination. However, the geogenic source for the PTEs was under investigation (Figure 7).

4. Conclusions

One of the most significant barriers to sustainable development and food security is the assessment of soil PTEs, which is the focus of the current study. This investigation showed that the spatial distribution maps of heavy metals in the study area were successfully predicted by this study using semi-variogram models. Moreover, the combination of PCA and HCA produced unique outcomes by dividing the study region into two zones, each with a different PTE content and pattern. The mC_d values in the research region often fall into three classes, with areas of 148.20, 26.16, and 0.495 km², respectively: nil to very low degree of contamination, low degree of contamination, and moderate degree of contamination. The findings revealed notable PTE variations in the research area’s eastern region. This is due to the abundance of PTEs and other metallic minerals in the basement rocks of Saudi Arabia’s Arabian Shield. To sum up, it is critical to regularly evaluate the condition of the soil in order to stop the study area from becoming more contaminated. Human activity may cause future contamination to increase, especially in the 0.5 km² of the research zone with moderate contamination and the extensively cultivated study area. It is recommended that the sampling size be increased to enhance the geographic interpretation of soil contamination.