Simulation of Sudden Water Pollution Accidents in Hunhe River Basin Upstream of Dahuofang Reservoir

: Dahuofang Reservoir is an important drinking water source for Shenyang, Fushun, Anshan, Liaoyang, and other cities. The water quality of the upstream inﬂow river directly affects the water supply safety. When a sudden water pollution accident occurs upstream of the reservoir, the pollution risk to the water source can be minimized if the variation rule of pollutant concentration along the course can be accurately simulated in time. Therefore, based on Mike 21, this paper established a hydrodynamic water quality coupling model of the Hun River basin upstream of Dahuofang Reservoir, and determined and veriﬁed the relevant parameters of the hydrodynamic model. In establishing the water quality model, the improved empirical frequency curve method was adopted to divide the high-ﬂow period, the level period, and the low-ﬂow period, so that the hydraulic conditions in each period were more reasonable. By a hypothetical scenario and working condition design, the suspended iron concentration and COD concentration along the course of a sudden water pollution accident were simulated. The diffusion rules of pollutants in different periods and under different working conditions were obtained. The most important objective was to obtain the six early warning index values in different hydrological periods, which allows the prediction of the scope and extent of the accident and provides a basis for ensuring the safety of the water supply at the water source.


Introduction
In recent years, frequent water pollution accidents have caused great threats to the safety of water sources everywhere, such as ship oil spills, nuclear leakage, pollutant discharge caused by industrial accidents, vehicles dumping dangerous goods, and so on. Because sudden water pollution accidents have uncertainty and other characteristics, once accidents occur, the relevant authorities are often caught off guard. For example, in 2005, an explosion at the Jilin Petrochemical Company's benzene plant caused about 100 tons of benzene to spill into the Songhua River, creating an 80 km-long pollution belt that cut off Harbin's water supply for five days; in February 2009, the water source of Yancheng City, Jiangsu Province was polluted by phenolic compounds, which led to the closure of two waterworks and affected the drinking water of hundreds of thousands of citizens; in 2020, a leak occurred in the tailing pond of the Luming mine in Yichun, Heilongjiang, where the concentration of molybdenum exceeded the standard by nearly 10 times, causing economic losses of 44,204,500 yuan. It can be seen that the occurrence of sudden water pollution accidents has a great impact on people's normal water safety [1]. Therefore, it is of great significance to simulate sudden water pollution accidents at water sources and analyze changing pattern of pollutant concentrations along their route, which is the basis for the relevant departments to develop emergency disposal measures [2,3]. 2 of 13 At present, for the simulation and analysis of sudden water pollution accidents, water quality models combined with simulation models and mathematical models are usually used to simulate the change in the concentration of pollutants in river water at home and abroad. Nowadays, there are many widely used water quality numerical model software, such as QUAL-2k, WASP, EFDC, Delft-3D, and MIKE [4][5][6][7][8][9][10]. Among them, Mike has become the most widely used software in water environment numerical simulation. Developed by the Danish Institute for Hydraulic Research, the Mike software has been widely used in the simulation of hydrodynamic, water quality, sediment erosion and other processes of major rivers, and has a good user interface and powerful front-andback treatment functions [11,12]. Yu et al. (2013) simulated the sudden water pollution accident of Daning Reservoir in Beijing based on Mike 21 FM. Taking the leakage accident of 5t nitrobenzene as an example, the study was carried out to obtain the diffusion speed, the proportion of excess area, and the emergency response time of pollutants under five different scenarios. Finally, two suggestions were put forward to reduce the probability of emergency. In this simulation, the wind direction was set as the only variable, and the two dominant wind directions in spring and summer, and autumn and winter were taken as the basis. However, in fact, the hydraulic conditions in different seasons of spring, summer, autumn, and winter were also different, which was not fully considered in the paper [13]. Wang et al. (2015) established a water quality simulation model for the Jing-Shi section of the South-North Water Transfer Central Line based on the hydrodynamic module and convective diffusion module of Mike 21, and made it achieve a better simulation effect through the rate and verification of the model. The model simulated the migration and diffusion of pollution masses, the time and concentration of pollution masses reaching the next water diversion, and the retreat gate after the over-tuning accident of dangerous goods transportation on the cross-drain highway bridge, guaranteeing the safety of the water supply by judging whether to open the retreat gate [14]. Zhai et al. (2019) used a two-dimensional water volume and quality model to simulate sudden water pollution accidents in the mainstream of the Yangtze River and the mainstream of the Jialing River in the Chongqing section of the Three Gorges Reservoir. The diffusion rate is the largest in the flood season, followed by the impoundment period, the declining period, and the dry period. This study provides a reference for choosing an appropriate water resource operation plan to reduce the peak concentration of river pollutants [15]. However, in the simulation of sudden water pollution accidents, the confirmation of hydraulic conditions was vague, and hydraulic conditions were important factors to ensure the accuracy of the model simulation [14]. Previous studies have focused on conventional pollution in the upstream of Dahuofang Reservoir, but the study of sudden water pollution accident is lacking.
Therefore, this study established a coupled model of hydrodynamic-water quality for the Hun River basin based on Mike 21 FM, and used an improved empirical frequency curve method to determine the hydraulic conditions and simulate the changes in pollutant concentrations when a sudden water pollution accident occurs, in order to provide a basis for the development of emergency plans for sudden water pollution accidents in the Hun River basin.

Study Area
Dahuofang Reservoir is the main drinking water source of Shenyang, Fushun, Anshan, Liaoyang, Yingkou, Panjin, and other cities in central Liaoning Province. The Hun River, Suzi River, and She River in the upstream of the reservoir are the main water sources of the reservoir, accounting for 52.7%, 37.1%, and 10.2%, respectively [16]. The geographical location of the reservoir and upstream basin is shown in Figure 1. The Hun River basin above the Dahuofang dam site is 169 km and controls a basin area of 2752 km 2 , which is the river with the largest amount of water entering the Dahuofang Reservoir, and the presence of chemical mines such as Qingyuan Limestone Mine and Fushun Hongtushan Mining Company (Fushun, China) Limited in the Hun River basin increases the possibility Water 2022, 14, 925 3 of 13 of sudden water pollution accidents. Therefore, it is of great significance to establish a hydrodynamic water quality model in Hunhe River basin and analyze the variation law of pollutant concentration when sudden water pollution accidents occur, so as to ensure the safety of the water supply in Dahuofang Reservoir. The depth of the Hunhe River in the research section is less than ten meters, which is surface water and can be studied with a two-dimensional mathematical model.

Construction of Mike 21 FM Model
In this study, a two-dimensional coupled model of hydrodynamic-water quality based on MIKE 21 FM was established to simulate the real-time change of pollutant concentration during the occurrence of water pollution incidents in the study area.

Fundamentals of Mike 21 FM Model
The average water depth of the study area was less than 10 m, and there was no obvious stratification phenomenon, so the two-dimensional model simulation can meet our needs. MIKE 21 FM uses unstructured mesh, and its numerical calculation method is the finite volume method, which has good conservation properties and can accurately deal with jet and discontinuous solutions. Therefore, this study used the MIKE 21 FM model to simulate.
The basic principle of the MIKE 21 FM two-dimensional hydrodynamic model is based on the Reynolds average stress equation of two-dimensional incompressible fluid, that is to say, the two-dimensional shallow water equation, which obeys the Boussinesq hypothesis and the hydrostatic pressure hypothesis [17].
The plane two-dimensional continuous equation of water flow is as follows: ∂h ∂t The momentum equation of the plane two-dimensional flow is as follows: In the formula, t is the time (s); x, y, and z are the Cartesian coordinate system; h is the total head (m); u is the velocity component in the direction of x (m/s); v is the velocity component in the direction of y (m/s); g is the acceleration of gravity (m/s 2 ); ρ is the density of water (g/m 3 ); ρ 0 is the relative density of water (g/m 3 ); p a is the atmospheric pressure (Pa); s xx , s xy , s yx , s yy are the components of radiation stress; S is the magnitude of the point source flow (m 3 /s); u s and v s are the flow velocities of the source-sink term (m/s).

Construction and Verification of Hydrodynamic Model
The construction of a hydrodynamic model is the basis of the coupled model of hydrodynamic water quality. First, we used the Mesh Generator tool in Mike Zero to generate the topography file, and then determined the calculation parameters and boundary conditions of the hydrodynamic model (HD) in Mike 21 FM. Finally, we analyzed the absolute and relative errors between the measured and simulated water levels at the Cangshi hydrological station to determine whether the model accuracy meets the requirements.

Grid Generation
The calculation range of this model was the Hunhe River from the upper reaches of Qingyuan limestone mine to Beizamu Bridge, with a length of about 38.8 km. In order to improve the simulation accuracy, the grid was encrypted within 1.5 km of the sewage outlet, as shown in Figure 2. Through the satellite maps and cross-section survey maps, the river boundary line and underwater terrain elevation can be obtained.
In the formula, t is the time (s); x , y , and z are the Cartesian coordinate system; h is the total head (m); u is the velocity component in the direction of x (m/s) ; v is the velocity component in the direction of y (m/s); g is the acceleration of gravity (m/s 2 ); ρ is the density of water (g/m 3 ); 0 ρ is the relative density of water (g/m 3 ); a p is the atmospheric pressure (Pa); xx s , xy s , yx s , yy s are the components of radiation stress; S is the magnitude of the point source flow (m 3 /s); s u and s v are the flow velocities of the source-sink term (m/s)

Construction and Verification of Hydrodynamic Model
The construction of a hydrodynamic model is the basis of the coupled model of hydrodynamic water quality. First, we used the Mesh Generator tool in Mike Zero to generate the topography file, and then determined the calculation parameters and boundary conditions of the hydrodynamic model (HD) in Mike 21 FM. Finally, we analyzed the absolute and relative errors between the measured and simulated water levels at the Cangshi hydrological station to determine whether the model accuracy meets the requirements.

Grid Generation
The calculation range of this model was the Hunhe River from the upper reaches of Qingyuan limestone mine to Beizamu Bridge, with a length of about 38.8 km. In order to improve the simulation accuracy, the grid was encrypted within 1.5 km of the sewage outlet, as shown in Figure 2. Through the satellite maps and cross-section survey maps, the river boundary line and underwater terrain elevation can be obtained.

Boundary and Initial Conditions
According to the simulation range of this model, the upstream boundary of the model was determined to be 1600 m upstream of Qingyuan Limestone Mine, and the daily average flow from 1 January 2014 to 31 December 2018 was used as the flow boundary, while the downstream boundary was Beizamu Bridge, and the average flow was used as the water level boundary.
The initial conditions of the hydrodynamic model include initial water level and initial velocity, which can be set to ensure the smooth start of the model. For this model, the initial water level of the model was set to the lowest water level value of 137.62 m at the downstream boundary, and the initial flow velocity was zero. (3) Eddy viscosity coefficient: the Smagorinsky formula was selected for control, and the constant was 0.3. (4) Rainfall and evaporation: the daily average rainfall and evaporation during the simulation time were used.

Simulation Results
The hydrodynamic model was validated based on the measured water level data from Beikouqian hydrological station (2014) and Cangshi hydrological station (2015-2018) (Before 2014, it was Beikouqian hydrological station, and later moved to Cangshi hydrological station). The simulated water level and measured water level of Beikouqian hydrological station (2014) are shown in Figure 3. The average absolute error was 0.144 m and the average relative error was 0.085%, of which the maximum absolute error was 0.192 m and the maximum relative error was 0.113%.The simulated water level of Cangshi hydrological station (2015-2018) and the actual measured water level were shown in Figure 4, with the average absolute error of 0.096 m and the average relative error of 0.063%, of which the maximum absolute error was 0.198 m and the maximum relative error was 0.129%. It can be seen that the calibration result was good and can be simulated.

Calculation Parameters
(1) Simulation time: the specific simulation time of the hydrodynamic model was from 1 January 2014 to 31 December 2018, with a simulation time step t Δ = 86,400 s and a total of 1825 steps. (2) River bed roughness: controlled by Manning coefficient, which was verified to be 41 m (1/3)/ s. (3) Eddy viscosity coefficient: the Smagorinsky formula was selected for control, and the constant was 0.3. (4) Rainfall and evaporation: the daily average rainfall and evaporation during the simulation time were used.

Simulation Results
The hydrodynamic model was validated based on the measured water level data from Beikouqian hydrological station (2014) and Cangshi hydrological station (2015-2018) (Before 2014, it was Beikouqian hydrological station, and later moved to Cangshi hydrological station). The simulated water level and measured water level of Beikouqian hydrological station (2014) are shown in Figure 3. The average absolute error was 0.144 m and the average relative error was 0.085%, of which the maximum absolute error was 0.192 m and the maximum relative error was 0.113%.The simulated water level of Cangshi hydrological station (2015-2018) and the actual measured water level were shown in Figure 4, with the average absolute error of 0.096 m and the average relative error of 0.063%, of which the maximum absolute error was 0.198 m and the maximum relative error was 0.129%. It can be seen that the calibration result was good and can be simulated.

Construction of Water Quality Model
A water quality model was established based on the parameters of hydrodynamic model calibration, and the ECO Lab Heavy Metal module and Transport module were

Construction of Water Quality Model
A water quality model was established based on the parameters of hydrodynamic model calibration, and the ECO Lab Heavy Metal module and Transport module were used to simulate the sudden heavy metal iron and chemical oxygen demand pollution accidents, respectively [18,19]. In addition, it was necessary to select reasonable simulation times for the water quality model, and to determine the initial conditions of pollution source through scenario setting and operating condition setting. According to the simulation, exceeding the multiple of maximum concentration, length of exceeding range, area exceeding the standard, time of pollutants reaching downstream boundary, time of pollutant concentration reaching its peak at the downstream boundary, and the exceedance of the multiple of pollutant concentration reaching its peak at the downstream boundary were selected as the pre-warning indexes that can be obtained under different working conditions.

Simulation Time
Due to the occurrence of sudden water pollution accidents, early warning is necessary to take relevant measures in time to ensure the safety of the reservoir water supply. Therefore, a reasonable choice of water quality model simulation time, that is to say, a reasonable determination of hydraulic conditions, is of great significance to the accurate simulation of sudden water pollution accidents.
In order to select reasonable hydraulic conditions, the improved empirical frequency division method was adopted by drawing the empirical frequency curve of annual runoff of Beikouqian hydrological station (changed to Cangshi hydrological station after 2014) from 1976 to 2018. The years with p = 10%, p = 50%, and p = 90% were selected to represent the abundant water year, the flat water year, and the dry water year, respectively. Then, according to the selected empirical frequency curves of daily average flow in the years of abundant water, flat water, and dry water, flows of p = 10%, p = 50%, and p = 90% were selected as the flows of the abundant water period, the flat water period, and the dry water period, respectively. Finally, the water quality model was established under the hydraulic conditions of the three periods, respectively.
The year of annual runoff p = 10% from 1976 to 2018 was calculated as 2013, the year of p = 50% was 2017, the year of p = 91% was 1978, and the year of p = 86.4% was 2018. To ensure the accuracy of the terrain, 2018 was chosen as the year of dry water. The empirical frequency p = 10% of the daily average flow in the year of abundant water 2013 was July 28, the empirical frequency p = 50% of the daily average flow in the year of flat water 2017 was September 21,and the empirical frequency p = 90% of the daily average flow in the year of dry water 2018 was May 4. Therefore, 28 July 2013, 21 September 2017, and 4 May 2018 were adopted as the hydraulic conditions for the abundant water period, flat water period, and dry water period.

Simulation of Sudden Pollution Accident
A Qingyuan limestone mine is located in Douhutun Village, Fushun. Iron is one of the elements with high content in the wastewater discharged from five nearby mines. If a sudden pollution accident occurs, a large amount of wastewater will be discharged into the river, which will have a huge impact on the water quality of Hunhe River, and even affect the normal water supply of Dahuofang Reservoir. Therefore, the accident simulation design is the diffusion and attenuation of the pollution mass after iron-containing wastewater is discharged beyond the standard. Because most heavy metal iron in wastewater exists in a suspended state, the heavy metal iron was simulated in this simulation. The location of the pollution source is shown in Figure 5. the river, which will have a huge impact on the water quality of Hunhe River, and even affect the normal water supply of Dahuofang Reservoir. Therefore, the accident simulation design is the diffusion and attenuation of the pollution mass after ironcontaining wastewater is discharged beyond the standard. Because most heavy metal iron in wastewater exists in a suspended state, the heavy metal iron was simulated in this simulation. The location of the pollution source is shown in Figure 5. Based on the simulation results of suspended iron concentrations during the abundant water, flat water, and dry water periods, the selected early warning index values were derived as shown in Table 1. In this study, suspended iron is simulated, which is not applicable to the heavy metal iron in the Environmental Standards for Surface Water GB3838-2002 [20]. Therefore, the suspended concentration requirements specified in the discharge standard for pollutants from Urban Wastewater Treatment Plants GB 18918-2002 [20] were used to implement the primary standard of 20 mg/L. Based on the simulation results of suspended iron concentrations during the abundant water, flat water, and dry water periods, the selected early warning index values were derived as shown in Table 1. In this study, suspended iron is simulated, which is not applicable to the heavy metal iron in the Environmental Standards for Surface Water GB3838-2002 [20]. Therefore, the suspended concentration requirements specified in the discharge standard for pollutants from Urban Wastewater Treatment Plants GB 18918-2002 [20] were used to implement the primary standard of 20 mg/L.

Simulation of Sudden Chemical Oxygen Demand (COD) Pollution Accident
Hongtoushan Mining Co., Ltd. (Fushun, China) is responsible for the sewage discharge of three companies, which is discharged into the drinking water source area in front of Hunhe Beikou by means of dark pipe discharge, with an average annual discharge of about 950,000 tons. The concentration of COD in wastewater was selected as the target pollutant to simulate the dispersion and attenuation of pollutant mass. The location of the source of pollution is shown in Figure 6.
Water 2022, 14, x FOR PEER REVIEW 9 of 12

Simulation of Sudden Chemical Oxygen Demand (COD) Pollution Accident
Hongtoushan Mining Co., Ltd.(Fushun,china) is responsible for the sewage discharge of three companies, which is discharged into the drinking water source area in front of Hunhe Beikou by means of dark pipe discharge, with an average annual discharge of about 950,000 tons. The concentration of COD in wastewater was selected as the target pollutant to simulate the dispersion and attenuation of pollutant mass. The location of the source of pollution is shown in Figure 6. According to the simulation results of COD concentration during the abundant water period, normal water period, and low water period, the selected warning index values are shown in Table 2. Based on the Environmental Standards for Surface Water GB3838-2002, the COD concentration is 20 mg/L for Class III water.  According to the simulation results of COD concentration during the abundant water period, normal water period, and low water period, the selected warning index values are shown in Table 2. Based on the Environmental Standards for Surface Water GB3838-2002, the COD concentration is 20 mg/L for Class III water.

Results and Discussion
According to Tables 1 and 2, by analyzing the values of warning indicators for two pollutants at different periods, different flows, and different concentrations, the following results were obtained.
(1) The higher the discharge concentration is, the greater the exceeding multiple of maximum concentration and the time of pollutant concentration reaching its peak at the downstream boundary, and the more serious the pollution to the reservoir, but there is little relationship with the discharge. For example, when the sudden COD pollution accident occurs, the maximum concentration exceeding multiple of the simulated discharge flow 50 m 3 /s, discharge concentration 50 mg/L, discharge flow 50 m 3 /s and discharge concentration 100 mg/L in the abundant water period were 1.26 and 3.596, respectively, while the maximum concentration exceeding multiples of the simulated discharge flow 50 m 3 /s, discharge concentration 50 mg/L, discharge flow 100 m 3 /s and discharge concentration 50 mg/L in abundant water period was 1.26 and 1.373 respectively. This may be because the dilution effect of the river is more effective for low-concentration pollution, and the dilution capacity for high-concentration pollution is limited/ The higher the initial concentration of pollutants, the more serious the pollution of the reservoir when the river water enters the reservoir, which is consistent with the conclusion reached in the numerical simulation of sudden environmental events in the Fen River reservoir by Yang et al. (2017) [21]. Therefore, when an accident occurs, the primary task is to grasp the concentration of pollutants in time, which is conducive to understanding the severity of the accident.
(2) The length and area of the exceeding standard range of suspended iron concentration were 36.88 (±0.262) km and 3.475 (±0.018) km 2 , respectively. The length and area of excess COD concentration were 9.29 (±0.16) km and 1.06 (±0.07) km 2 , respectively. On the whole, the length and area of the excess concentration range are only related to the location of the pollution source under the assumption that the pollution mass reaching the downstream boundary concentration exceeds the limit; they have little to do with the time of pollutant discharge, discharge volume, and discharge concentration. Specifically, the length of the exceeding standard range of pollution source in the same position is from large to small: dry water period > flat water period > abundant water period. For example, when a sudden heavy metal iron pollution accident occurs, with a simulated discharge flow rate of 50 m 3 /s and discharge concentration of 50 mg/L, the length of the exceedance range in the abundant water period, the flat water period, and the dry water period are 36.633 km, 36.687 km, and 36.923 km, respectively. This conclusion is consistent with the conclusion reached by Bai et al. (2013) in the early warning of sudden water pollution incidents in the Yellow River [22]. In the future, when dealing with accidents, the location close to the pollution source will be more affected by pollution, and more resources should be invested.
(3) With the same pollution source location and the same time period, the time of pollutants reaching the downstream boundary decreases with the increase of discharge flow and concentration. However, the time of pollutant concentration reaching its peak at the downstream boundary is only related to the discharge flow, decreasing with the increase of the discharge flow. With the same pollution source location and different time conditions, the time of pollutants reaching the downstream boundary and the time of pollutant concentration reaching its peak at the downstream boundary are respectively, from large to small: dry water period > flat water period > abundant water period. For example, when a sudden COD pollution accident occurs, the time of pollutants reaching the downstream boundary is three and a half hours, four hours, and five hours, respectively, when the simulated discharge is 50 m 3 /s and the discharge concentration is 50 mg/L in abundant water period, the flat water period, and the dry water period. The pollutant diffusion rate is related to the period of the pollution accident, and the pollutant diffusion rate from fast to slow is: abundant water period > flat water period > dry water period,. Therefore, the time of pollutant reaching the downstream boundary and time of pollutant concentration reaching its peak at the downstream boundary, that is to say, the time of early warning response, from long to short is: dry water period > flat water period > abundant water period. This conclusion is consistent with that obtained in the simulation of sudden pollution incidents in Nanchang section of Ganjiang River by Shu et al. (2019) [23]. However, in the simulation of COD in this paper, the pollutants reach the downstream boundary at the same time under the same conditions (conditions: discharge flow rate: 100 m 3 /s, discharge concentration: 50 mg/L; discharge flow rate 100 m 3 /s, discharge concentration 100 mg/L) during the flat water period and the dry water period, which is mainly due to the fact that the location of the pollution source is closer to the downstream boundary, at only 9.04 km. The pollutants reach the downstream boundary at the same time because the discharge flow rate is larger and the pollutant diffusion rate is basically the same. Therefore, if the accident occurs during the abundant water period, the degree of danger is greater, and the management department should strengthen water quality management during this period.

Conclusions
By establishing a coupled hydrodynamic water quality model based on Mike 21 software for the upper reaches of Dahuofang Reservoir in the Hun River basin, the numerical simulation of sudden suspended iron pollution and COD pollution under different working conditions was carried out. According to the simulation results, the index values of the selected six early warning indexes were obtained and the pollutant diffusion law was found, which provides a basis for the safe water supply of Dahuofang Reservoir.
The study shows that pollutant concentration gradually decreases with time during the diffusion process. The pollutant emission concentration can affect the exceeding multiple of maximum concentration and the exceeding multiple of pollutant concentration reaching their peaks at the downstream boundary. The higher the discharge concentration is, the larger the exceedance multiplier is, though this has little relationship with the discharge flow rate. Under the premise that the default pollutant mass reaches the downstream boundary with excessive concentration, the length of the exceedance range and the area of the exceedance range are only related to the location of the pollutant source, and have little relationship with the pollutant discharge period, discharge flow rate, or discharge concentration. The time of the pollutant reaching the downstream boundary and the time of the pollutant concentration reaching its peak at the downstream boundary from large to small are as follows: dry water period > flat water period > abundant water period.
In an actual watershed accident, there may be more than one pollutant, but the two pollutants in this paper need to be calculated by two independent modules: Transport and Eco Lab. Therefore, in the future, we can focus on establishing a more complete mathematical model that can calculate the simultaneous pollution of multiple pollutants, in order to manage and ensure water quality.