Numerical Simulation of Bohai Oil Spill in the Winter Sea Ice Period

The Bohai Sea is a seasonal icy sea area that has the lowest latitude of any sea experiencing icing in the northern hemisphere, and simulation studies on oil spills during its sea ice period are the key to analyzing winter oil spill accidents. This study applied the three-dimensional free surface to establish a high-resolution hydrodynamic model and simulate tidal distributions in the Bohai Sea. Then, the oil spill model of the open sea area and thermodynamic model were combined to establish a numerical model for the Bohai oil spill during the winter sea ice period. The hydrodynamic model and sea ice growth and melting model were verified, and the parameters were adjusted based on the measured values, which indicate that the numerical model established in this paper is of high accuracy, stability and ubiquity. Finally, after checking the calculations repeatedly, the diffusion coefficient for the Bohai Sea was determined to be 1.0×10–7 m2/s. It is better that the comprehensive weathering attenuation coefficient is lower than that of a non-winter oil spill, with 1.3×10–7 m2/s being the most appropriate coefficient. This study can provide the reliable technical support for the operational safety and reduction in losses caused by winter oil spill accidents for the petroleum industry.


Introduction
Icing is a typical natural phenomenon that occurs every winter in the northern seas of China. The Bohai Sea in China has the lowest latitude of any sea experiencing icing in the northern hemisphere. Every winter, icing occurs at different degrees, and in severe cases, sea ice disasters may occur, threatening oil platform and vessel safety. For example, in 1969 and 2010, the most serious sea ice disasters ever recorded occurred in the Bohai Sea. The ice disaster in 1969 occurred during a severe ice period when the edge of the drifting ice approached the Bohai Strait, and the entire Bohai Sea surface was almost covered by sea ice. It was the most serious sea ice disaster ever recorded in China. The 1969 Bohai Sea ice disaster brought sea ice research and forecasting to the attention of Chinese scholars and marks the starting point of China's sea ice research. Therefore, it is of great practical significance for the disaster prevention and mitigation to predict ice conditions in the Bohai Sea, study the law of sea ice growth and melting combined with met-eorological conditions, simulate sea ice evolution with numerical models, and establish and improve sea ice forecasting models. The ice disaster mentioned above affected the exploration and exploitation of offshore oil, blocked the navigation channel, froze over the ports, interrupted transportation, and even damaged the port facilities. It even crashed and destroyed vessels, offshore oil platforms and other maritime structures, posing a great threat and causing destruction to shipping and maritime engineering operations. The above only describes the direct impacts of the sea ice disaster. An oil spill during the sea ice period will produce more serious consequences. For example, in 1989, the largest oil spill in American history occurred in the icy sea area near Alaska. Exxon's giant tanker Exxon Valdez failed to return to its normal channel in time when trying to avoid drift ice on the sea and was stranded in Alaskan waters, resulting in an overflow of 37600 tons of crude oil (a quarter of the ship's load). The oil spill then spread into an approximately 1300 km 2 oil spill area and gradually drifted and ex-panded over time due to hydrodynamic forces. The direct victims of the accident were coastal wildlife and fishery resources. Some experts predicted that the long-term environmental impacts of the oil spill during the ice period will last for more than a decade. After the accident, Exxon cleaned and recovered the oil spill at a cost of USD 2 billion in just one year, but the effect was not significant, with only 13% of the total oil spill being recovered. It can be seen that an oil spill in an icy sea area is different from that in an open sea area (Yu et al., 2016) in that the former will cause more harm to the sea area. Specifically, due to the existence of sea ice, oil spills diffuse at relatively slower rates but last longer. In addition, an oil spill in an iced area is very difficult to be cleaned and requires large amounts of manpower and material resources. In that way, it poses key technical problems for the vigorous development of the petroleum industry in the 21st century, which primarily regard how to ensure the safe operation of the oil industry during winter ice periods and to timely grasp the dynamics of oil film drift expansions resulting from oil spills in ice area. However, scientific and reliable oil spill prediction models for sea areas during winter sea ice periods are also the key to finding a solution to this problem.
The environmental conditions in sea areas with ice differ from those in the ice-free sea area. This is mainly because the temperatures of the sea water and the atmosphere in the sea area with ice are lower than those in the ice-free sea area, and the sea water temperatures are close to or equal to the freezing point of sea water. Additionally, due to the existence of the medium of sea ice, oil spills interact not only with the atmosphere and sea water but also with sea ice, which complicates various processes of the behaviours and fates of oil films. In recent years, numerous researchers in foreign countries have studied oil spill behaviours in icy sea areas by focusing on the interactions between oil and ice (Wilkinson et al., 2007;Stanovoy and Eremina, 2012). Canadian scientists are leading the observation and simulation of oil spill behaviours in icy sea areas. Ross and Dickins (1987) conducted three controlled oil spill tests in the sea near Breton Island on the eastern coast of Canada in March 1986. Metge and Telford (1980) studied the interactions of oil mixed with pancake ice and sticky ice through laboratory experiments. To investigate an oil spill in an ice area, Getman and Schultze (1976) studied an ice trough with an equilibrium thickness of #2 fuel oil and Prudhoe Bay crude oil in a stationary ice rink, with the research mostly adopting the experimental method. Domestic research on the behaviours and fates of oil spills in icy sea areas began relatively late. In addition to observational studies, oil spill accidents are mainly studied through laboratory or offshore field tests, and applications of numerical methods remain limited to numerical qualitative descriptions at present (Chen et al., 2011;Fu et al., 2013;Guo et al., 2013;Zhao et al., 2014;Diao et al., 2016). Few have carried out oil spill simulation and prediction research for the Bohai Sea during the winter ice period based on hydrodynamic models and sea ice and oil spill models. To undertake timely and effective emergency measures for oil spill accidents during the ice period, provide accurate and effective quantitative information for various government departments, and study the drift ranges and expansion directions of oil film pollution after oil bodies leak into the sea during the winter sea ice period or even the ice disaster period, this paper establishes a three-dimensional free surface capture high-resolution hydrodynamic model of the Bohai Sea and reports on simulations of the tide-current distribution in the Bohai Sea. Then, using the model as the background dynamic field and based on the open sea area oil spill model and thermodynamic processbased sea ice growth and melting model, this paper constructs a numerical prediction model for oil spills in the Bohai Sea during the winter sea ice period. The model is then verified, and the parameters calibrated by the measured values, fully demonstrating that the numerical model of oil spills during the winter ice period established in this paper is of high accuracy, stability and ubiquity.
where u (x, y, z, t), v(x, y, z, t), and w(x, y, z, t) are the flow velocity components in the horizontal x and y and vertical z directions, respectively; t is time, and p is the fluid pressure determined by the constant reference density. Not considering the atmospheric baroclinic and dynamic pressures, p=g(η-z), where η is the free surface rising value or tidal level, f is the Coriolis force coefficient, g is the gravitational acceleration, and γ h and γ v are the horizontal and vertical eddy viscosity coefficients, respectively. These coefficients are obtained by the dual equation turbulent flow model (Wang, 2009), also meeting the requirement to close the Reynolds-Average Navier-Stokes equation.
Mass conservation of incompressible fluid can be expressed by the continuity equation Integration of the continuity equation at the interval of [-h, η] between the bottom and free surface in combination with kinematic boundary conditions lead to the equation of the free surface capture and tidal level, presented as follows: where h(x, y) is the depth of still water, and H(x, y, t)=h(x, y)+η(x, y, t) is defined as the total depth of the water. The detailed elevation transform relation is shown in Fig. 1.

Boundary conditions
Boundary conditions of the free surface (z=η) are defined with wind effects, shown as: where u a and v a are the components of the wind speed in the x and y directions, the wind speed is ; γ T is the wind stress coefficient as determined by the wind speed.
Boundary conditions of the seabed surface (z=-h) are defined with the bottom friction as shown in Eq. (8): where γ B is the nonnegative bottom friction coefficient, which can be obtained through the turbulent boundary layer theory.

Thermodynamic model
Simulation results of the hydro-dynamical field were taken as the background field. In addition, sensible heat and latent heat flux were calculated by the bulk method, given the heat exchange between the sea ice and atmosphere. Additionally, the amount of solar radiation absorbed by the sea was calculated via an empirical formula. All three aspects above were integrated to simulate the growth and melting process of sea ice. The main heat budget forms of the sea surface include solar shortwave radiation entering the sea, Q sw , the net infrared radiation from the sea surface, Q lw , the sensible heat flux (transmitted from the sea to atmosphere), Q s , the latent heat flux (evaporated from the sea water), Q L , and the advection heat flux generated by heat transfer during ocean circulation, Q V .
The heat balance formula of air-sea interface is shown in Eq. (9): where Q is the total heat budget via the sea surface. Note that Q is zero when the global ocean is taken in its entirety in the long run but usually remains unequal to zero in short time series when only local sea areas are studied. Therefore, when Q>0, it indicates that the sea gained heat or vice versa. Absorption of solar radiation refers to the intake and reflection of solar radiation by sea water, which differ largely. Sea surface reflectivity of solar radiation ranges from 0.07 to 0.12, and most solar radiation is absorbed in that range. Additionally, the reflection coefficient of solar radiation for a frozen sea increases to approximately 0.6. The latent heat flux, Q L , and sensible heat flux, Q s , of the air-sea interface are collectively known as the turbulent heat flux of the airsea interface, which works as an important mechanism for the interaction between the air and sea to affect climate. This is not only a necessary parameter to represent underlying surface forcing and atmosphere interactions but is also a key factor affecting the oceanic mixed layer and seasonal thermocline changes. The sea exchanges moisture and heat with the atmospheric boundary layer in the form of latent and sensible heat to change the stability of the atmosphere and affect atmospheric circulation form. In this paper, solar radiation, sensible heat flux and latent heat flux are all considered for the heat budget.

Computational process for solar radiation
Solar radiation is the main energy imparted via the earth's atmosphere to the sea. Many scholars have conducted detailed analyses and studies of solar radiation and have proposed various empirical estimation methods for its calculation (Zhang et al., 2001;Wu and Xu, 2003;Liu and Yin, 2006;Ji and Yue, 2000). Here, it is calculated via the empirical formula (10) proposed by Sellers (1965): where Q * is the solar constant, which has a value of 1368 W/m 2 (Dobson and Smith, 1988), and A and B are the cloud scale coefficients determined by specific cloud scales. Cloud scales are divided into nine levels according to the WANG Kun et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 185-197 187 number of cloudy days in each month. Cloud scales and coefficients for the Bohai Sea in 2011-2012 were used in Eq. (10) and are presented in Table 1.
where is the local latitude, is the sun declination, and is the hour angle. changes little in a day and can be regarded as a constant. The sun declination is δ=arcsin[sin23.5 π180sin(2πKD)] (Wu and Xu, 2003), where D=365 or 366, and K is the number of days since March 21 in the year of calculation. The equation for calculating F 1 is integrated over one day to obtain the empirical equation for calculating the total amount of solar radiation days: (2ω 0 sin 2 φ sin 2 δ + 4 sin φ cos φ sin δ cos δ sin ω 0 + ω 0 cos 2 φ cos 2 δ + 0.5 cos 2 φ cos 2 δ sin 2ω As shown in Table 1, the cloud scale coefficients and cloud scales for the solar radiation of the Bohai Sea in a common ice year from October 1, 2011, to March 31, 2012, were calculated. The comparison between the solar radiation calculated with the selected cloud scale and the observed values of the average solar radiation of the Bohai Sea over the years provided by the website is shown in Table 2. It can be seen that the data fit well and are reliable.

Calculation process of oceanic heat flux
Three methods can be used to determine the oceanic heat flux, which are the eddy correlation, bulk and residual methods (Ji and Yue, 2000). In this work, the most convenient bulk method is adopted, which employs direct calculation through the water temperature, flow velocity and seawater salinity. The key here was to determine the heat transfer coefficients. The coefficient of the heat flux transmission between ice and water directly affected oceanic heat flux, the calculation of latent heat and sensible heat fluxes (Gao et al., 2003;Chen et al., 2006). The selection of the transmission coefficients for the two terms varied slightly with the number of hierarchies and dimensions in the simulation model, with the values ranging from 0.002 to 0.003. Determination of the values in the numerical simulation of the Bohai Sea ice considered sea ice characteristics, meteorology, and hydrology, and this paper accounted for the effects of the temperature, relative humidity, wind field characteristics, and the coefficients of sensible heat and latent heat fluxes on oceanic heat flux.
According to the formula of the bulk method (Lu, 1988), the sensible heat and latent heat fluxes are calculated by Eqs. (13) and (14), respectively: where Q H and Q E are sensible heat and latent heat fluxes, respectively, and C H and C E are the coefficients of sensible heat and latent heat fluxes. T a is the potential temperature of air near the sea surface and is calculated by Eq. (15): where z r is the observed altitude, which is usually 2 m. The temperature used was obtained from difference values among regions close to the Bohai Sea. For instance, hourly measured temperatures in coastal cities within the Bohai Rim Region (like Dalian and Yantai) and cities close to inland (like Shenyang) were taken as splashes for interpolation to obtain the temperatures at each grid node and calculate the real-time temperature field in the area, and the temperature measured points are shown in Fig. 2. T is the temperature at z r ; T s is the sea surface temperature; L is the latent heat of water evaporation; ρ a is air density, which is usually the moist air density given complex situations at the air-sea interface; C p is the air specific heat at constant pressure; V a , q a and q s are the wind speed, specific humidity and specific humidity of the ice surface at 10 m, respectively.
With regard to the selection of the wind speed value, the average wind speed in 10 years of 50 m above the sea surface was found on the NASA website, and the sea surface  wind speed showed a logarithmic distribution with height, as seen in Eq. (16): where z=50 m, K z =[ln(10/z 0 )/ln(z/z 0 )], and z 0 is the roughness height. When the wind speed of 50 m above the sea surface exceeds 7 m/s, z 0 =0.022 (Chen et al., 1989), and we obtain K 50 =0.79. When this wind speed is smaller than 7 m/s, z 0 =0.0023 (Chen et al., 1989), we obtain K 50 =0.84. We then obtain the wind speed 10 m above the sea surface according to the data found on the NASA website by Eq. (16). Control of sensible heat and latent heat fluxes was realized by adjusting the coefficients of sensible heat and latent heat flux. Table 3 shows the values of all the thermodynamic parameters used for the program calculation of the Bohai sea ice growth and melting model.
In Chinese offshore waters, sensible and latent heat reflected obvious seasonality in the spatial distribution, and their distributions changed with time. Generally, latent heat exchange has obvious seasonal changes and remains positive all year round, representing the heating of the atmosphere by the sea. Its value is also much higher than the airsea sensible heat exchange. Comparison found that the effects of relative humidity on latent heat were also quite significant. Excessive relative humidity could cause reductions in heat evaporated from seawater and thus lead to corresponding decreases in the thickness and icing range of sea ice. Relative humidity also took the average values over many years and months, as shown in Table 4.

Oil spill model
By taking the oil slick thickness d as the main unknown quantity, marine oil spill and diffusion can be described through the following equation: where d is the oil slick thickness, u i is the component of all flow velocities solved through the hydrodynamic equation, and the oil slick flow velocity on grid points in the oil spill model results from superposed effects of the water flow velocity and wind speed; γ oil is the diffusion coefficient of the oil slick; and S is the point source of the oil spill determined by the synthetic attenuation coefficient, flow rate at the oil spill point and oil concentration. The value of the synthetic attenuation coefficient for the Bohai Sea was smaller than that for a non-winter oil spill.

Equation discretization
The grid is the basis of discreteness. This paper adopted the Delaunay triangulation method, which is based on the principle of"Circumcircle Property"(CP) to divide the computational domain plane with body-fitted non-orthogonal triangular grids, and applied the direct layering method in the vertical direction. Therefore, the three-dimensional computational domain in fact adopted the tri-prism grid. Semiimplicit finite-difference discretization was undertaken for the momentum equation at the edge of grid cells. Finitevolume discretization was undertaken for the continuity equation at the unit centre. The definitions of the variables on the grid and specific discrete equation were adapted from previous work Liu et al., 2017)   WANG Kun et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 185-197 189 as the turbulent equation for discreteness (Wang, 2009).

Model verification
4.1 Verification of the hydrodynamic model The Bohai tidal model was used as the basis for the transport models of the sea ice and the oil spill in the sea, and the simulation results achieved reliable precision and ensured the simulation precision of subsequent models. To verify the accuracy of tides predicted by the hydrodynamic model, measured tidal levels at ports in Yantai were selected for comparison with simulated values using continuous 20-day and 10-day measurements from June 26, 2015, and June 01, 2015, at the Yantai Port and Bayuquan Station, respectively. The computational domain, grid distribution and results are shown in Figs. 3-5.
Since verification data for the flow rates in the Bohai Sea area have not yet been obtained, the flow rate was verified by the numerical simulation results of the Bohai Sea and the Yellow Sea large computational domain. Specifically, it is verified by the continuously measured sea current data for the three stations in this area for 25 hours starting at 12 o'clock on June 2, 2015. The station positions are shown in Fig. 6. The comparisons between the simulated and measured values are shown in Figs. 7-9.
Through the comparative verification of the abovemeasured tidal levels, sea current values and simulated values, the calculation results of the hydrodynamic model in this paper agree with the measured values, except for individual gaps that may mainly derive from the impact of the treatment of dynamic dry/wet boundaries and the consideration of open boundary tidal wave systems. However, in general, the calculation accuracy of the hydrodynamic model in this paper can meet the needs of the simulation calculations of sea ice growth and melting and oil spill models.        Fig. 11. Comparison of the ice thicknesses during the severe ice period.
WANG Kun et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 185-197 191 influence of dynamic factors such as sea ice collision. On the other hand, it could be due to the inaccuracy of the initial available thermodynamic factors. However, generally speaking, the accuracy of the simulation results can be guaranteed, and the sea ice growth and melting model in this paper can be used to simulate sea ice growth and melting in the Bohai Sea in serious ice years.

Light ice year
Taking the light ice year from 2006-2007 for example, the sea ice growth and melting condition was simulated with the sea ice model in this paper, the results of which are shown in Figs. 13-15. Compared with the monitoring data, the results are found to be basically in good agreement with the maximum deviation of approximately 6 days in the initial ice period, which indicates that the sea ice growth and melting model can also be used to simulate sea ice growth and melting in the Bohai Sea in light ice years.

Numerical simulation of Bohai oil spill during the winter sea ice period and result analysis
The winter oil spill model for the Bohai Sea surface was established by coupling the thermodynamic model and oil spill model on the basis of the hydrodynamic model. Most original data adopted relevant data in the Bohai Sea ice condition model. Topographic data were obtained from 1:10000 chart data, to which were added the oil spill transport quantity on the basis of the ice condition model, provided boundary conditions, process and relevant parameters of the oil spill during the winter ice period (mainly including the oil spill diffusion coefficient and the synthetic attenuation coefficient). The grid scale was the same as for the previous model, with a maximum scale of 5000 m, a minimum scale of 3000 m, totally 4994 nodes and 9648 units generated. All values took the initial value of 0 and roughness coefficient of 0.015. The time step was 0.05 h, there were 86880 counted steps, and the output interval was 480 steps (24 h). A total of 4344 h of the water levels and the oil spill process were simulated. For the Bohai oil spill model during the winter sea ice period, it could be found through calibration that the diffusion attenuation coefficients of 1.0×10 -7 m 2 /s and 1.3×10 -7 m 2 /s, respectively, were proper.
(1) Oil spill boundary Since the oil spill point source was in the middle of the Bohai Sea, no oil spill into or out of the boundary, thus, the boundary condition of the oil spill was zero.
(2) Oil spill point source The oil spill process was provided by referring to that of the Penglai 19-3 oil spill accident. It was assumed that there were two point sources, the positions of which are shown in Table 5. The first point source was provided with the flow  process of the oil spill at Platform B, where it continued from January 2 to January 17, 2012. The second point source was provided with the flow process of the oil spill at Platform C, where it continued from January 13 to March 13, 2012. The point source process is shown in Fig. 16. Simulation results of the ice conditions, oil spill field, and oil spill area, which vary with time, are shown in Figs. 17-26.
The distribution diagrams of the ice field and oil slick indicated that the initial oil spill volume and sea ice coverage were zero. From December 30, 2011, to January 3, 2012, freezing initially appeared in the Liaohe Estuary in the northern Bohai and gradually developed into the Bohai Sea. Subsequently, the area of the ice region was obvious and was enlarged in the Liaohe Estuary. Ice conditions began to occur in the Laizhou Bay, the ice thickness and area gradually increased, with the area reaching the peak on February 10, 2012, and then decreased gradually. As the area of the ice region in the Laizhou Bay declined to zero, yet a relatively large ice area continued in the Liaohe Estuary on approximately February 28, 2012. The sea ice in the Liaohe Estuary melted completely on approximately March 18, 2012, and the ice period for Bohai ended. The simulation results are in good agreement with the satellite imagery remote sensing data from the National Satellite Ocean Application Service. On January 3, 2012, a distributed oil slick appeared near the Laizhou Bay and diffused in all directions. Since the point sources were near the shore, the oil film arrived at the coastline southwest of the oil spill point on approximately January 11, 2012. Meanwhile, the forma-   WANG Kun et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 185-197 193 tion of sea ice impacted the motion of the oil film. With the increasing areas of ice and the oil film, the impact of sea ice on the oil film on approximately January 18, 2012, was obvious. The impacts of sea ice on the movement of the oil film lasted until February 28, 2012. Subsequently, the movement of the oil film shed the effects of sea ice due to the melting and disappearance of the latter. The oil film gradually moved towards the right side of the Laizhou Bay and decayed with the motions of the tides and the wind field. On March 31, 2012, part of the oil film remained along the Laizhou Bay. The time sequence of the distribution of the oil spill volume and area showed that the movement and change    process of the oil film were obviously affected by sea ice. Since the air temperature was relatively low in winter, the weathering velocity of the oil spill diminished, and the attenuation of the oil spill therefore was slower than that of the oil spill under the hypothetical ice-free condition from the same period. In winter, sea ice affected the movement of the oil film and decreased its diffusion velocity. Meanwhile, the weathering velocity of the oil spill in winter also declined, resulting in a longer retention time of the oil spill and a larger diffusion area of the oil film.

Conclusions
This paper adopted the three-dimensional free surface to   WANG Kun et al. China Ocean Eng., 2019, Vol. 33, No. 2, P. 185-197 195 establish the high-resolution hydrodynamic model and simulate tidal distribution in the Bohai Sea. By establishing an oil spill model for an open sea area and sea ice growth and melting model, the oil spill model for an open sea area and thermodynamic model were combined to establish the numerical model for the Bohai oil spill during the winter sea ice period. Measured values were used to verify the hydrodynamic model and calibrate the roughness and eddy viscosity coefficients to prove the relatively high precision, stability and resolution of the proposed numerical model. Finally, the most appropriate diffusion coefficient and comprehensive weathering attenuation coefficient of the Bohai oil spill during the winter sea ice period were determined. The model not only considered the sea ice growth and melting model for the temperature transport quantity but also coupled the ice and oil film field based on studying the oil spill model for the open sea area. This could reflect the effects of environmental conditions such as air temperature and sea ice on the motions of oil films and prediction of oil spill positions and pollution ranges to forecast accidental oil spills rapidly. In conclusion, this study provides reliable technical means for emergency studies on winter oil spill accidents.