The modern coal chemical industry has been developed rapidly with the increasing energy demand due to the economy development (Xia and Han 2017). However, it has drawn great attention that the coal chemical industry uses tremendous energy and water, and produces lots of pollution (Liu et al. 2017). Heavy metals metal (loid)s (HMs) are major pollutants in the coal chemical industry, and they can enter the environment with wastewater, gas emission, and solid residuals (Wang 2012). Soil is the main media of HM deposit, and the accumulated HM in the soils may pose human health risks through chronic exposures (Wang et al. 2012). Previous studies reported that HMs can be toxic, carcinogenic and teratogenic, and about 90% of human carcinogens were caused by the uptake of or exposure to high concentrations of HMs and other chemicals for a long period of time (Yi et al. 2011; An et al. 2016). China’s Prevention and Control of Heavy Metal Pollution in the 12th Five-Year Plan claimed that Pb, Hg, Cd, Cr and As were classified as main pollutants in China. Therefore, it is of great significance to know the concentrations of the five HMs and their changing tendencies in coal chemical plants and make the corresponding preventive measures. Particularly, Cr mainly exists as Cr(III) and Cr(VI) in the soil (Stewart et al. 2003), and Cr(VI) can cause worse adverse effects on humans than Cr(III) (Wu et al. 2016; Liang et al. 2000). Therefore, the present study mainly focused on Cr(VI).

The commonly used models to predict HM concentrations can be roughly divided into the flux prediction model, the two-stage cumulative rate estimation model, the spatial distribution prediction model at different times and the mathematical prediction model based on Hazen probability curve (Li et al. 2014). The soil environmental capacity model, which is one of the commonly used models for predicting HM accumulation in the soil, uses the material balance principle to deduce the kinetic equation of the soil pollution process (Wang et al. 1993). Zhao et al. (1999) analyzed the parameters of the model and applied the model to investigate the changing tendencies of HMs under different conditions. Lin (2009) used the soil environmental capacity model to forecast the HM concentrations in the soil of a developed industrial city in the south of Guangdong province and obtained the trends of the HM concentrations in the city under the optimistic, normal, and pessimistic scenarios. Chen (2014) used the model to predict the HM concentrations in the soil and obtained the HM spatial distributions in the study area in 2020. However, current studies mainly aimed at urban or agricultural soils and other large sites, few reported heavy metal prediction in a single industrial area, especially in the coal chemical industry.

In the present study, based on the soil environmental capacity model, the Kriging method was adopted to construct spatiotemporal graphs of the HM distributions in the soil and determined the changing tendencies of HM concentrations at different production units. The study can provide the theoretical and practical basis for the prevention and control of the HM pollution in coal chemical plants.

2.1 Research area profile

The coal chemical plant is located in the northeastern part of Ningxia Autonomous Region (Ningxia) (Fig. 1a). The plant trial operation was conducted at the end of 2010, and the plant was officially put into operation at the end of 2011. The terrain of the plant location is relatively smooth and flat with a slope less than 10°. The soil type in the study site is mainly sierozem (China National Environmental Monitoring Centre, 1990). The plant is in the hinterland of the mainland and has the typical continental monsoon climate. The NW wind is prevalent during the whole year with a frequency of 25.8%, followed by the SSE wind with a frequency of 25.5%. The NW wind speed is significantly higher than the SSE wind, thus the NW wind has greater transportation capacity than the SSE wind.

The study site: a Geographical location. b Remote sensing aerial map of the plant. c Schematic diagram of the plant layout

A top view of the study area was obtained by the remote sensing technique (Fig. 1b).The top view can be divided into ten functional units (Fig. 1c). During the operation, the coal is added in the coal corridor from the southeast corner and then transported into the power unit and the gasification unit. Separate slag dumps are in service for the power unit and gasification unit in the north. Three chimneys are in operation at the southeast corner of the power unit. The product is mainly transported by railway and a small amount by the highway.

2.2 Soil sampling

Soil samples were collected in July 2014. The chessboard method (Liu 1987) was used to sample the soil with the spacing of 50 m × 50 m, due to the gentle terrain of the study area and the difference in the concentration of soil pollutants in production units. Sampling points were offset or skipped when they were located in the pools, buildings, slag dumps and other locations inconvenient for sampling. The coordinates of sampling points were recorded using the global positioning system (GPS) (Fig. 2a, b). Topsoil samples (0–20 cm) were collected and air-dried at room temperature. After removing rocks and foreign debris, all samples were ground and passed through a 0.15 mm nylon sieve. The prepared soil samples were sealed in sampling bottles for later analysis.

Distribution of sampling points: the factory layout

2.3 The leaching experiment

A leaching experiment was conducted to obtain the residual rates of the HMs in the soil, which provided important parameters for the soil environmental capacity model in the later section. The pH of the rainwater of the target area is between 6.37 and 7.67, according to the weather bureau of Ningxia Province. Therefore, acid leaching experiment was not involved in this study.

The topsoil from sampling site 1 was used as the uncontaminated reference sample in the leaching experiment since this site was far away from the main production unit. Three Plexiglas leaching columns, with a length of 35 cm, a diameter of 4 cm, and a soil loading of about 0.8 kg, were used. The procedure was as follows: (1) Lead nitrate, mercuric chloride, cadmium chloride, potassium dichromate and arsenic oxide were added to 0.3 kg reference soil at the concentrations of 10, 5, 5, 50 and 10 mg/kg soil for lead (Pb), mercury (Hg), cadmium (Cd), chromium (VI) (Cr(VI)) and arsenic (As), respectively (Shangguan et al. 2015). The homogenized mixture was allowed to stand for 10 days and then used as the source of pollution. (2) The bottom of each leaching column was coated with a layer of absorbent cotton. Each column was packed with 0.5 kg uncontaminated soil, and 0.3 kg of the synthetic polluted soil obtained in the previous step was added. Finally, a layer of absorbent cotton was added on the top of the leaching column. (3) A peristaltic pump was connected to the column, and the down-flow rate was set at 0.4 ml/min, with the experimental period of 15 days when all the water was estimated to flow out from the column. (4) Three soil samples (2 g) were taken out of each column after well mixed at the end of leaching. The HM concentrations in the soil were measured. The soil residual rates of HMs were calculated as (Shangguan et al. 2015):

where K is the soil residual rate (%), C₁ is the HM concentration after leaching (mg/kg), C₀ is the HM concentration before leaching (mg/kg), M₀ is the mass of the soil (kg), and M₁ is the mass of HM addition (mg).

2.4 Sample analysis

The prepared soil samples were digested with HF, HNO_3, and H₂O₂ in a microwave digestion system. The concentrations of Pb, Cd, Cr, and As in the soil samples were then measured by inductively coupled plasma mass spectrometry (ICP-MS, Varian X series II). The Hg concentration was measured using a DMA-80 mercury analyzer (Milestone S.r.L).

2.5 The soil environmental capacity model

The soil environmental capacity model is the mathematical expression of the HM concentration change in soil environmental system (Xia 1989). The core formula is:

where n is the target time period, Q₀ is the initial soil concentration (mg/kg), Q_n is the soil concentration in the nth year (mg/kg), K is the soil residual rate (%), Q is the annual input (mg/kg), and Z is the annual output (mg/kg).

The formula (2) can be summarized as follows: the amount of soil concentration in the nth year = the initial residue + n-year total amount of input − the n-year total amount of output. The following formula is obtained by merging the input and output:

where n is the target time period, Q_n is the soil concentration in the nth year (mg/kg), B is the soil background concentration (mg/kg), K is the soil residual rate (%), and R is the net annual input (mg/kg).

Based on assumptions that the plant annual operation is stable, the annual discharge of pollutants (i.e. the annual input R) is constant, and soil properties in the plant do not change significantly (i.e., K is constant), the annual input R at the sampling point is:

After the R at each sampling point was obtained (with the measured value in 2014 and the background value) and assumed to be constant, the Q_n value at the nth year can be calculated by the formula (3). However, long-term forecasts are of little practical significance to industrial plants due to the variation in the coal quality and production process. Therefore, n was given 10, 15, and 20 (i.e., the year 2021, 2026, and 2031), with the year 2011 as the offset.

2.6 Kriging method

The Kriging method, also known as the spatial local interpolation method, is based on spatial autocorrelation. The unknown sampling points of regional variables were estimated with linear, unbiased, optimal estimation of the interpolation method using the original data of regional variables and the structure of the semivariogram function (Tang and Yang 2006).The core formula is expressed as:

where Z (x₀) is the estimated point concentration (mg/kg), n is the total number of sampling points, Z (x_i) is the concentration of the i point (mg/kg), and λ_i is a set of weight coefficients. The core of the algorithm is to calculate the minimum unbiased estimate of Z (x₀).

More details about the principal and process of this method were shown in a previous study (Zhang et al. 2018).

2.7 Software platform

SPSS19.0 software was used for descriptive statistical analysis and multiple stepwise regression analysis. GS+7.0 was used for semi-variance analysis, interactive semi-variance analysis, and Kriging interpolation. ArcGIS10.0 was used for space autocorrelation analysis, map digitalization, mapping and cropping of interpolation results, and cross-validation.

3.1 Analysis of heavy metal concentration in the soil

The HM background values were lower than those of Ningxia before the plant was put into operation (Table 1). The average concentrations of Hg, Cd(VI) and As exceeded their Ningxia background values after plant operation for 3.5 years. The averages of Hg and Cd (0.092 and 0.553 mg/kg, respectively) were more than four times of the Ningxia background values (0.021 and 0.112 mg/kg, respectively). Thus, the HMs, except for Pb, were enriched in the coal chemical plant. In addition, Pb, Cr(VI) and As had weak variation (coefficient of variation < 0.1), while Hg and Cd had medium variation (0.1 < coefficient of variation < 1) (Zhu 1999). Therefore, Pb, Cr and As had low degrees of dispersion in the soil, indicating that they had relatively uniform distributions. Hg and Cd, particularly Hg, showed a higher degree of dispersion, indicating that the pollution areas of the two elements were more concentrated.

3.2 Longitudinal migration characteristics of HMs in the soil

The results of leaching experiments and the five HM residual rates in the soil are shown in Table 2. The residual rates of Pb, Hg, Cd and As in the soil were all above 90%, in the order as Pb > As > Hg > Cd. However, the residual rate of Cr(VI) was less than 50%, which was mainly due to its strong dissolution morphology and the strong migration capacity (Yang et al. 2014).

3.3 Temporal changes of HM concentrations in the soil

The annual inputs (R) at the 153 sampling points were calculated by the derived formula (4) with the measured concentrations of the samples and the local background values (Table 3). The input (R) of Cr(VI) had a weak variation, and the inputs (R) for the other four HMs had a medium variation, which indicated that the five HM emissions varied greatly in different process units of the plant. The annual input of Cr(VI) was over 10 mg/kg, which was much higher than that of the other elements. Therefore, the Cr(VI) concentration in the soil increased rapidly at the initial operation stage of the plant. However, the Cr(VI) concentration in the topsoil would be quickly stabilized due to its low residual rate in the soil. At the same time, a large quantity of Cr(VI) could enter the deep soil, causing a threat to groundwater.

The HM concentrations in the year of 2021, 2026 and 2031 were estimated with the soil environmental capacity model (Table 4). The concentrations of Hg, Cd, Cr(VI) and As were all below the corresponding Ningxia background values in 2011, but over 90% of them would exceed their background values in 10 years (Table 4), indicating that HMs were enriched quickly in the soil by the coal chemical production activities.

The enrichment rate of Pb was the lowest among the five HMs. The 10-year (from 2011 to 2021) enrichment rate was predicted to be 75.0%, and the Pb concentration in more than 80% of the sampling sites would exceed the Ningxia background value. However, the Pb concentration in 98% of the sampling sites would reach 29.028 mg/kg on average in 2031, and the highest concentration would be about two times of the Ningxia background value. On the whole, the enrichment rate of Pb in the soil is relatively low, but Pb has a great enrichment tendency in the topsoil due to the high residual rate.

The concentrations of Hg and Cd was estimated to be increasing rapidly, and the average values in 2021 were predicted to be significantly higher than their corresponding Ningxia background values. The Hg concentration would increase by about 28.7 times from 2011 to 2021, and Cd concentration would increase by nearly 10 times. The enrichment rates of the two elements were 29.8% and 23.5%, respectively, from 2021 to 2026. However, the enrichment rates from 2026 to 2031 were predicted to be decreased to 16.7% and 12.5%, respectively. Overall, Hg and Cd were enriched rapidly in the soil.

The topsoil concentration of Cr(VI) at all sampling sites exceeded the background value in five years, due to the high annual net input. However, the average concentration was predicted to be stabilized at 9.72 mg/kg in the topsoil because of the low soil residual rate of Cr(VI) (Table 2).

Similar to Pb, the enrichment rate of As from 2011 to 2021 was predicted to be 168.9%, and the percentage of sampling sites where the topsoil Pb exceeds the background value would increase from 0 to 100%. The enrichment rates would decrease to 21.2% and 13.5%, respectively, from 2021 to 2026 and from 2026 to 2031.

3.4 Spatial trend simulation and risk warning

3.4.1 Spatiotemporal change of Pb concentration and early warning

The enrichment of Pb would first occur on both sides of the coal corridor and remain afterward (Fig. 3). In 2026, high level of Pb enrichment was predicted to occur on both sides of the trunk road, around the power unit chimney, and near the entrances and exits. In 2031, the enrichment would become more pronounced, concentrating on both sides of the corridor and near entrances and exits. Moreover, there would be a certain enrichment around the chimney of the power unit. Zhang et al. (2017) studied the Pb-enriched area in the coal chemical field and found that Pb emissions from coal chemical industry were relatively low, and the main enrichment area was located on both sides of the corridor and the highway. Zhu et al. studied the HM distribution on both sides of the highway suggesting that the use of leaded gasoline and Pb-containing lubricants contributed to a large amount of Pb emissions, which could indicate that the main sources of Pb in the soil were the lubricant application in the corridor and the automobile exhaust (Zhu et al. 2015).

Spatiotemporal distribution of Pb concentrations in the soil of the coal chemical factory (left to right for the year 2021, 2026 and 2031)

For the prevention and control of Pb, vehicle exhaust emissions along the traffic should be well controlled, lubricating oil should be used reasonably, and the friction of the transmission belt should be reduced. In addition, it is necessary to take measures to protect the staff working in the coal corridor and personnel commuting through the entrance and exit.

3.4.2 Spatiotemporal change of Hg concentration and early warning

The first enrichment area of Hg was observed in the dynamic unit and its south side, and the high Hg concentration in the soil would spread out around this area (Fig. 4), which indicated the source of Hg in the plant area is relatively monotonic. However, the enrichment rate of Hg is relatively fast. The concentration of Hg at the core region would increase steadily at around 0.050 mg/kg every five years, while the enrichment rate of the surrounding area would decrease with the distance from the core region.

Spatiotemporal distribution of Hg concentration in the soil of the coal chemical factory (left to right for the year 2021, 2026 and 2031)

A study suggested that more than 90% of the Hg was discharged with the flue gas, due to the high temperature in the coal-fired boiler (Sun and Bai 2008). Another study found that almost all Hg in the coal was released at a temperature of 815 °C (Zhang et al. 2003). Therefore, the soil Hg in the present study was predicted to mainly come from the gas emission of the power unit.

In terms of prevention and control of Hg in the plant area, the management of flue gas emissions in the plant area should be strengthened. Measures should be taken to protect staff members working under the flue gas discharge ports such as power units and logistics outlets.

3.4.3 Spatiotemporal change of Cd concentration and early warning

The core regions of the Cd enrichment were predicted to be located on both sides of the line from the northwest to the southeast of the plant and spread to farther distance (Fig. 5). These regions cover the main road in the factory, the downward direction of the flue gas discharge, and the coal corridor. In addition, the highest concentration of Cd would occur at the personnel entrance and exit and the coal corridor.

Spatiotemporal distribution of Cd concentrations in the soil of the coal chemical factory (left to right for the year 2021, 2026 and 2031)

The main soil pollutant on both sides of the roads in Beijing was found to be Cd, and the main sources were supposed to be the motor vehicle exhaust and dust generated by the abrasion and braking of motor vehicle tires (Wang 2014). A study suggested that a certain amount of Cd was emitted through flue gas during coal combustion (Xu et al. 2004). As a result, motor vehicles and coal-fired gas emissions were deducted to be the main Cd sources in the present study. By 2021, the highest concentrations of Cd would be concentrated around the personnel entrance and exit. By 2026, the first-line of Cd enrichment along the roadway-coal-coal corridor would be more noticeable, and the high concentration around the personnel entrance and exit will expand, and a core area will be formed along the coal corridor. By 2031, Cd high concentration area will cover both sides of the main road, the power unit, and the coal corridor. In general, Cd enrichment will be relatively fast and concentrated along the line of the personnel entrance and exit, the downwind of the flue gas vent-coal corridors.

To prevent and control Cd in the plant, the management of flue gas emission and the pollution control from traffic in the plant should be strengthened. Working along the coal corridors, near the personnel entrance and exit, and along the main roads, staff members should be protected from the hazard of Cd pollution.

3.4.4 Spatiotemporal change of Cr(VI) concentration and early warning

The highest soil Cr(VI) concentration in the plant was predicted to be located at the gasification unit and the power unit, and the concentration would decrease in a stepwise manner with the distance from the two units (Fig. 6). The estimated Cr(VI) concentration in the soil will exceed the Ningxia background value in the whole plant in 2021, and the maximum exceedance will reach 184.32%.

Spatiotemporal distribution of Cr(VI) concentrations in the soil of the coal chemical factory (left to right for the year 2021, 2026 and 2031)

A study found that the Cr concentration in the gasification slag dumps in the coal chemical industry was significantly higher than that in raw coal and flue gas (Zhang et al. 2017). Another study resulted in finding that the volatilization temperature of Cr was over 1500 °C, and the volatile amount was very low (Yang Li 2011). Therefore, the main source of Cr(VI) in the present study was proposed to be the slag dumps of the power unit and gasification unit.

To control Cr(VI) pollution in the plant, the slag dump management of gasification units and power unit should be strengthened to control the dust and make good seepage prevention. Measures should be taken to protect staff members working in gasification unit and power unit from Cr(VI) hazards.

3.4.5 Spatiotemporal change of As concentration and early warning

In 2021, the high As concentration in the soil will cover the eastern part of the gasification unit, power units, main roads in the west–east direction, and the coal corridor. The As concentration is relatively uniform between 24.590 and 27.762 mg/kg (Fig. 7). In 2026, there would be a droplet-like area of As around the flue gas discharge and its southeastern side, where the concentration would be significantly higher than the high concentration in 2021. In 2031, the droplet-like area would be further highlighted with a shift to the southeast, and the As concentration would be further increased.

Spatiotemporal distribution of As concentrations in the soil of the coal chemical factory (left to right for the year 2021, 2026 and 2031)

A study proposed that As was mainly converted into fly ash and gaseous oxides in coal combustion process (Lu et al. 2003). In contrast, the main source of soil As in this plant in the present study was predicted to be coal-combustion flue gas emissions. The management of flue gas emissions and the control of the As emission should be strengthened to prevent and control the As contamination. For staff working in the gasification unit and its downwind direction, necessary protective measures should be taken to prevent As pollution hazards.

1. (1)

   The soil environmental capacity model was used to predict the HM concentrations in the coal chemical industry with parameters obtained by the leaching experiment. With the HM baseline values in 2011, the soil concentrations of Pb, Hg, Cd, Cr and As in the year of 2021, 2026 and 2031 were predicted.

2. (2)

   The Kriging interpolation method was introduced with ArcGIS10.0 to achieve spatiotemporal simulations of the predicted HM concentrations in the soil. Spatiotemporal maps to predict the changing trends of HM concentrations were obtained. Prediction and early warning of HM pollution in the plant were also proposed.

3. (3)

   HMs can be rapidly enriched in the soil from discharge sources to the whole plant, due to the coal chemical industrial activities. The enrichment should be paid much attention both spatially and temporally.

4. (4)

   For the prevention and control of HMs in the plant, vehicle exhaust emissions along the traffic should be well controlled, lubricating oil should be used reasonably, and the friction of the transmission belt should be reduced. In addition, it is necessary to take measures to protect the staff working in the coal corridor and personnel commuting through the entrance and exit.

Also, the management of flue gas emissions in the plant area should be strengthened. Measures should be taken to protect staff members working under the flue gas discharge ports such as power units and logistics outlets.