Characteristics, Source and Risk Assessment of Soil Polycyclic Aromatic Hydrocarbons around Oil Wells in the Yellow River Delta, China

Abstract

The Yellow River Delta (YRD) is the most complete wetland ecosystem in the warm temperate zone of China and is rich in oil resources. However, with petroleum extraction and the development of the economy, pollution of the YRD has been paid increasing attention, in particular, pollution via polycyclic aromatic hydrocarbons (PAHs), as they have caused great harm to human health and the ecosystem balance. Based on the investigations of a research group in 2009, this study re-collected samples according to the same sampling points and analyzed the concentration, composition, source, ecological risk and health risk of PAHs in 2021. The concentration of ΣPAH₁₆ in the surface soil of YRD in 2009 ranged from 2.6 to 8275.46 ng/g, with an average of 1744.41 ng/g. The concentration of ΣPAH₁₆ in 2021 ranged from 56.25 to 582.56 ng/g, with an average of 149.63 ng/g. Therefore, the pollution situation in the YRD in 2021 was significantly improved compared with 2009. The composition of PAHs in soil in 2009 and 2021 was similar, which was dominated by low-ring PAHs. The evaluation results of the toxicity equivalent factor method showed that there was no potential ecological risk in the soil in 2009 and 2021. The evaluation results of the lifetime cancer risk increment model showed that the incremental lifetime cancer risk models (ILCRs) of soil PAHs in 2009 and 2021 were lower than the safety threshold of 10⁻⁶; therefore, there was no carcinogenic risk. The existing management measures for oil wells need to be further promoted to protect the regional ecological environment in the YRD.

1. Introduction

The Yellow River Delta (YRD) is the most primitive, youngest and most complete wetland ecosystem in the warm temperate zone of China and even the world [1]. There are abundant mineral resources dominated by oil in the Yellow River Delta. With the in-depth development of oil resources and the production and construction of other industries, the ecological environment of the Yellow River Delta has also undergone great changes [2]. The discharge of a large number of pollutants seriously threatens the safety of the ecosystem in the region [3,4], among which the problem of oil pollution is particularly prominent. Previous studies have shown that PAH pollution in surface water and groundwater was low, and the soil and sediment pollution was moderate. The former pollution was mainly from oil input and oil combustion, and the latter was mainly from wood, coal combustion and oil input [5].

The pollution of polycyclic aromatic hydrocarbons (PAHs) in the YRD has attracted the attention of scholars. PAHs are a class of persistent organic pollutants that are widely present in various environments [6]. The sources of PAHs in the environment are mainly divided into natural sources and anthropogenic sources [7]. The main cause of PAH pollution at present is anthropogenic sources, including petroleum and pyrolysis sources. The large quantity of PAHs produced by human activities usually enter the atmospheric environment first and then eventually converge into the soil, sediment and water environment through surface runoff, dry and wet deposition, etc. [8]. PAHs can also be bioenriched, accumulating and transferring into the food chain, thus increasing the risk level [9]. PAHs are not only carcinogenic, teratogenic and mutagenic but also phototoxic [10,11], causing great harm to human health.

Previous pollution investigations in soil and sediment have been performed in the YRD [10,12,13,14]. For example, the concentration of PAHs in 61 surface soil samples from the YRD has been studied. Results have shown that the total concentration of PAHs ranged from 27 to 735 ng/g, with the highest concentration occurring in the central and southern YRD, which was related to oil exploration. The PAHs in the surface soil of the YRD were mainly from coal and biomass burning, oil leakage and/or vehicle emissions (Yuan et al. 2014). In addition, it was found that the concentration of PAHs in soil around oil wells at different ages increased with the increase in mining time. From the source point of view, the soil near the new oil wells was mainly polluted by oil sources, while the old oil wells were mostly related to combustion sources [15]. The distribution of coastal sediments in the YRD is different. The concentration of PAHs in surface sediments near the Yellow River estuary is lower than that in the northern and southern coastal wetlands of the YRD [16]. Other studies have found that the PAH pollution level in the Yellow River Delta Nature Reserve is relatively low, and the potential ecological risk level is extremely low, but there is a potential carcinogenic risk. In addition, the spatial distribution characteristics of PAH concentration and PAH health risk value are similar, which indicates that the carcinogenic risk of PAHs is proportional to its content [17]. These research results have important theoretical significance for the protection of the ecological environment and the management of the Yellow River Delta.

At present, the study of oil pollution in the YRD is mainly the result of single sampling observations, which renders it difficult to reflect the long-term change characteristics of pollution in the YRD. There are few studies on the distribution of PAHs and their changes in different years in this area. Based on the research group’s survey in the Yellow River Delta in 2009, this paper collected samples from the same sampling points to analyze the concentration, composition, source, ecological risk and health risk of PAHs and to understand the change in PAH pollution in the region in the recent 10 years.

2. Materials and Methods

2.1. Study Area

The Yellow River Delta is located in East China (118°07′–119°18′ E, 36°55′–38°12′ N) (Figure 1). About 96% of its area is distributed in Dongying City, Shandong Province. It is one of the most important wetland protection sites in the world. It is the breeding habitat and migration transit site of rare birds in Northeast Asia and the Western Pacific. The research area of this paper is mainly located in Xianhe Town, Gudao Town and Huanghekou Town of Dongying City, with a total area of 2152.7 km². Affected by the Eurasian continent and the Western Pacific, the study area belongs to the warm temperate semi-humid continental monsoon climate, with four distinct seasons and a hot summer and cold winter. The annual average temperature is between 11.7 and 12.8 °C. The average annual precipitation is about 530~630 mm, of which the summer rainfall accounts for more than half of the year. The rain and heat are in the same period, which is conducive to the growth of plants and crops and can easily cause flood disasters [18,19]. According to the statistics of soil census data, the study area can be divided into five soil types: cinnamon soil, saline soil, fluvo-aquic soil, mortar black soil and paddy soil, among which fluvo-aquic soil and saline soil account for 95% of the total soil area. The YRD is rich in natural resources, with a large number of brine, oil and gas, geothermal, mineral and other resources. It is an important energy base in China [20].

2.2. Sample Collection

The soil samples in this study were collected in July 2021, and the collection areas were mainly concentrated in Xianhe Town, Gudao Town and Huanghekou Town of Dongying City. The sampling points were collected in accordance with the sampling points of the field survey conducted in 2009 around the wells. At each sample point, soil samples of 0–20 cm of surface layer were taken according to the plum sampling method. The basic physical and chemical properties such as pH and temperature were measured. The samples were mixed evenly; then, about 1kg of the evenly mixed samples were taken by the quartering method and put into a Ziplock bag for marking. A total of 15 soil samples were collected and transported back to the laboratory for freezing (−20 °C) preservation at low temperature and in darkness. The notations O1–O15 are soil points in 2009 and 2021 (Figure 1).

2.3. Laboratory Analysis

2.3.1. Determination of PAHs

According to the HJ805-2016 standard, polycyclic aromatic hydrocarbons in soil were determined by gas chromatography–mass spectrometry. The determination process included sample preparation, extraction, concentration, purification, elution and constant volume determination as follows. After removing impurities such as gravel and roots in the sample, 20 g was weighed and an appropriate amount of anhydrous sodium sulfate was added to dehydrate. The sample was ground into fine particles. About 5 g of samples was weighed and extracted by the Soxhlet extraction method. An amount of 100 mL acetone-n-hexane mixed solvent was added to the sample and extracted for 16~18 h, refluxing 4~6 times per hour. The extract was concentrated to about 2 mL using a parallel concentrator. Then, the samples were purified by a silica gel chromatography column. An amount of 25 mL dichloromethane-pentane mixed solvent was added for elution and the eluent was collected. After the eluent was concentrated again, an appropriate amount of acetone-n-hexane mixed solvent was added, and the volume was fixed to 1 mL for determination. Gas chromatography–mass spectrometry with a quartz column (DB-5MS, 30 m × 0.25 mm × 0.25 μm) was used to quantify the PAH concentrations. The carrier gas type is high-purity helium. The sample was injected without a shunt. The injection volume was 1 μL and the inlet temperature was 280 °C. The oven temperature program was as follows: the initial temperature was 80 °C, held for 2 min, then heated at 20 °C/min to 180 °C, then heated at 10 °C/min to 290 °C, and held for 5 min.

2.3.2. Quality Control

In order to ensure the accuracy and reliability of the data, strict quality control was carried out throughout the experiment. In the determination of PAHs in soil and sediment, not only was o-fluorobiphenyl added to each sample for recovery detection but blank experiments and parallel sample tests were also carried out. The spiked recoveries of 16 PAHs were between 40% and 150%. The blank experiment was performed once every 20 samples, and the relative standard deviations of parallel samples were all within 20%. The external standard method was used to quantitatively analyze PAHs, and the correlation coefficients of the standard curves were all above 0.995.

2.4. PAH Source Apportionment Method

2.4.1. Positive Definite Matrix Factorization Model (PMF)

The most critical step in the source apportionment of the PMF model is the determination of the number of factors. In the robust model, the number of factors is set to 3~6, and 20 iterations are carried out, respectively. Considering the factors such as the residual and fitting effect, the 4-factor model was used to analyze the source of PAHs across the 21 years, and the 6-factor model was used to analyze the source of PAHs in the soil in 2009.

2.4.2. Diagnostic Ratios

Diagnostic ratios are commonly used in the source apportionment of PAHs [21]. Based on the characteristics of similar kinetic mass transfer coefficients and thermodynamic distribution coefficients between isomeric PAHs, the heat source and petroleum source can be distinguished by their concentration ratio [22]. To differentiate the PAHs at our sampling sites with respect to petroleum vs. combustion sources, we applied two isomeric ratios.

2.5. Ecological Risk Assessment Method of PAHs

2.5.1. Toxic Equivalent Factor Method (TEQ)

Among the 16 US EPA priority PAHs, Benzo[a]pyrene (BaP) is used as a marker for estimating the carcinogenic potency of PAHs [23]. To evaluate and quantify the individual PAH and ∑PAH_7c carcinogenic toxic potency, the toxicity equivalency factors (TEFs), which define the carcinogenic potency of the individual PAH relative to Benzo[a]pyrene (BaP), were used to estimate the toxic equivalent of PAH monomer relative to BaP (TEQ_BaP) in this study [24]. This evaluation approach pays attention to the estimation of carcinogenic PAHs. The toxic equivalent factor (TEF) of BaP was set to 1, and the remaining PAHs were assigned to the corresponding TEF according to the difference in toxicity (

Table 1. Toxicity equivalent factors of individual PAHs.

Monomer PAHs	Ring Number	Toxic Equivalent Factor	Monomer PAHs	Ring Number	Toxic Equivalent Factor
Naphthalence (Nap)	2	0.001	Pyrene (Pyr)	4	0.001
Anthracence (Ant)	3	0.01	Fluoranthene (Fla)	4	0.001
Phenanthrene (Phe)	3	0.001	Dibenzo[a,h]anthracene (DahA) *	5	1
Fluorene (Flu)	3	0.001	Benzo[a]pyrene (BaP) *	5	1
Acenapthene (Ace)	3	0.001	Benzo[k]fluoranthene (BkF) *	5	0.1
Acenaphthylene (Acpy)	3	0.001	Benzo[b]fluoranthene (BbF) *	5	0.1
Chrysene (Chr) *	4	0.01	Benzo[g,h,i]perylene (BghiP)	6	0.01
Benzo[a]anthracene (BaA) *	4	0.1	Indenol [1,2,3-cd]pyrene (IcdP) *	6	0.1

Note: * means carcinogenic compound.

). The TEFs listed in Table 1 are the most reasonable values for the relative titers of PAHs currently available and using them would significantly reduce the uncertainty in risk assessments involving PAHs [24]. The calculation formula is as follows:

$${\mathrm{T}\mathrm{E}\mathrm{Q}}_{\mathrm{B}\mathrm{a}\mathrm{P}}=\sum {\mathrm{C}}_{\mathrm{i}}\times {\mathrm{T}\mathrm{E}\mathrm{F}}_{\mathrm{i}}$$

(1)

where C_i is the exposure concentration of the ith PAHs; TEF_i is the toxic equivalent factor corresponding to the ith PAHs. The greater the toxic equivalent factor of monomer PAHs, the greater the toxicity.

2.5.2. Effect Interval Low Median Method (ERL/ERM)

The evaluation criteria proposed by Long et al. [25] are an effective tool to identify the potential ecological risks of PAHs [26,27]. The tool uses the effects range low (ERL) and effects range median (ERM) as thresholds; if the PAH concentration is less than ERL, the probability of adverse toxic effects of PAHs on organisms is less than 10%, and ecological risks rarely occur. ERL greater than or equal to ERM indicates that PAHs have more than 50% probability of causing negative effects on the ecological environment, and risks may occur frequently. Between the two, the PAH concentration will occasionally harm the ecological environment; that is, the concentration has potential ecological risks.

In order to quantitatively predict the comprehensive ecological toxicity of various PAHs, the ecological risk can be further evaluated by the mean quotient of ERM (M-ERM-Q). The calculation formula is:

$$\mathrm{M}-\mathrm{E}\mathrm{R}\mathrm{M}-\mathrm{Q}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{C}}_{\mathrm{i}}}{{\mathrm{E}\mathrm{R}\mathrm{M}}_{\mathrm{i}}}/\mathrm{n}$$

(2)

where C_i is the exposure concentration of individual PAHs; ERM_i is the median effect range (ERM) corresponding to individual PAHs; and n is the type of PAHs. When M-ERM-Q < 0.1, the possibility of ecological risk is low; when 0.1 < M-ERM-Q < 0.5, there is a 30% probability of toxicity; when 0.5 < M-ERM-Q < 1.5, the probability of medium and high toxicity is about 50%; when M-ERM-Q > 1.5, the ecological risk is high, and the probability of high toxicity is about 75%.

2.6. Health Risk Assessment Methods of PAHs

Incremental Lifetime Cancer Risk Models (ILCRs)

The incremental lifetime cancer risks (ILCRs) model was used in this paper to quantitatively evaluate the health risk of PAHs in soil, that is, the incidence of cancer caused by a certain dose of carcinogens at a certain time. There are three main ways for humans to be directly exposed to soil PAHs: ingestion (ILCR_ingestion), inhalation (ILCR_inhalation) and dermal contact (ILCR_dermal) [28]. The cancer risk calculation formula for each pathway is as follows:

$$\mathrm{C}\mathrm{S}=\sum ({\mathrm{P}\mathrm{A}\mathrm{H}}_{\mathrm{i}}\times {\mathrm{T}\mathrm{E}\mathrm{F}}_{\mathrm{i}})$$

(3)

$${\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\mathrm{C}\mathrm{S}\times \left({\mathrm{C}\mathrm{S}\mathrm{F}}_{\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}\times ∛\frac{\mathrm{B}\mathrm{W}}{70}\right)\times {\mathrm{I}\mathrm{R}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}\times \mathrm{E}\mathrm{F}\times \mathrm{E}\mathrm{D}/\mathrm{B}\mathrm{W}\times \mathrm{A}\mathrm{T}\times {10}^6$$

(4)

$${\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\mathrm{C}\mathrm{S}\times \left({\mathrm{C}\mathrm{S}\mathrm{F}}_{\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}\times ∛\frac{\mathrm{B}\mathrm{W}}{70}\right){\mathrm{I}\mathrm{R}}_{\mathrm{a}\mathrm{i}\mathrm{r}}\times \mathrm{E}\mathrm{F}\times \mathrm{E}\mathrm{D}/\mathrm{B}\mathrm{W}\times \mathrm{A}\mathrm{T}\times \mathrm{P}\mathrm{E}\mathrm{F}$$

(5)

$${\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}=\mathrm{C}\mathrm{S}\times \left({\mathrm{C}\mathrm{S}\mathrm{F}}_{\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}\times ∛\frac{\mathrm{B}\mathrm{W}}{70}\right)\times \mathrm{S}\mathrm{A}\times \mathrm{A}\mathrm{F}\times \mathrm{A}\mathrm{B}\mathrm{S}\times \mathrm{E}\mathrm{F}\times \mathrm{E}\mathrm{D}/\mathrm{B}\mathrm{W}\times \mathrm{A}\mathrm{T}\times {10}^6$$

(6)

$${\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{s}}={\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}+{\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}+{\mathrm{I}\mathrm{L}\mathrm{C}\mathrm{R}}_{\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}$$

(7)

where CS is the toxic equivalent concentration (mg/kg) of monomer PAHs; PAH_i is the exposure concentration of the i th PAHs; TEF_i is the toxic equivalent factor of monomer PAHs (

Table 2. Parameters of incremental lifetime cancer risks.

Parameter	Unit	Adult	Child
BW	Kg	61.5	15
IRsoil	mg·d−1	100	200
EF	d·a−1	350	350
ED	a	24	6
AT	d	25,550	25,550
IRair	m3·d−1	20	10
PEF	m3·kg−1	1.36 × 109	1.36 × 109
SA	cm2·d−1	5700	2800
AF	mg·cm−2	0.07	0.2
ABS	zero dimension	0.13	0.13

); CSF is the carcinogenic slope factor (mg/(kg·d)), which was determined by the carcinogenic ability of BaP. The values of CSF_ingestion, CSF_inhalation and CSF_dermal were 7.3, 3.85 and 25, respectively. BW is body weight, IR_soil is soil inhalation rate, EF is exposure frequency, ED is exposure time, AT is the average life of carcinogens, IR_air is inhalation rate, PEF is soil particle emission factor, SA is skin surface exposure area, AF is skin adhesion factor, and ABS is skin adsorption parameter. The specific parameter values are shown in Table 2. ILCRs are the sum of the carcinogenic risks of the three exposure pathways: when ILCRs < 10⁻⁶, the carcinogenic risk is negligible; when ILCRs > 10⁻⁴, this indicates a high carcinogenic risk; between the two indicates a potential carcinogenic risk.

3. Results and Discussion

3.1. PAH Concentration Characteristics

In 2009, the concentration of PAHs in the samples varied greatly, ranging from 2.6 to 8275.46 ng/g, with an average of 1744.41 ng/g, according to the analysis of surface soil samples collected by the research group in the YRD (

Table 3. Descriptive statistics of PAHs in soils in 2009 (ng/g).

Monomer PAH	Ring Number	Range	Mean Value	Standard Deviation	Coefficient of Variation	Detection Rate
Nap	2	ND~218.67	86.82	58.33	0.67	93
Ant	3	ND~188.81	25.22	46.65	1.85	73
Phe	3	ND~1002.96	117.00	243.43	2.08	93
Flu	3	2.6~7651.99	1409.20	1794.76	1.27	100
Ace	3	ND	ND	ND	ND	ND
Acpy	3	ND~2.16	0.14	0.54	3.74	7
Chr	4	ND~230.93	29.71	58.16	1.96	67
BaA	4	ND~227.56	19.40	56.38	2.91	47
Pyr	4	ND~137.97	24.59	44.82	1.82	80
Fla	4	ND~52.52	16.20	14.82	0.91	80
DahA	5	ND	ND	ND	ND	ND
BaP	5	ND~22.71	2.51	5.94	2.36	40
BkF	5	ND~14.78	2.03	4.80	2.36	20
BbF	5	ND~19.93	4.44	6.67	1.50	60
BghiP	6	ND~86.77	5.82	21.64	3.72	20
IcdP	6	ND~10.77	1.32	3.37	2.56	13
ΣPAH7		ND~481.35	59.42	118.88	2.00	73
ΣPAH16		2.6~8275.46	1744.41	1941.10	1.11	100

Note: ND means “Not Detected”.

). Compared with ΣPAH₁₆ concentrations in the YRD studied in other years, such as Yuan et al. (27–753 ng/g, with an average of 118 ± 132 ng/g) [29]; Yuan et al. (79.2–311 ng/g, with an average of 119 ng/g) [30]; and Fu et al. (278.7–733.5 ng/g, with an average of 382.5 ± 128.4 ng/g) [15], PAH concentrations in soil samples investigated in 2009 were higher on average. Among soil samples, four samples exceeded the mean value, accounting for 26.67% of total soil samples, and eleven samples were below the mean value, accounting for 73.33% of total soil samples. The concentration of seven carcinogenic PAHs ranged from ND to 481.35 ng/g, with an average of 59.42 ng/g, accounting for 3.41% of ΣPAH₁₆. The BaP concentration ranged from ND to 22.71 ng/g, with an average of 2.51 ng/g, respectively, accounting for 4.23% and 0.14% of ΣPAH₇ and ΣPAH₁₆. The other PAH monomers were detected to varying degrees, except for Ace and DahA, and the detection rate was between 7% and 100%. Flu had the highest detection rate, followed by Nap and Phe, and Acpy had the lowest detection rate (Table 3). According to Maliszewska-Kordybach’s [31] definition of soil PAH pollution, only one sample had a concentration below 200 ng/g, accounting for 6.67% of total soil samples. There were two samples with a concentration range of 200–600 ng/g, two samples with a concentration range of 600–1000 ng/g, and ten samples with a concentration higher than 1000 ng/g, respectively, accounting for 13.33%, 13.33% and 66.67% of total soil samples. The above results show that more than half of the soil samples in the study area reached the level of heavy pollution, which suggests the pollution situation is serious and corresponding measures need to be taken to repair and control the situation.

In 2021, the concentration of PAHs in soil samples ranged from 56.25 to 582.56 ng/g, with an average of 149.63 ng/g, collected from the YRD corresponding to the location information in 2009 (

Table 4. Descriptive statistics of PAHs in soils in 2021 (ng/g).

Monomer PAH	Ring Number	Range	Mean Value	Standard Deviation	Coefficient of Variation	Detection Rate
Nap	2	45~291.28	112.79	59.64	0.53	100
Ant	3	ND~10.2	1.30	3.32	2.55	13
Phe	3	ND~10.2	1.30	3.32	2.55	13
Flu	3	ND~8.48	1.63	3.26	2.00	20
Ace	3	ND	ND	ND	ND	ND
Acpy	3	ND	ND	ND	ND	ND
Chr	4	ND~21.58	1.44	5.38	3.74	7
BaA	4	ND~32.36	2.89	8.34	2.88	13
Pyr	4	ND~11.25	1.37	3.51	2.56	13
Fla	4	ND~12.63	2.94	4.91	1.67	27
DahA	5	ND	ND	ND	ND	ND
BaP	5	ND~107.88	13.19	26.18	1.98	53
BkF	5	ND	ND	ND	ND	ND
BbF	5	ND~53.94	5.04	13.58	2.69	20
BghiP	6	ND	ND	ND	ND	ND
IcdP	6	ND~75.52	5.74	18.83	3.28	13
ΣPAH7		ND~291.28	28.30	71.04	2.51	53
ΣPAH16		56.25~582.56	149.63	122.33	0.82	100

Note: ND means “Not Detected”.

). According to the study of PAH pollution in the soil around oil wells in the Yellow River Delta in other years, the PAH concentration decreased with the increase in years, which indicates that the PAH pollution has been effectively managed [15,32] (

Table 5. Concentrations (ng/g) of PAHs in sediments or soils from other research on the YRD.

Sampling Points	Sampling Year	∑PAHs, ng/g	Concentration Mean, ng/g	Reference
61 Stations within the YRD	2006–2008	27–753	118 ± 132	[29]
The Experimental, Buffer and Core area of the YRDNR	2013	79.2–311	119	[30]
Around Oil Wells with Different Extraction Histories in the YRD	2015	278.7–733.5	382.5 ± 128.4	[15]
The Oil Field Soil and the Estuary Soil of the YRD	2018	157.8–481.7	274	[32]
Around the wells in Xianhe Town, Gudao Town and Huanghekou Town of Dongying City	2009 2021	2.6–8275.46 56.25–582.56	1744.41 149.63	This study This study

Note: YRDNR means “Yellow River Delta Natural Reserve”.

). From the perspective of spatial distribution, the PAH concentration in the soil of the northern oil field area of the Yellow River Delta is higher than that in the southern area of the Yellow River, the Yellow River estuary and the former route [29,30]. Among them, four samples exceeded the mean value, accounting for 26.67% of total soil samples, and eleven samples were below the mean value, accounting for 73.33% of total soil samples. The concentration of seven carcinogenic PAHs ranged from ND to 291.28 ng/g, with an average of 28.3 ng/g, accounting for 18.91% of ΣPAH₁₆. The BaP content ranged from ND to 107.88 ng/g, with an average of 13.19 ng/g, respectively, accounting for 46.61% and 8.82% of ΣPAH₇ and ΣPAH₁₆. The other PAH monomers, except for Ace, Acpy, DahA, BkF and BghiP, were detected to varying degrees, and the detection rate was between 7% and 100%. The detection rate of Nap was the highest, followed by BaP. The detection rate of Chr was the lowest (Table 4). In previous reports, Nap and Phe dominated ΣPAH₁₆ in the soil of the YRD [33], and Chr, Phe and Pyr dominated ΣPAH₁₆ in the soil of the YRD [15]. Among the 15 samples, only two samples were slightly polluted, accounting for 13.33% of total soil samples, and the remaining samples were not polluted. Previous research showed that 60.1% of the soil samples were slightly polluted by ΣPAH₁₆, and 18.1% of the soil samples were moderately polluted [15]. Combined with our research results, the pollution situation in the study area was significantly improved compared with 2009 (Figure 2).

The main reasons may be related to the improvement in production conditions, the enhancement of staff’s awareness of environmental protection and the improvements in technology. According to the research, the water source of the YRD was mostly surface runoff, which caused serious soil and water loss [34,35]. Severe soil erosion may also carry away some PAH pollutants from the soil into rivers and oceans [20,36,37]. This may also be part of the reason for the drop in pollution levels. In addition, there were also studies on oil pollution and soil microbial remediation. Two strains selected by previous researchers from oil-contaminated soil can deal with phenanthrene pollution in a salinity environment, among which strained FM6-1 metabolizes phenanthrene through the “phthalic acid” pathway and strained FM8-1 metabolizes phenanthrene through the “naphthalene” pathway [38]. This may also explain a decrease in the concentration of PAHs in the soil near the well in 2021. However, the proportion of seven carcinogenic PAHs and the most carcinogenic BaP were increasing. The management of PAHs in this area was still worthy of attention.

3.2. PAH Component Characteristics

The composition of ΣPAH₁₆ in 2009 and 2021 is shown in Figure 3. In 2009, the soil samples were dominated by low-ring PAHs (2–3 rings). The proportion of each ring in 2009 was 3 rings (88.94%) > 4 rings (5.15%) > 2 rings (4.98%) > 5 rings (0.52%) > 6 rings (0.41%). Among the low-ring PAHs, Flu contributed significantly, accounting for 86.01%. The proportion of each ring in the soil in 2021 was similar to that in 2009, and the low-ring PAHs were dominant. The proportion of each ring in 2021 was 2 rings (75.38%) > 5 rings (12.19%) > 4 rings (5.78%) > 6 rings (3.83%) > 3 rings (2.82%). The difference is that Nap contributed the most to 96.39% of the low-ring PAHs in the soil in 2021, and Flu accounted for only 1.39%. Low-ring PAHs mainly exist in the environment in a non-persistent gaseous form [39,40]. In contrast, high-ring PAHs are relatively resistant to microorganisms and photolysis [41,42]. In addition, the overall decrease in low cyclic PAHs was 15.72%. With the change in time, the proportion of low-ring PAHs gradually decreased [15,39,43]. The proportion of high-ring PAHs gradually increased from 6.08% to 21.8%. It is speculated that this may be related to the rapid economic development in the past ten years. With the increasing number of chemical factory enterprises, vehicles and the increasing intensity of human activities, oil pollution is still the long-term and main source of PAHs in this area.

3.3. PAH Source Apportionment

Now, the main cause of PAH pollution is anthropogenic sources, including oil sources and pyrolysis sources. Specifically, oil source refers to oil spill events that occur in the process of oil exploration and development, transportation and storage and offshore ship operations; pyrolysis source refers to the production activities of chemical enterprises such as petroleum processing, coking gas, steel making and iron making, the incomplete combustion of fuel oil, coal, gas, wood and other biomass in daily life, as well as waste incineration, smoking, cooking and frying, and the emission of various vehicle exhaust [44,45]. In general, an individual PAH ratio is used to identify the source of PAHs in the soil [33,46], different PAH molecular weights or different aromatic ring numbers are used to identify different PAH sources [47,48].

Figure 4 shows the source composition of PAHs in soil in 2009 and 2021. The results show that soil PAHs in 2009 were mainly affected by six sources (Figure 4a). The contribution rate of factor 1 to high-ring PAHs is higher, such as BaA (77.7%), Chr (56.36%) and BbF (35.51%). Gasoline and diesel produce BaA and Chr when they are not completely burned [49]. BbF is a typical substance emitted by diesel vehicles [50], so this factor can be identified as traffic pollution source 1. The main load of factor 2 is Ant, which represents the biomass combustion source [51]. Factor 3 has a large load on 2~3 ring PAHs, including NaP, Flu, Ant and Phe. These low molecular weight components are closely related to oil pollution [52,53,54], so factor 3 is considered to be an oil pollution source. BghiP (87.88%), IcdP (67.49%) and Pyr (48.34%) were found to be the higher-loading compounds in factor 4. BghiP is a typical substance of gasoline combustion emissions [55], so factor 4 is a mixed source of traffic pollution and fossil fuel combustion. The load of factor 5 was mainly distributed on BkF (77.98%) and BaP (76.31%), followed by BbF (37.59%). These high molecular weight PAHs represent traffic pollution source 2. The highest contribution rate of factor 6 to Fla was 66.32%. Phe and Flu were medium loads. So, factor 6 could be considered as the combustion of coal and biomass [56].

In 2021, soil PAHs were mainly affected by four sources (Figure 4b). Factor 1 not only has a higher load on low-ring PAHs such as NaP and Flu but also has a greater contribution to the high-ring components BaA, BaP and BkF. Therefore, factor 1 can be identified as a mixed source of oil pollution and traffic pollution [57,58,59]. Factor 2 is dominated by Pyr load, and the contribution rate reaches 98.92%, which is much higher than other compounds, reflecting the source characteristics of fossil fuel combustion [52,60]. Factor 3 has the highest contribution rate to Fla, which is 92.41%, representing the source of coal combustion [52,53,59]. The main loads of factor 4 are Phe (99.38%) and Ant (98.78%), which can be attributed to the combustion of coal and firewood. Therefore, factor 4 can be identified as the combustion source of coal and biomass [51].

Figure 5 shows the contribution rate of each pollution source of soil PAHs in different periods via PMF model analysis. The positive matrix factorization model (PMF) is recommended by the United States Environmental Protection Agency for pollution. A factorization model for source apportionment was first proposed by Paatero et al. [61]. In recent years, it has been widely used in soil, sediment, atmospheric particles and so on. The model can avoid the negative value of factor loading and scores in the application process and further optimize the data by using uncertainty for the model output results to achieve strong reliability.

Common isomeric ratios such as Fla/(Fla + Pyr), Ant/(Ant + Phe), BaA/(BaA + Chr) and Fla/(Fla + Pyr) were used to analyze the possible sources of PAHs in the soil around oil wells in the Yellow River Delta. According to the ratio analysis of Fla/(Fla + Pyr) and Ant/(Ant + Phe) [62,63], the main source of PAHs in the soil of samples A9, A11, A12, A13 and A15 in 2009 was biomass and coal combustion, and the main source of PAHs in the soil of sample A2 in 2009 was petroleum combustion. The PAHs in the soil of other plots were mainly derived from petroleum (Figure 6a). According to the ratio analysis of BaA/(BaA + Chr) and Fla/(Fla + Pyr) [64,65], the main source of PAHs in the soil of samples A11 and A13 in 2009 and B12 in 2021 was biomass and coal combustion, and the main source of PAHs in the soil of plots A2 in 2009 was petroleum combustion. The main source of PAHs in other soil samples was petroleum pollution (Figure 6b).

On the whole, oil pollution, traffic pollution, coal and biomass combustion are the long-term and main sources of PAHs in this area. This is consistent with previous research reports [12,66]. The source contribution rates of PAHs in soil in 2009 and 2021 are slightly different (Figure 5). In 2009, the contribution rate of oil pollution sources was the highest, reaching 39.01%, followed by coal and biomass combustion sources, and traffic pollution and fossil fuel combustion sources (1.53%) contributed the least. In 2021, oil pollution and traffic pollution were the main sources. The other three pollution sources and their contribution rates were coal combustion (29.73%), coal and biomass combustion (2.4%) and fossil fuel combustion (3.99%). Zhao found that PAHs in the surface sediments of the YRD were mainly derived from oil leakage and fossil fuel combustion and had a moderate risk to the local ecosystem [12]. Xu found that PAHs in sediments were mainly derived from oil leakage and fossil fuel combustion and had low to moderate risks to local ecosystems [66]. Studying the sources of PAHs in different periods can not only effectively prevent and control pollution but also provide a theoretical basis for relevant departments to formulate corresponding laws and regulations.

3.4. Ecological Risk Assessment of PAHs

3.4.1. Toxic Equivalent Factor Method (TEQ)

The ecological risk assessment was performed via the toxic equivalent of PAH monomer relative to BaP (TEQ_BaP) [24]. The TEQ concentrations of ΣPAH₁₆ in soil ranged from 0.003 to 34.76 ng/g in 2009, with an average of 7.49 ng/g. The TEQ of seven carcinogenic PAHs accounted for 73.79% of the total TEQ. The monomer with the largest contribution to the total TEQ was BaP (33.53%), followed by BaA (25.89%). The TEQ concentrations of ΣPAH₁₆ in soil in 2021 ranged from 0.06 to 124.57 ng/g, with an average of 14.71 ng/g, which was about twice the TEQ of 2009. The TEQ of seven carcinogenic PAHs accounted for 99.1% of the total TEQ, an increase of 25.31% from 2009. The carcinogenic risk was also increasing. In general, the TEQ values of different years showed a trend of high-ring PAHs being greater than low-ring PAHs, indicating that high-ring PAHs had high toxicity. Compared with the reference value (33 ng/g), there was basically no potential ecological risk in 2009 and 2021. It is worth noting that although the pollution level in 2009 was much more serious than that in 2021, the ecological risk did not show the same results. On the contrary, the TEQ value in 2021 was larger, which was mainly due to the dominant position of the Flu monomer in soil PAHs in 2009, but its toxic equivalent factor was smaller. In 2021, BaP in the soil was not only the most abundant but also the most toxic, accounting for 89.7% of the total TEQ. Therefore, the TEQ was larger, especially in the O10 plot, which exceeded the target value by about 3.8 times (Figure 7), and may have been affected by PAHs, with potential ecological risks. Based on the analysis of sampling points and human activities, the sampling point is located beside the road. The emissions of motor vehicle exhausts and the lifestyles of the surrounding residents are conducive to the increase in PAH toxicity in the soil. Therefore, the traffic flow should be controlled to reduce the ecological risk. The ecological risk of the O4 plot has been significantly reduced (Figure 7), which may be related to the reduction in the use of fuels such as coal and biomass, which effectively prevents and controls Flu pollution.

3.4.2. Effect Interval Low Median Method (ERL/ERM)

Figure 8 and Figure 9 show the ecological risk assessment results of individual PAH, low-ring PAHs, high-ring PAHs and ΣPAH₁₆ in soil in 2009 and 2021. According to Figure 8, Flu concentrations detected in 12 of the 15 soil samples in 2009 were higher than the ERM value, indicating that 80% of the samples had frequent ecological risks. In addition, the concentrations of Nap in two sites (13%), Ant in one site (7%), Phe in two sites (13%) and Flu in two sites (13%) were between ERM and ERL, indicating that these monomeric PAHs occasionally produced ecological risks in their corresponding sampling points. At the same time, IcdP was also detected at two sampling points, which was harmful to the ecological environment. The ecological risk of other PAHs was small. In 2021, the ecological risk of each monomer PAH in the soil was significantly improved, and only three sampling points (20%) had Nap concentrations between ERM and ERL, with potential ecological risks. The concentrations of other PAHs in all sampling points were lower than ERL, and there was no negative ecological effect.

Compared with the low-ring PAHs and ΣPAH₁₆ in some samples in 2009 that will cause harm to the ecological environment, the concentrations of low-ring PAHs, high-ring PAHs and ΣPAH₁₆ in soil samples in 2021 will not only decrease but also be less than the ERL values, and there is basically no ecological risk (Figure 9). It shows that the ecological environment quality of the study area is developing well, and the environmental protection system and policy formulated on the prevention and control of PAH pollution have achieved initial success. The future protection in the YRD will continue to be aimed at regulating and reducing PAH pollution and promoting ecological protection and high-quality development in the Yellow River basin.

3.5. Health Risk Assessment of PAHs

Domestic and foreign researchers mainly use the incremental lifetime cancer risk model to assess the health risks of PAHs in soil. Halfadji et al. studied the health risks of PAHs in soil in northwestern Algeria and found that adults and children had the highest potential cancer risk and adolescents had the lowest risk [67]; Zheng et al. combined the PMF model with the lifetime cancer risk increment model to analyze the source of cancer risk and found that biomass burning was the main factor causing the carcinogenicity of soil PAHs in the Chengdu Economic Zone [68].

Table 6. ILCR values of PAHs in soils under three exposure pathways in 2009 and 2021.

Exposure Pathways	2009		2021
	Adult	Child	Adult	Child
ILCRingestion	8.56 × 10−9	4.28 × 10−9	1.68 × 10−8	8.4 × 10−9
ILCRinhalation	6.64 × 10−13	4.74 × 10−20	1.3 × 10−12	9.31 × 10−20
ILCRdermal	1.52 × 10−8	5.34 × 10−9	2.99 × 10−8	1.05 × 10−8
ILCRs	2.38 × 10−8	9.62 × 10−9	4.67 × 10−8	1.89 × 10−8

shows the carcinogenic risk values of soil PAHs under three exposure pathways in 2009 and 2021.

In 2009, the carcinogenic risks of soil PAHs to adults and children through three exposure pathways of ingestion of soil, respiratory intake and skin contact were, respectively, 8.25 × 10⁻¹²~1.1 × 10⁻⁷ and 3.34 × 10⁻¹²~4.46×10⁻⁸, and the average values were, respectively, 2.38 × 10⁻⁸ and 9.62 × 10⁻⁹ (Table 6). The carcinogenic risks of PAHs in soil to adults and children in 2021 were 1.78 × 10⁻¹⁰~3.95 × 10⁻⁷ and 7.22 × 10⁻¹¹ and 1.6 × 10⁻⁷, respectively, with an average of 4.67 × 10⁻⁸ and 1.89 × 10⁻⁸, respectively (Table 6). In terms of time, ILCRs in both years were lower than the safety threshold of 10⁻⁶, and there was no carcinogenic risk. Compared with other studies on health risk assessment of soil PAHs in the YRD [14,17], the health risks of soil PAHs in this study were lower in 2009 and 2021 on average (

Table 7. ILCR values of PAHs in soil samples reported in the previous literature.

Study Area	Study Year	Mean of ILCRs (Adult)	Mean of ILCRs (Child)	Reference
Yellow River Delta Natural Reserve	2017	1.14 × 10−6	1.23 × 10−6	[17]
Farmland soil in the YRD	2020	9.00 × 10−6	3.60 × 10−5	[14]
Xianhe Town, Gudao Town and Huanghekou Town of Dongying City in the YRD	2009	2.38 × 10−8	9.62 × 10−9	This Study
Farmland soil in the YRD	2021	4.67 × 10−8	1.89 × 10−8	This Study

). Consistent with the trend of ecological risk, the health risk in 2021 was slightly higher than that in 2009 mainly because the toxic equivalent factor of BaP was large. Despite this, the average ΣPAH₁₆ content in the surface soil of the Yellow River Delta in 2009 was ten times that of the average ΣPAH₁₆ content in 2021. The soil PAH status in 2021 was at the unpolluted level. Therefore, compared with 2009, the soil pollution situation in 2021 is significantly better. It is worth noting that in many organizations, a lifetime cancer risk in the vicinity of 10⁻⁶ and 10⁻⁴ connotes a potential hazard. For instance, Tarafdar et al. (2018) found that fly ash PAH ILCR values of adults and children near oil refineries in India were, respectively, 1.823 × 10⁻⁵ and 1.854 × 10⁻⁵ [69]. The simulated risk exceeded the negligible risk of 10⁻⁶ but fell below the acceptable risk benchmark of 10⁻⁴, which meant that people from both age groups have a plausible threat of cancer occurring [69].

Skin contact and ingestion of soil are the main means of soil PAH exposure, which was consistent with previous studies [13,14,70]. Zheng et al. used the lifetime cancer risk increment (ILCR) to evaluate the polycyclic aromatic hydrocarbons in farmland soil in Ningde. The results showed that the ILCR value ranged from 7.1 × 10⁻⁴ to 1.1 × 10⁻³, which mainly caused moderate to high cancer risk to human health through soil intake and skin contact [13]. Xie found that among the three exposure pathways, skin contact was the most important exposure mode of soil PAHs, followed by oral ingestion, and inhalation exposure was minimal [14]. Studies have shown that long-term exposure to high concentrations of PAHs in the environment through skin contact, oral and nasal intake would cause damage to the respiratory system, nervous system, liver and kidneys, causing lung cancer, nasopharyngeal carcinoma, skin cancer, gastric cancer and other diseases, which would cause great harm to health [71,72,73,74]. Therefore, it is essential to clarify the concentration, content and health risk of ΣPAH₁₆ [43,75] and actively study and improve the pollution status of PAHs. From the perspective of the population, the health risk of soil PAH pollution in adults is higher than that in children, because ILCR_dermal contributes the most to ILCRs, which was inconsistent with previous studies [14,17,67,69]. The specific causes need to be further studied, so caution should be taken when referring to the management measures of the threat of PAHS pollution to adults and children in the YRD.

4. Conclusions

Based on the investigation of the research group in 2009, this study re-collected samples according to the same sampling points and analyzed the concentration, composition, source, ecological risk and health risk of PAHs. The concentration of ΣPAH ₁₆ in the surface soil of the Yellow River Delta in 2009 ranged from 2.6 to 8275.46 ng/g, with an average of 1744.41 ng/g, which reached the level of severe pollution. The concentration of ΣPAH 16 in 2021 ranged from 56.25 to 582.56 ng/g, with an average of 149.63 ng/g, which was unpolluted. The composition and distribution characteristics of PAHs in the soil in 2009 and 2021 were similar, and low-ring PAHs are dominant. The main source of soil PAHs in 2009 was oil pollution. In 2021, soil PAHs were mainly affected by oil pollution and traffic pollution, which further confirmed the analysis results of the proportion of rings. The evaluation results of the toxicity equivalent factor method showed that there was no potential ecological risk in the soil in 2009 and 2021, but the TEQ value in 2021 was large. It was mainly because the Flu monomer in the soil PAHs in 2009 was dominant, but its toxicity equivalent factor was small. The ecological risk of each monomer of PAH also improved significantly. Only 20% of the NaP samples had potential ecological risks, and the remaining PAHs had no negative ecological effects. The evaluation results of the lifetime cancer risk increment model showed that the ILCRs of soil PAHs in 2009 and 2021 were lower than the safety threshold of 10⁻⁶, and there was no carcinogenic risk. From the perspective of exposure pathways, the highest carcinogenic risk was skin contact, followed by ingestion of soil, and the lowest carcinogenic risk was respiratory intake.