Literatures were searched from full-text databases including Web of Science (WOS), Google Scholar, China National Knowledge Infrastructure (CNKI), and Wanfang Database in the last 2 decades (2003–2022), by using the following key words “China”, “PAHs”, “soil”, and specific city names. Data collection and screening followed the following principles: (1) 16 priority controlled PAHs were reported, i.e., naphthalene (Nap), acenaphthylene (Acy), acenaphthene (Ace), fluorene (Flu), phenanthrene (Phe), anthracene (Ant), fluoranthene (Fla), pyrene (Pyr), benzo [a] anthracene (BaA), chrysene (Chr), benzo [b] fluoranthene (BbF), benzo [k] fluoranthene (BkF), benzo [a] pyrene (BaP), dibenzo [a,h] anthracene (DahA), benzo [g,h,i] perylene (BghiP), and indeno [1,2,3-cd] pyrene (IcdP); (2) data that appeared directly in the literature (individual PAH or total PAH concentrations) were used for inclusion in the statistics, and literature that did not give information on sampling time, spatial location of sites, and sampling data for individual and total PAHs content in the search results were discarded; (3) if soil samples were collected from contaminated sites, such as chemical plants, sewage irrigation areas, and oil pollution irrigation areas, they were deleted; (4) for location information of sites in the literature, if only latitude and longitude ranges or all latitude and longitude coordinates were given, the latitude and longitude of the center point of the described range was selected as a representative.

A total of 155 papers was screened, from which PAH concentrations from 245 sampling sites from topsoil at a depth of 0–20 cm were compiled. The specific information about the references, including sampling time, location, size, and characteristics are described in the Supporting Information (Supplementary Table S1). All the soil sampling sites were distributed in 22 provinces (excluding Taiwan Province), 4 municipalities, 5 autonomous regions, and 2 special administrative regions of China, and were divided into seven geographical regions based on location and climate: Northwest China (Shaanxi, Gansu, Qinghai, Ningxia, Xinjiang, and western Inner Mongolia), East China (Jiangsu, Shandong, Fujian, Shanghai, Zhejiang, Anhui, and Jiangxi), Central China (Henan, Hubei, and Hunan), North China (Hebei, Tianjin, Shanxi, Beijing, and central Inner Mongolia), Northeast China (Liaoning, Jilin, and Heilongjiang, Eastern Inner Mongolia), South China (Shenzhen, Hainan, Guangdong, Guangxi, Hong Kong, Macao), and Southwest China (Chongqing, Sichuan, Guizhou, Yunnan, and Tibet) (Figure 1).

In addition, energy and economic data were collected from municipal statistical yearbooks of those provinces from 2003 to 2022. These data included the following variables: coal consumption, coke consumption, gasoline consumption, kerosene consumption, diesel oil consumption, fuel oil consumption, natural gas consumption, electrical power consumption, liquefied petroleum gas consumption (LPG), petroleum consumption, and base oil consumption, the number of car ownership, population and gross domestic product (GDP).

Positive matrix factorization (PMF) is one of the most widely used source resolution methods (Zhang et al., 2019), which is applied to quantify and separate the contributions of pollutant sources in different environmental media, including soil, sediment, and air (Souto-Oliveira et al., 2021; Sakizadeh and Zhang, 2021). PMF mainly quantifies the contribution of different factors to pollutant concentrations by mathematical statistical methods. PMF decomposes the data into two matrices: source contribution (G) and factor spectrum (F), and obtains the results by minimizing the objective function Q using the weighting of the concentration and uncertainty data of the samples, and the formula for Q is as follows: Q = ∑ i = 1 n ∑ j = 1 m x i j − ∑ k = 1 p g i k f k j u i j (1) In Equation 1: where Q is the cumulative residual, and i is the number of samples, and j is the type of component measured; p is the number of factors resolved by the PMF model; f is the component matrix for each source; g is the contribution matrix of each pollutant in the sample; u _ij is the uncertainty of the components in the sample (Sun et al., 2023). EPA PMF 5.0 software was used to analyze the sources of PAHs in the surface soil of seven regions in China, and the default mode was used to load the concentration and uncertainty concentration matrices of PAHs from different regions into EPA PMF 5.0, respectively, and the detection limits of the statistical literature were uniformly organized as 0.8–3 μg kg⁻¹, and the number of factors was determined by finding the smallest Q-value in order to ensure the reasonableness of the analysis. The Q value was minimized and kept stable with R ² > 0.75 when using 4 factors and 20 iterations of the factorial operation were performed.

Geographically weighted regression (GWR) is a combination of spatial correlation and linear regression (Brunsdon et al., 1996), which has certain advantages over traditional models such as the least squares (OLS) method because it supports parameter estimation of local changes in the associations between the independent variables and the dependent variable, so that the relationship between the variables can change with geographic location and can reflect neglected local characteristics (Zhang and Ren, 2012; Ji et al., 2014), and its model is structured as follows (Mcmillen, 2004): y i = β 0 u i , v i + β 1 u i , v i x i 1 + β 2 u i , v i x i 2 + ⋯ + β k u i , v i x i k + ε i (2) in Formula 2: yi represents the explained value of the dependent variable for sample point i; u i , v i represents the geographic coordinates of sample point i ;   x i k ( k = 1 , ⋯ , 12 ) represents the explained value of the independent variable for sample point i; β k u i , v i represents the regression coefficient at sample point u i , v i ; ε i refers to random error term which shows a positive correlation between independent and dependent variables when ε i >0, while negative correlation when ε i < 0 (Wang et al., 2020). The essence of GWR lies in selecting an appropriate spatial weighting function to achieve spatial correlation between data through a spatial weight matrix (Liu et al., 2022). Common methods include the Gauss function method and bisquare function method, etc. This study used the Gauss function method to determine spatial weights.

The Shapiro-Wilk test for lognormal distribution of concentrations of individual PAHs and total PAHs showed significant results (p-value less than 0.05), indicating that none of the study data obeyed a normal distribution, and therefore the median of the concentration was used to represent its mean. Statistical analysis and data mapping were done using IBM SPSS Statistics 26.0 (SPSS Inc., Chicago, United States) and Origin 2022 (Origin Lab Inc., Massachusetts, United States), geostatistical analysis and GWR were done using ArcMap 10.2 (ESRI, Redlands, CA), and the PMF modeling was done using EPA PMF5.0 (EPA, Washington, United States), and the national political zone map was obtained from the website of standard maps of the Ministry of Natural Resources (http://bzdt.ch.mnr.gov.cn/), the original review map number is GS (2020) 4,619, and the base map is not modified.

Concentrations of total PAHs (Σ₁₆PAHs) in China soil ranged from 1.05 to 6,778 μg kg⁻¹ or four orders of magnitude, with the arithmetic mean concentration and median value were 710.49 μg kg⁻¹ and 407.05 μg kg⁻¹, respectively. The coefficient of variation was up to 144.19%, which belongs a strong variation, indicating that the status of PAH pollution varied greatly throughout China. This may have been due to the impacts of anthropogenic activities on PAH occurrence. Differences in economic development, energy consumption, climatic characteristics, and land use could synthetically affect the dissimilarity in PAH pollution (Wang et al., 2015; Peng et al., 2016a; Li et al., 2023). According to Maliszewska-Kordybach’s soil PAH pollution classification criteria, Σ₁₆PAH contamination levels can be divided into 4 categories: not contaminated (<200 μg kg⁻¹), weakly contaminated (200–600 μg kg⁻¹), contaminated (600–1,000 μg kg⁻¹), and heavily contaminated (>1,000 μg kg⁻¹) (Maliszewska-Kordybach, 1996). PAH pollution levels across China soil varied with different regions, about 17.45% of the soil were considered heavily contaminated, 16.04% and 40.57% were classified as contaminated and weakly contaminated, respectively, and only 25.94% were not contaminated.

Compared with the studies in other countries, the mean Σ₁₆PAH concentrations in China soil were generally higher than that reported in Italian (84–333 μg kg⁻¹, Thiombane et al., 2019), Japan (243 μg kg⁻¹, Hong et al., 2020), and rural soils in the United States (138 μg kg⁻¹, Yuan et al., 2014), and were comparable to those of surface soil in Korea (960 μg kg⁻¹, Kwon and Choi, 2014). While, mean Σ₁₆PAHs in China soil were significantly lower than those reported in Iran (1,241 μg kg⁻¹, Ashjar et al., 2022), and Germany (1,448 μg kg⁻¹, Aichner et al., 2015). The above comparison results showed that most of the PAH contamination levels in China soil were considered moderate, but soil PAH pollution may pose potential environmental and human health risks.

Obvious differences in the concentrations and ratios of PAHs with different ring numbers were observed (Figure 2). The 4-ring PAHs were generally the most abundant composition with the median concentrations of 163.04 μg kg⁻¹, accounting for 38.55% of the total PAHs. In particular, Phe and Fla exhibited the highest contents, with median values of 42.9 μg kg⁻¹ and 39 μg kg⁻¹, respectively. This was consistent with the previous national-scale investigation on surface soil of China (Ma et al., 2015). Subsequently, 3-ring (median concentration of 106.03 μg kg⁻¹) and 5-ring (median concentration of 80.76 μg kg⁻¹) PAHs followed, constituting 25.07% and 19.09% of the total concentration, respectively. The 2-ring PAHs recorded the lowest concentration (median concentration of 29.62 μg kg⁻¹), comprising 7% of the total amount. Ace and Acy exhibited the lowest monomer concentrations, with median concentrations of 3.54 μg kg⁻¹ and 3.85 μg kg⁻¹, respectively, corresponding to the results of Cao et al. (2012). Therefore, it is concluded that the relative proportion of 3 or 4-ring soils PAHs in China is high, primarily originating from combustion sources (Li and Duan, 2015). Additionally, traffic sources, emitting PAHs with 5 or 6-rings, also make a major contribution (Zhang and Chen, 2017).

The composition of PAHs varied across different regions (Figure 3), leading to different contamination levels and exhibiting diverse transport and accumulation behaviors in soil (Yu et al., 2019). Particulate PAHs with 5 or 6- rings, owing to their low saturation vapor pressure values, exhibited lower transportability, higher toxicity, aggregation, and degradation resistance, exhibiting a significant increasing trend over time. Conversely, PAHs with 2–3 rings exhibited greater volatility and were photodegradable in soil, with their concentration decreasing over time (Jones et al., 1989; Wild and Jones, 1995; Wilcke et al., 2005). Figure 3 depicts the ratios of PAHs with different ring numbers in soil across different regions of China. Notably, 4-ring PAHs dominated in all regions, accounting for 29.26%–46.5%. Studies have shown that 4 rings PAHs are mainly derived from coal combustion sources (Li and Duan, 2015), indicating a prevalence of coal-fired gas as the source of PAHs in China soils. Furthermore, comparing low-ring and high-ring PAHs in each region using 4-rings as the boundary, North China, Southwest China, South China, and Central China exhibited a higher proportion of PAHs with 5 or 6-rings, accounting for 36.26%, 34.86%, 34.44%, and 33.48%, respectively. Automobile exhaust primarily releases 5-6 rings PAHs (Zhang and Chen, 2017), indicating traffic sources as the primary contributors to soil PAHs. The concentrations of PAHs with 5 or 6-rings in North China were significantly higher than in other regions, corresponding to the results of Zhang Y. et al. (2019). In Northwest China and East China, PAHs with 2–3 predominated, constituting 48.92% and 41.76%, respectively, indicating petroleum spills, petrochemicals, incomplete combustion of biomass, and coke, as the main sources (Zhang Y. et al., 2019). Given the higher frequency and concentrations of PAHs with 5 or 6-rings (Yu et al., 2019), effective prevention and control measures are necessary to manage PAHs pollution in North China.

The model validation of the semi-variance function identified the best-fitting and most accurate theoretical model for interpolation, and the results are shown in Table 1. The ratio of the nugget (C₀) and the Sill (C₀ + C) is regarded as the nugget coefficient (C₀/C₀ + C). This coefficient is indicative of the spatial autocorrelation of PAH contents in China. When the nugget coefficient <25%, the variable exhibited strong spatial autocorrelation, and the variations were attributable to structural factors such as topography, precipitation, climate, parent material type, and other non-human factors. When the nugget coefficient was within the range of 25%–75%, the variable exhibited moderate spatial autocorrelation. Conversely, a coefficient >75% indicated weak spatial autocorrelation, with variations primarily arising from random factors such as releases from anthropogenic sources such as coal combustion, oil leakage, automobile exhaust, and more (Zhang et al., 2014). In Table 1, the nugget coefficient of Σ₁₆PAHs and 3-ring PAHs were less than 25%. PAHs were emitted from the surface into the atmosphere, and the spatial heterogeneity of PAHs has been weakened during atmospheric transport. Therefore, Σ₁₆PAHs and 3-ring soil PAHs exhibited strong spatial autocorrelation, primarily influenced by structural factors and less affected by human activities. The nugget coefficients of 2-ring PAHs and 4-6 ring PAHs were all higher than 25%, indicated that both 2-ring PAHs and 4-6 ring PAHs exhibited a moderate level of spatial autocorrelation. The spatial heterogeneity of PAH concentration was attributable to a combination of structural and random factors, with random factors inferred as the primary factors leading to the spatial heterogeneity of the soil PAHs. Therefore, the spatial autocorrelation of soil PAH concentrations in different regions was affected by long-term anthropogenic factors such as coal-fired gas, petrochemicals, and other forms of energy consumption.

The Kriging interpolation method in ArcGIS 10.2 geostatistical analysis was employed to interpolate the collected Σ₁₆PAH and different rings PAH concentrations. Given the spatial interpolation results (Figure 4), the spatial distribution of soil PAHs exhibited characteristic aggregation. The spatial distributions of mean Σ₁₆PAH concentrations descended in the order of Northeast China > North China > Central China > East China > Northwest China > Southwest China > South China, with different PAH composition profiles among regions affected by energy consumption. The northeastern region exhibited the highest Σ₁₆PAH concentration, which was an important old industrial base in China, corresponding to the results of Ma et al. (2015). There were several areas of high concentration, with high latitude, low annual average solar radiation, and a large amount of coal consumed for heating in winter, resulting in a low amount of radicalized PAHs. Distributions of PAH concentrations with 3–6 rings in Northeast China exhibited higher interpolation results. PAHs with 3 or 4-rings were typically thought to come from the incomplete combustion of coke, coal, and biomass (Li and Duan, 2015), while PAHs with 5 or 6-rings were evident in developed transportation areas, primarily from vehicle emissions (Zhang and Chen, 2017). Therefore, the massive combustion of coal, the production of the petroleum industry, and the heavy traffic pressure have led to serious pollution of the soil in Northeast China (Sun et al., 2018).

As an important economic hotspot in China, North China is also the location of the Beijing–Tianjin–Hebei industrial region. Distributions of PAH concentrations with 2–6 rings in North China exhibited a high interpolation value, and the 2-ring PAHs in the environment corresponded to petroleum pollution (Li and Duan, 2015). Therefore, the soil PAHs in North China were primarily derived from the combustion and the traffic sources, and the economic development exhibited increased energy consumption (Ma et al., 2015), leading to more serious soil contamination influenced by human factors.

The spatial distribution of 3–6 ring PAH concentrations in Central China was characterized by high concentrations, primarily in Jiangxi Province, corresponding to the results of Zeng et al. (2019). The 3-ring PAHs were primarily emitted from coke-oven combustion and petrochemical factories (Wang et al., 2016). Thus, it was inferred that the main sources of PAHs in the soils of Central China were industrial emissions and vehicle exhausts.

There were high concentration areas in the spatial distribution of the PAH concentrations with 3–4 rings and 6-rings in East China, primarily concentrated in the Yangtze River Delta region. The economy of East China is developed, and the per capita gross domestic product (GDP) of most of the provinces was 2–4 times higher than that of the western provinces. Thus, the difference in the regional economic pattern was attributable to the higher PAH concentrations in East China than in West China (Zhang and Chen, 2017). The Shanghai–Nanjing–Hangzhou industrial region in East China emitted more pollutants, resulting in a higher interpolation result for the 3-ring PAHs. The Yangtze River Delta, which was in the transition zone between North and South, also accounted for a considerable proportion of coal-fired gas, emitting the 4-ring PAHs. The 6-ring PAHs were derived from diesel combustion in East China, inferring that the main sources of PAHs in the soils of East China were combustion and traffic sources.

The spatial interpolation of distributions of PAH concentrations with 2–4 rings in Northwest China was relatively high. Xinjiang had the highest base reserves (NBS, 2016). The volatile pollution of PAHs with 2 or 3-rings existed in the process of oil exploitation and transportation. Urumqi had a developed petrochemical and metallurgical industrial base, which increased the emission of PAHs by 3 or 4-rings. Additionally, the altitude and latitude of Xinjiang are higher, with a lower annual average temperature, reducing the volatilization and degradation of PAHs in the soil (Zhang et al., 2019).

In Southwest China, the spatially interpolated high levels of 2–4 rings PAHs were attributable to petroleum and petrochemical pollution, coal-fired gas, and incomplete combustion of biomass. Nevertheless, the median value of Σ₁₆PAHs data in Southwest China was low, while the interpolation result was relatively high due to fewer sample points in this area, with most located in industrial areas and basins where air circulation was inconvenient, leading to regional contamination. The Pearl River Delta (PRD) in South China was characterized by the development of light industry, exhibiting low emissions of PAHs (Zhang et al., 2019).

The distribution of PAHs in South China not only affected the level of economic development and energy consumption but also the climatic conditions. South China is located at a low latitude, high temperature, and significant ultraviolet radiation and volatilization, photolysis, and microbial degradation of PAHs (Zhang and Chen, 2017), resulting in the lowest Σ₁₆PAH concentrations in the soil, corresponding to the results of Yu et al. (2019). However, the concentration of PAHs with 3–4 rings and 6-rings in South China exhibited a high spatial interpolation, suggesting that PAHs in South China were primarily derived from biomass, incomplete combustion of coal and natural gas, and diesel combustion.

To explore the possible sources of PAHs in seven major regions’ soils of China, the PMF model was first used to apportion the potential source based on collected PAH concentrations. The specific information related to PFM modeling are described in the Supporting Information (Supplementary Figure S1; Supplementary Table S2). In the northwest region, Factor 1 was dominated by Flu and Phe with over 70% contributions. Flu and Phe were dominant release substances from the steel industry and coke production (Ravindra et al., 2008). Therefore, factor 1 revealed the origin from industry production. Factor 2 was dominated by Ant and most 4-ring PAHs, including Fla, Pyr, BaA, which was the indicator of coal combustion (Li and Duan, 2015). Factor 3 was dominated by HMW PAHs, which was the indicator of gasoline and diesel combustion (Khalili et al., 1995; Simcik et al., 1999; Yang et al., 2021). Factor 4 had significant correlations with Nap, which mainly originated from petroleum spill (Hu et al., 2014). In addition to extracting the dominant sources, the percentage contributions to PAH source were also determined from PMF model (Figure 5). The percentage contributions of the four factors are21.47%, 33.72%, 14.13% and 30.69%. Therefore, the primary PAH source in northwest was coal combustion, followed by oil spill pollution. Similarly, four factors were also extracted for the other six major regions. In the East region, the percentage contributions of factor 1 is 29.28%, factor 1 was dominated by Chr, BbF and Bap, which was the indicator of gasoline combustion (Wang et al., 2008; Xing et al., 2011). The main contribution factor in the central region is a mixture of origins from industrial emission and diesel combustion with 3-ring PAHs, Bap and IcdP as dominant release substance, accounting for 33.09%. The factor with the highest contribution rate in the north region is the coal and natural gas combustion source, which was dominated by the 4-ring PAHs. The highest contribution factor in northeast region is factor 4 dominated by Fla, Pyr, BkF, BaP, BghiP and IcdP. Fla and Pyr were dominant release substances from coal burning, while BkF, BaP, BghiP and IcdP are mainly emitted from automobile exhaust. Therefore, soil PAHs in the northeast region are mainly caused by coal combustion and traffic pollution. The dominant factor in the south region is biomass combustion source with Nap and Phe as the main load. PAHs in southwest region are mainly derived from coal combustion sources and transportation sources with MMW and HMW as the main contribution.

PAHs present in soil surfaces across Northeast China, North China and Northwest China are primarily derived from coal combustion, with coal-fired heating in the northern region emitting a significant amount of PAHs in winter (Zhao et al., 2008; Liu et al., 2016). Additionally, Northeast China and North China, acting as major industrial and economic hubs in China, faced heightened challenges due to the production demands of the petroleum industry, ongoing modernization of transportation networks, and substantial energy consumption. These factors collectively contributed to a significant increase in soil pollution (Gao et al., 2012; Wang XT. et al., 2015; Peng et al., 2016b; Zhang and Chen, 2017; Sun et al., 2018). In economically vibrant regions like East China and Central China, the primary sources of PAHs in surface soil were traffic-related, particularly emanating from vehicle emissions, indicating the role of economic development and urbanization in increasing pollution levels (Yu et al., 2019). Conversely, Southwest China and South China exhibited PAHs sourced primarily from biomass combustion, indicating a lower level of industrial pollution. Instead, PAHs in these regions primarily originated from household activities such as cooking, corresponding to the results of Yu et al. (2019). In the Northwest China, besides coal combustion, another major source of soil PAHs was petroleum spill. The Northwest China possesses abundant oil reserves, with petroleum chemical plants also established in major cities (NBS, 2016). Therefore, soil pollution was easily caused during the processes of oil extraction, processing, and transportation. Consequently, PAHs in China surface soil predominantly emanated from combustion and traffic sources, with a minor contribution from petroleum and petrochemical sources, which was consistent with previous studies.

Most PAHs released in the environment related to production activities and energy utilization (Gong et al., 2022); therefore, Pearson’s correlation and hierarchical cluster analysis were performed to determine the relationship between PAH pollution and anthropogenic activities (Figure 6). Three major clusters were identified based on PAH concentrations and proportions, as well as socio-economic parameters. On the whole, total PAHs, 2-ring and 3-ring PAHs, natural gas, and electricity were grouped in a group. Economic indicators (population and gross domestic product) and energy consumption indicators (coal, coke, gasoline, diesel, liquefied petroleum gas, petroleum, fuel oil, and crude oil), the concentrations and proportions of 4–6 rings PAHs, kerosene, and car volume, were well clustered in the same category, respectively, which clearly indicating that the close relationship between energy utilization and middle and high molecular weight PAHs. It should be noted that due to the influence of climate, geographical location, energy structure and utilization mode, there is no significant correlation between PAHs in different regions and specific source indicators on the national scale. Therefore, it is necessary to analyze the relationship between PAHs and energy indicators at regional scale.

Geographically weighted regression model was employed to demonstrate the spatial heterogeneity of PAH concentrations in surface soil in China under the influence of energy consumption structure in different regions. Visualization of the regression coefficients corresponding to energy consumption indexes in ArcGIS allowed exploration of the spatial heterogeneity of combined explanatory factors affecting PAH concentrations in surface soil. Explanatory variables include the consumption of coal, coke, gasoline, kerosene, diesel, fuel oil, natural gas, electricity, liquefied gas, petroleum, and crude oil. To ensure significant estimation results, a covariance test for these above indicators was conducted. Among them, diesel, electricity, and crude oil consumption exhibited multicollinearity with other dependent variables and were therefore excluded. The remaining indicators exhibited covariance expansion factor values <10, with conditional indices of <30, indicating the absence of covariance relationships among the analyzed indicators (Table 2).

In the current study, both the multiple linear regression model and the GWR model were employed to analyze the energy consumption and the concentrations of Σ₁₆PAHs in surface soil. The overall prediction effect of the multiple linear regression model was relatively satisfactory, the interpretation (regression R ²) was 38.2%, and the indicators with significant results and positive correlation coefficients were coal, natural gas, and petroleum consumption, indicating a close positive correlation with the contents of Σ₁₆PAHs in the surface soil. The regression coefficients for the three indicators on Σ₁₆PAH contents were in descending order: natural gas, oil, and coal consumption. Regression analyses of the three specific types of energy consumption with the PAH concentrations of different ring numbers showed that coal consumption had positive correlations with the Σ₁₆PAH concentration and 2–4 rings PAHs, and negative correlations with 5–6 rings PAHs. Natural gas consumption was positively correlated with the Σ₁₆PAH concentration and 2–4 ring PAHs, and negatively correlated with the 5–6 rings PAHs. Petroleum consumption has a positive correlation with the concentrations of 4–6 rings PAHs and a negative correlation with the others. The GWR model results were consistent with those of the multiple linear regression model, passing significance tests (p < 0.05). Examined indicators were coal, natural gas, and petroleum consumption, with coal and natural gas consumption exhibiting a significant positive correlation with the content of 4-ring PAHs. Petroleum consumption and the Σ₁₆PAH content exhibited a significant positive correlation effect. Furthermore, petroleum consumption positively correlated with the PAH contents with 5 or 6-rings. GWR model regression coefficient spatial distribution is shown in Figure 7.

The impact of coal energy consumption on the regional contents of the 4-ring PAHs in the surface soil exhibited significant differences between the northern and the southern regions. Notably, the regression coefficient was significantly higher in the northern region than in the southern region, exhibiting a positive correlation effect (Figure 7A). The distribution of 4-ring PAHs in surface soil in China indicated higher levels in North China, attributed to central heating, large coal consumption, and a focus on heavy industry (Liu et al., 2016). The old industrial base of the northeast, the location of the Beijing–Tianjin–Hebei industrial zone, and their large-scale coal-fired industrial emissions contributed to PAH deposition and accumulation (Gao et al., 2012). Coal was the primary source of the 4-ring PAHs (Li and Duan, 2015), with the northern region characterized by high latitude, low temperatures, and less annual solar radiation. This region is mostly inland, which is not conducive to the volatilization, degradation, and water vapor transport of PAHs in soil. Additionally, the terrain is relatively flat, providing optimal diffusion conditions for soil pollution. As a result, PAHs were more prone to deposition and accumulation in this area, particularly in most provinces in the northeast, north, and northwest regions.

There was a strong positive correlation between natural gas energy consumption and 4-ring PAHs emissions, particularly in Southwest China, South China, and Xinjiang, China, with a decreasing trend to the northeast (Figure 7B). China is the third largest natural gas consumer after the United States and Russia. In 2022, total natural gas consumption in China was 372.6 billion cubic meters (NBS, 2022), contributing to PAHs pollution. Furthermore, based on statistics reported in 2022, the top five provinces with the highest natural gas consumption in China were: Guangdong, Jiangsu, Sichuan, Shandong, and Beijing (NBS, 2022), accounting for 9.77%, 8.42%, 7.19%, 6.35%, and 5.83% of the total national consumption, respectively. From Figure 7B, it was observed that the trend of natural gas consumption in each province corresponded to the GWR model correlation coefficient trend, indicating regional variations in the impact of natural gas on PAH pollution.

Fossil energy consumption significantly impacted ΣPAHs in most provinces north of the Huaihe River in the Qinling Mountains, gradually decreasing southward (Figure 7C). The emissions of PAHs with 2 or 3-rings were primarily attributable to petroleum leakage and petrochemical sources in fossil energy consumption. Conversely, PAHs with 5 or 6-rings stemmed predominantly from traffic exhaust emissions, that is, gasoline and diesel combustion in fossil energy (Li and Duan, 2015). Notably, variations existed in the spatial distribution of PAHs with different ring numbers, impacting the fitting accuracy of the GWR model. To further explore the effect of fossil energy consumption type on PAH concentrations in surface soil, a GWR analysis was conducted for fossil energy consumption and the concentrations of PAHs with 2–3, 4, and 5–6 rings. The results indicated that only the contents of PAHs with 5–6 rings were significantly affected (p < 0.05), exhibiting a positive correlation. This indicated that the quantity of traffic fuel in fossil energy consumption promoted the emission of PAHs with 5–6 rings, predominantly in Northeast China and North China (Figure 7D). These regions experienced high latitudes and low winter temperatures, with most vehicles relying on fuel as the primary driving force (Glaser et al., 2005). The challenge of popularizing clean energy in transportation contributed to traffic emissions as a main source of PAHs in surface soil. Conversely, the southern region, where vehicles predominantly run on electricity, exhibits a lower significance, attributed to the widespread use of clean energy (Zhang Y. et al., 2019).

In summary, Σ₁₆PAH concentration from the perspective of spatial distribution, Northeast and North China showed the largest pollution degree, with several high concentration areas, and this trend was consistent with the spatial distribution of the regression model for coal and petroleum energy consumption, indicating that the combustion of coal and petroleum is an important influencing factor affecting the spatial pattern of Σ₁₆PAHs in China’s topsoil. Meanwhile, the Σ₁₆PAH concentration in the southwestern region also shows higher results, which is consistent with the spatial distribution of the regression model for natural gas energy consumption, indicating that natural gas is also one of the important factors contributing to the pollution of ΣPAHs. In terms of spatial distribution of 2–3 rings PAH content, there are large areas of regional pollution in Northwest, Northeast and Central China, which is consistent with the trend of the regression model of petroleum consumption, indicating that petroleum consumption increases the emission of low-ring PAHs, especially oil spills and petrochemical sources mainly cause pollution of low-ring PAHs. The spatial distribution of 3–4 rings PAH content is consistent with that of the spatial regression model of coal consumption, both of which are high in the north and low in the south, and coal combustion is the main source of 3–4 rings PAHs. 5–6 rings PAHs distribution is obviously similar to that of the regression model of petroleum consumption, but unlike the low-ring PAHs, the type of petroleum consumption that causes the pollution of high-ring PAHs mainly originates from gasoline and diesel combustion, and the emission of high-ring PAHs by automobile emissions. In addition, the regression spatial model of natural gas consumption shows that natural gas combustion mainly affects the PAH content in Southwest China and Xinjiang province. The spatial distribution of 3 rings PAHs in Southwest China and Xinjiang province shows a higher level in agreement with the regression model, which suggests that natural gas consumption affects the spatial pattern of 3 rings PAHs. China should focus on managing coal and oil consumption in the northern region and natural gas consumption in the southern region in terms of environmental pollution control. The impact of energy consumption on the concentrations of PAHs in surface soil in China was nationally significant, revealing apparent spatial heterogeneity primarily linked to coal, natural gas, and fossil energy consumption. This indicated that combustion and traffic sources were the primary contributors to PAHs in surface soil in China, with coal, natural gas combustion, and oil consumption identified as the principal energy factors influencing PAH concentrations.