Introduction

On August 12, 2015, a fire erupted in the arrival area of the hazardous goods warehouse of Ruihai Company in Binhai New District of Tianjin, China [1]. A subsequent investigation discovered approximately 129 types of chemicals, including 1300 metric tons of ammonium nitrate and potassium nitrate, 500 metric tons of sodium metal and magnesium metal, 43.4 tons of potassium dichromate, 315.8 tons of butanone, and hundreds of tons of organic compounds stored in the warehouse. The incident led to the pollution in the local area and threatened human health after some of the chemicals leaked into the soil–water system.

Hazards, such as nitrate, heavy metallic salt, and various toxic organic compounds [2], which were likely to bring about significant casualties and loss of property once they leaked into the soil and groundwater, highlighted the severity of this explosion accident. Potassium dichromate is a toxic and carcinogenic chemical that is classified as primary carcinogen by the International Agency for Research on Cancer. It displays stronger cytotoxicity and mouse acute toxicity than does trivalent chromium [3]. Butanone is irritating to the eyes, nose, larynx, and mucous membrane, and it leads to the death of aquatic organisms [4].

Efforts have been made to better understand the hazard transport in various unsaturated or saturated conditions. Li et al. [5] researched methyl tert-butyl ether (MTBE) dissolution in saturated porous media to obtain MTBE dissolution data with different groundwater velocities, initial MTBE saturation, and grain size of porous medium. Dev et al. [6] discussed the potentiality of using anaerobic packed-bed bioreactor for the treatment of acid mine drainage. Multiple process parameters, such as pH, hydraulic retention time, concentration of marine waste extract, total organic carbon, and sulfate, were optimized together using Taguchi design. The study also established the future scope for developing an efficient treatment process for sulfate and metal-rich mine wastewater on a large scale.

Kim et al. proposed a new approach for integrating the quasi-distributed watershed model called SWAT with the fully distributed groundwater model MODFLOW [7]. The verified model demonstrated that an integrated SWAT–MODFLOW is capable of simulating the spatial–temporal distribution of groundwater recharge rates, aquifer evapotranspiration, and groundwater levels. Kuznetsov et al. presented a new approach for a numerical scheme for 3D variably saturated flow using the quasi 3D Richards equation and finite-difference scheme [8]. Computationally, the quasi-3D code was more efficient by an order of 10–300%, while being accurate with respect to the benchmark fully 3D variable saturation code, when the capillary fringe was considered.

Some accident investigations have been conducted to identify human and organizational factors. Zhou et al. developed a modified version of the Human Factors Analysis and Classification System (HFACS), namely HFACS-Hazardous Chemicals, to identify the human factors involved in Chinese HCAs [9]. Zhou et al. selected the Tianjin 8·12 accident as an investigation case to present a simple and effective compound model in view of the advantages of different cause analysis models [10]. Most previous studies focused on the facts and lessons related to the accident, but not on the transport process of specific chemicals entering the soil–groundwater system.

This work aims to investigate the migration of potassium dichromate and butanone in the coastal soil–groundwater system of Binhai through a coupling unsaturated–saturated numerical model by incorporating the HYDRUS model into the MODFLOW/MT3D model. The soil properties were tested, and the area hydraulic conditions were explored. The verified numerical models were used to explore the substances’ migration distance and boundaries during a period of time. On the basis of the results, a migration boundary was proposed for delimiting contaminated zones. The obtained migration results are expected to be useful in instructing environment remediation in similar accidents.

Materials and Methods

Research Area

The explosion site (39.2N 117.4E) is located in Tanggu of Binhai, China (Fig. 1). It is situated 2.7 km away from the coastline in the east and about 0.3 km away from Binhai Road in the west. As a result of floodplain, the upper soil is dredger filled, consisting of fine-grained sediment [11, 12]. The elevation of the ground surface is 2–3 m above the mean sea level, and the aquifer table is located about 5 m below the soil surface. Binhai has a temperate monsoon climate, an annual average temperature of 7.5–12.3 °C, and an average annual precipitation of about 645 mm in the coastal region.

Fig. 1
figure 1

Location of the explosion site

Full size image

Properties of Soil

The soil profile information obtained from geology borehole-G2 is listed in Table 1.

Table 1 Geology borehole-G2 description
Full size table

To determine the soil properties, soil samples were collected from a few points near the explosion site. The soil was naturally air-dried, crushed, and sieved through the American Society for Testing and Materials standard sieves. Soil particles, which are finer than 2 mm, were applied in the test experiments. A number of soil sample properties, such as mechanical composition, pH, total salt content, organic matter, and hydraulic conductivity, were evaluated in this study. The physical and chemical properties are presented in Table 2.

Table 2 Physical and chemical properties of soil
Full size table

Hazardous Chemicals

The physical and chemical properties of potassium dichromate [13] and butanone [14] are listed in Table 3.

Table 3 Properties of potassium dichromate and butanone
Full size table

Numerical Simulation

After leaking into the soil from the surface, the transport process of hazards consisted of infiltration through the vadose zone and migration in the groundwater after reaching the water table. Given the great difference between the two zones, establishing an integrated numerical model is of great significance to investigate water and solute transport in the whole soil–water system [15]. One-dimensional or two-dimensional HYDRUS models have been proven to be efficient methods to describe the process in an unsaturated zone, considering the thickness and significant hydraulic gradient of such a zone [16]. In the saturated zone, MODFLOW and MT3D packages are applied to describe water and solute transport in a three-dimensional model, respectively. Being fully incorporated into the MODFLOW program, the HYDRUS package provides MODFLOW with recharge fluxes at the water table, while MODFLOW provides HYDRUS with the position of the groundwater table that is used as the bottom boundary condition in the package [17]. The schematic of the coupling model is presented in Fig. 2.

Fig. 2
figure 2

Schematic of the coupling model

Full size image

Mathematical Model

HYDRUS-2D numerically solves the Richards equation for saturated–unsaturated water flow and Fickian-based advection–dispersion equations for heat and solute transport [18].

MODFLOW and MT3D are modules integrated into the Groundwater Modeling System, which is a comprehensive graphical user environment for performing groundwater simulations. MT3D is a modular three-dimensional transport model for the simulation of advection, dispersion, and chemical reactions of dissolved constituents and is used in conjunction with MODFLOW in flow and transport simulation [19].

Water Flow Equation

The relation between water flux and soil pressure head in both unsaturated and saturated zones is governed by the Richards equation, which has the following form [20]:

$$q = - k(\theta )\left( {\frac{\partial h}{\partial x} + \frac{\partial h}{\partial y} + \frac{\partial h}{\partial z} + 1} \right),$$
(1)

where q is the water flux, L/T; h is the soil pressure head, L; k(θ) is the hydraulic conductivity as a function of the soil water content, L/T; and θ is the volumetric water content, L3/L3.

Van Genuchten model [21] is widely used to describe the relationship between soil pressure head h and water content θ, which is defined by the following equation:

$$\theta = \theta_{\text{r}} + \frac{{\theta_{\text{s}} - \theta_{\text{r}} }}{{[1 + (\alpha h_{\text{s}} )^{n} ]^{m} }},$$
(2)

where θr is the residual water content, L3/L3; θs is the saturated water content, L3/L3; hs is the pressure head (hs = |h| for h < 0 and hs = 0 for h ≥ 0), L; and α, n, and m are parameters (m = 1 − 1/n).

The water flow is defined by the following equation [22]:

$$\frac{\partial }{\partial x}\left( {K_{x} \frac{\partial h}{\partial x}} \right) + \frac{\partial }{\partial y}\left( {K_{y} \frac{\partial h}{\partial y}} \right) + \frac{\partial }{\partial z}\left( {K_{z} \frac{\partial h}{\partial z}} \right) + q_{\text{s}} = S_{\text{s}} \frac{\partial h}{\partial t},$$
(3)

where Kx, Ky, and Kz are the hydraulic conductivity in the direction of the x, y, and z coordinates, respectively, L/T; qs is the sinks/sources, 1/T; and Ss is the specific storage of the porous material, 1/L.

Solute Transport Equation

The solute transport equation is defined by the following advection–dispersion equation [22]:

$$\frac{{\partial \left( {\theta RC} \right)}}{\partial t} = \frac{\partial }{{\partial x_{i} }}\left( {\theta D_{ij} \frac{\partial C}{{\partial x_{j} }}} \right) - \frac{\partial }{{\partial x_{i} }}(q_{i} C) + q_{\text{s}} C_{\text{s}} - \lambda_{1} \theta C - \lambda_{2} \rho_{\text{b}} \overline{C} ,$$
(4)

where C is the liquid-phase concentration, M/L3; \(\overline{C}\) is the sorbed concentration, M/L3; qs is the source or sink term in the water flow equation, 1/T; Cs is the concentration of source or sink term, M/L3; λ1 and λ2 are the liquid-phase reaction rate coefficient and adsorption phase reaction rate coefficient, 1/T, respectively; D is the dispersion coefficient, L2/T; ρb is the bulk density; and M/L3; and R is the retardation factor.

Results and Discussion

Unsaturated Simulation

Numerical Model

The model area is presented in Fig. 3. The lateral flow is very weak in the unsaturated zone. Thus, the 2D model in the x and z directions is capable of simulating hazard migration without considering the process in the y direction [23]. The two-dimensional vertical section is 3 m in depth [24] from the soil surface to the water table and about 500 m in width from west Haibin Road to east Yuejin Road.

Fig. 3
figure 3

Position of two-dimensional profile at the accident site

Full size image

The two-dimensional finite-element contaminant transport model is presented in Fig. 4. The solute transport boundary conditions are as follows: EF is the pollutant discharge boundary condition, and AB, BC, and A′C are the no-flux boundary conditions. The water flow boundary conditions are as follows: AA′ is the atmospheric boundary condition, and BC is the constant head boundary. At the initial time, the pressure head above the water table obtained through linear interpolation was negative [25,26,27], and the initial contamination concentrations were zero. The rainfall data, as obtained from the Tianjin Water Affairs Bureau, are presented in Fig. 5. The model parameters are listed in Table 4.

Fig. 4
figure 4

Contaminant transport model in the unsaturated zone

Full size image
Fig. 5
figure 5

Rainfall data of the latest 2 years

Full size image
Table 4 Parameters of water flow and solute transport
Full size table

The model was verified through a field experiment conducted at a site located about 800 m near the explosion site because the actual explosion site was not available. The hydrogeological properties of the experiment site were similar to those of the explosion site. The experiment site is shown in Fig. 6.

Fig. 6
figure 6

Field experiment site

Full size image

The prepared potassium chloride solution was leaked on the soil surface for a duration of 6 h. After the solution leakage, soil solution samples were taken every 24 h at depths of 1.0, 2.0, and 2.5 m under the leakage surface until the concentration reached a steady state. The K+ concentration of the samples was tested.

A comparison between the experimental and simulation results at the monitoring points is presented in Fig. 7. Judging from the K+ breakthrough curve (BTC), the variation trend of the simulation is approximately consistent with that of the experiment. However, the errors were up to about 30% at some stages, and the experiment transport rate was higher, which may be due to the strong driving force in filling the solution and the preferential flow [28,29,30] in the soil. According to the comparison results, the model was assumed to be reliable in simulation research.

Fig. 7
figure 7

BTCs of experiments and simulation results at monitoring points

Full size image

Analysis on Potassium Dichromate

The simulation result of Cr2O 2-7 transport over 730 days is presented in Fig. 8, and the BTCs at the monitoring points are presented in Fig. 9. Judging from the contamination plume, Cr2O 2-7 reached a depth of 1.2 m at the beginning of 2 months because of the heavy rain in this period. After 6 months, the migration depth was up to 1.7 m, and the concentration in liquid phase began to decrease because of soil adsorption. In the horizontal direction, the migration process was very weak compared with that in the vertical direction as a result of weak lateral flow [31] in the vadose zone.

Fig. 8
figure 8

Contamination plume of Cr2O 2-7 in the unsaturated zone

Full size image
Fig. 9
figure 9

BTCs of Cr2O 2-7 at the monitoring points

Full size image

In Fig. 9, BTCs showed that the variation trends of concentration vary from the depths and time stages. At depths of 0.1, 0.3, and 0.5 m, the concentration reached peak values at 80, 130, and 210 days, respectively. The solute migration rate was in line with the trend of precipitation intensity. Cr2O 2-7 reached the water table at a depth of 3 m after about 330 days. According to the mass balance information, after 2 years, about 12% of the total Cr2O 2-7 entered the aquifer, while the rest remained in the soil at a depth of 3 m. The Cr2O 2-7 entering the aquifer transported in the groundwater system [32], which is different from the unsaturated zone.

Analysis on Butanone

The simulation results of butanone migration are presented in Fig. 10. The transport rate of butanone was much slower than that of Cr2O 2-7 . Two months after the leakage, butanone reached a depth of only about 0.35 m and then reached a depth of about 0.67 m after 6 months. Obviously, the largest difference in the transport processes of butanone and Cr2O 2-7 was that in 2 years, butanone reached a depth of 1.5 m instead of the water table. To predict the migration distance of butanone in the long term, a simulation of a 5-year period was conducted. Compared with the distance in the first 2 years, butanone transported about 0.7 m farther in another 3 years, but still did not reach the aquifer.

Fig. 10
figure 10

Contamination plume of butanone in the unsaturated zone

Full size image

According to the BTCs in Fig. 11, in the last 2 years, the concentration at the monitoring points changed smoothly without significant fluctuations and showed low sensitivity to precipitation. At depths of 0.1 and 0.3 m, the concentrations reached peak values at about 120 and 480 days, respectively. Under the same conditions, the soil adsorption to butanone was much stronger [33] than that to Cr2O 2-7 . In addition, the advection process in water was much weaker due to the low water solubility of the soil [34]. At a conservative estimate, a small amount of butanone would reach the water table due to the removal of the unsaturated zone.

Fig. 11
figure 11

BTCs of butanone at the monitoring points

Full size image

Saturated Simulation

Numerical Model

Considering the significant water flow in the saturated zone, a three-dimensional model was built to track the solute transport process after entering the aquifer. To simplify the simulation, the bottom boundary flux of the HYDRUS model was incorporated as the recharge flux of the MODFLOW model. The model boundary and the MODFLOW model are presented in Fig. 12.

Fig. 12
figure 12

MODFLOW model conducted on the research area

Full size image

The water flow boundary conditions are as follows: the east, north, and south are constant head boundary conditions; Nanfeng Road in the west is the specified flow boundary conditions; and the area surface is the recharge boundary conditions. The solute transport boundary conditions are as follows: all the boundaries are no-flux boundary conditions. The background concentrations of hazards are assumed to be zero in the research area. The flow simulation results computed by MODFLOW are read by MT3D and utilized as the flow field for the transport portion of the simulation [35].

Analysis on Potassium Dichromate

From the flow field presented in Fig. 12, we could determine that the water level in the west was higher than that in the east coast. The flow direction of the groundwater in the area is generally from the northwest to the southeast at a low flow rate.

The migration results of Cr2O 2-7 in the saturated zone are shown in Fig. 13. First, judging from the contamination plume, Cr2O 2-7 had a greater tendency to migrate in the east and south directions than in the west and north directions. Two years after the leakage, the farthest migration horizontal distance of Cr2O 2-7 was about 73 m in south. To predict the potential migration distance of Cr2O 2-7 in the long term, a simulation of 3-, 4-, and 5-year periods was conducted. The contaminant reached a horizontal distance of about 105, 127, and 161 m in groundwater, respectively. According to the concentration distribution, the contaminant tended to accumulate at the section under the leakage source area and reached a farther distance with a low concentration.

Fig. 13
figure 13

Contamination plume of Cr2O 2-7 in the saturated zone

Full size image

On the basis of the above results, the unsaturated soil zone played a vital role in the removal of contaminants that infiltrated from the soil surface. Organic matter showed different transport properties from that of inorganic matter in the soil–water system. Butanone transported much more slowly than Cr2O 2-7 in unsaturated soil without reaching the aquifer in 5 years. In comparison, a portion of Cr2O 2-7 reached the aquifer and transported with the groundwater flow to a much farther distance, especially in the horizontal direction. Considering the harm posed by Cr2O 2-7 to the environment, measures should be taken immediately to clean up the contamination in the potential covered areas and the migration boundary [36, 37]. Different renovation strategies for various contaminants should be applied because of their different migration rates and contamination plumes [38].

Conclusions

An estimation of the coupling two-dimensional HYDRUS model and three-dimensional MODFLOW/MT3D model to address contamination mobility in the coastal soil–groundwater system of Binhai showed that Cr2O 2-7 and butanone displayed different migration properties under the same conditions. Cr2O 2-7 migrated relatively faster in the unsaturated zone than did butanone and reached the water table in about 1 year. In comparison, butanone reached a depth of only 1.5 m after 2 years and was unlikely to enter the aquifer even 5 years later with a migration distance of about 2.2 m, according to the simulation prediction. Driven by underground water, the Cr2O 2-7 that entered the aquifer migrated about 161 m toward the southeast. The contamination plume covered mainly the area in the southeast direction because of the groundwater flow direction. Measures should be taken immediately after leakage to clean up the contamination in the upper soil in case of arrival at the water table.