Migration of Cr2O72- and Butanone in Soil and Groundwater System After the Tianjin Port 8·12 Explosion

After the Tianjin Port 8·12 explosion, an enormous amount of potassium dichromate (K2Cr2O7) and butanone (CH3COCH2CH3) leaked into the coastal soil–groundwater system, which potentially threatened the environment and human health. Determining the transport process of hazards is necessary to establish guidelines for remediating the contaminated area. This work aims to investigate the migration of potassium dichromate and butanone in the coastal soil‒groundwater system through a coupling unsaturated–saturated numerical model, incorporating the HYDRUS model into the MODFLOW/MT3D model. In the unsaturated zone, two-dimensional HYDRUS model was applied, and its recharge flux at the bottom boundary was utilized as the input of MODFLOW/MT3D model in the saturated zone. Results showed that Cr2O72- migrated much faster than butanone in the unsaturated zone and reached the water table in about 1 year. In comparison, butanone was unlikely to enter the aquifer even 5 years later with a migration depth of about 2.2 m. Driven by groundwater, the Cr2O72‒ that entered the aquifer migrated about 161 m toward southeast 5 years later. In the saturated zone, the contamination plume covered mainly the southeast area due to the groundwater flow direction.


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.

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.

Properties of Soil
The soil profile information obtained from geology borehole-G2 is listed in Table 1.
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.

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

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, MOD-FLOW 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.

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]: 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, L 3 /L 3 .
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: where θ r is the residual water content, L 3 /L 3 ; θ s is the saturated water content, L 3 /L 3 ; h s is the pressure head (h s = |h| for h < 0 and h s = 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]: where K x , K y , and K z are the hydraulic conductivity in the direction of the x, y, and z coordinates, respectively, L/T; q s is the sinks/sources, 1/T; and S s 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]: where C is the liquid-phase concentration, M/L 3 ; C is the sorbed concentration, M/L 3 ; q s is the source or sink term in the water flow equation, 1/T; C s is the concentration of source or sink term, M/L 3 ; λ 1 and λ 2 are the liquid-phase reaction rate coefficient and adsorption phase reaction rate coefficient, 1/T, respectively; D is the dispersion coefficient, L 2 /T; ρ b is the bulk density; and M/L 3 ; and R is the retardation factor.

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 Fig. 2 Schematic of the coupling model 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. 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.
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.
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.

Analysis on Potassium Dichromate
The simulation result of Cr 2 O 7 2transport 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, Cr 2 O 7 2reached 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.
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. Cr 2 O 7 2reached 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 Cr 2 O 7 2entered the aquifer, while the rest remained in the soil at a depth of 3 m. The Cr 2 O 7 2entering 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 Cr 2 O 7 2-. 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 Cr 2 O 7 2was 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.
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 Cr 2 O 7 2-. 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.

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

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 Cr 2 O 7 2in the saturated zone are shown in Fig. 13. First, judging from the contamination plume, Cr 2 O 7 2had 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 Cr 2 O 7 2was about 73 m in south. To predict the potential migration distance of Cr 2 O 7 2in 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.
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 Cr 2 O 7 2in unsaturated soil without reaching the aquifer in 5 years. In comparison, a portion of Cr 2 O 7 2reached the aquifer and transported with the groundwater flow to a much farther distance, especially in the horizontal direction. Considering the harm posed by Cr 2 O 7 2to 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 Cr 2 O 7 2and butanone displayed different migration properties under the same conditions. Cr 2 O 7 2migrated 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.