Numerical Simulation of Ship Oil Spill in Arctic Icy Waters

Abstract

Featured Application

The trajectory and geometry of oil film could be predicted at a certain future moment to provide an effective method for the emergency treatment of maritime oil spill.

Abstract

This paper presents a three-dimensional numerical simulation model of an oil spill for application in emergency treatment methods under icy water conditions. The combined effects of wind, wave, current and ice implemented in our model correspond to Arctic Ocean conditions. A discrete element method combined with an overset grid was adopted to track the trajectory movements of oil film with medium-density ice floes and simulate the flow field of moving ice of large displacement in six degrees of freedom (6DOF). The probability of oil spill area extensions were estimated by a response surface method (RSM). Results showed reduced risk of pollution in icy water conditions and greater drift action of oil film. Accordingly, the spraying location and quantity of oil-dispersant could be rapidly specified.

1. Introduction

In recent years, ice-free areas in the Arctic Ocean increased especially during summer months. This increased business utilizing deep waters in this region [1]. With the opening of the Arctic Shipping Route, oil spill accidents from irregular operations or ship collisions in icy Arctic waters can cause tremendous damage to the aquatic ecosystem [2]. The oil spill is a serious pollution problem for the marine environment, as well as a concern for the survival quality of fish and algae [3], particularly in fragile Arctic regions. The hydrocarbons contained in petroleum can seriously damage the physiological and biochemical processes of marine organisms. Therefore, it is imperative to study emergency methods to treat ship oil spills. The diffusion and drift of oil films on the water surface, as well as the spread of oil under ice or trapped by ice floes, can be simulated to protect the Arctic maritime environment.

Oil-spill behaviors in ice-free waters have been researched systematically since the 1970s. In this work, the convection diffusion equation was considered suitable for the prediction of oil-spill diffusion under water [4]. In this paper a buoyant jet model based on the Lagrangian control volume method was established. The method considered both the diffusion and dissolution processes of oil spills [5]. Additionally, the “particle flow” method was used to calculate the transport positions of oil particles at each time step to simulate the transport process of the oil film and associated oil-spill trajectory [6]. The trajectory of oil-spill diffusion on the sea surface was presented visually using a discrete phase model (DPM) in FLUENT software [7]. The Sea Track Web system was applied for oil drift forecasting [8].

The forecasting of oil spills in water is complex, and the presence of ice in cold waters increases further the difficulty of this research. This is because of the fundamental interaction among different ice blocks that possess different equilibrium thickness in cold and warm open waters [9]. Depending on the natural geographic condition of the Bohai Sea, a dynamic sea ice model has been developed to understand the effects of various forces on sea ice movement [10]. Based on continuum theory, the viscous-plastic constitutive relationship of sea ice was analyzed to forecast the movement of an oil spill with different ice densities [11]. The morphology of “oil particles” was used to establish the drift-diffusion model of the oil spill [12]. Afterwards, the evaporation and emulsification processes were added to the forecasting of oil spills in the pack ice area to correct the drift equation [13]. A surface current model and a sea ice model were combined to simulate the oil slick movement process, including drifting, spreading, evaporation and emulsification [14]. However, in its current model, the waves on the sea surface were neglected, and only the oil spill on the ice surface could be forecasted. The movement of an oil spill under ice was observed with the help of 3-D under-ice imagery from an autonomous underwater vehicle [15]. Nevertheless, the spread of an oil spill under the ice was difficult to predict precisely.

The objective of this work has been to consider synthetically the real environmental elements of wind, current as well as waves in the Icelandic Sea of the Arctic Ocean with the aim to simulate an authentic environment flow field for the prediction of oil-spill trajectory. To achieve this, the multi-phase volume of fluid method (VOF) [16] and discrete element methods (DEM) [17] were adopted to track the non-linear free surface and establish a three-dimensional icy waters model. The Eulerian–Lagrangian method was applied to simulate the movement of the oil spill in ice-infested and ice-free waters. To shorten the forecasting time, the response surface method (RSM) [18] was used to fit the functional relationship between the oil-spill variables and various responses.

2. Modelling Methods

Based on the unique characteristics of the Arctic sea, the combined influences of wind, current, wave and sea ice were considered with the aim to establish a three-dimensional numerical model of a ship oil spill. The oil, air and water were separated into three insoluble phases. The air phase assumed steady wind and the water phase assumed combination of current and non-linear waves. The VOF was used to track the free surface of the air–water phase interface and DEM helped idealize the interactions between ice floes. A Eulerian–Lagrangian discrete method was applied to simulate the movement of oil. Accordingly, instead of static analysis, a six degrees of freedom (6DOF) dynamic simulation code was introduced to control the pitching motion and horizontal displacement of ices. The oil trajectory dynamic variation was achieved by coupling the movement of floating ice and oil drift.

2.1. VOF Method

The sea surface with wave was presumed as an air–water interface, and the oil spill from underwater to water surface was a three-phase oil–water–air flow. Therefore, the VOF was adopted to track the free surface of the water [19]. The VOF method tracked the fluid flow in each control cell by constructing a fluid volume fraction function to inform the free surface shape with its function and derivative values [20]. When the specified phase fluid was located at the lower left of the grid element, the lower left vertex was the datum, while the fluid was in another orientation, and the datum was obtained by coordinate flipping. The outward normal in a free surface grid was taken as the unit of the outward normal vector. The horizontal direction was the X axis and the vertical direction was the Y axis. The fluid volume fraction function $F_q$ satisfied Equation (1) [21]. If $F_q = 1$, the entire cell consisted of the q-th phase fluid; if $F_q = 0$, there was no q-th phase fluid in the cell; and if $0<F_q<1$, the cell was considered as an interface cell. $F_q$ was a continuous function of time and space, and the sum of the volume ratios $\sum F_q$ of each phase q was 1. In the VOF method, the physical parameter $\phi$ was calculated by Equation (2).

(1)$\frac{\partial F_{\mathrm{q}}}{\partial t}+\frac{\partial \left(u{F}_q\right)}{\partial x}+\frac{\partial \left(v{F}_q\right)}{\partial y}=0,\sum F_q = 1$

(2)$\phi =\sum {\phi }_q{F}_q$

2.2. Discrete Element Method (DEM)

The interactions among sea ice cells were modeled by DEM. In the normal direction, the cells in DEM can be regarded as being connected by a damper and a spring in parallel. In the tangential direction, they could be assumed connected by a spring, a damper and a sliding friction device [22]. During the interaction process, the viscoelastic forces among cells resulting from their relative velocities and elastic deformations were considered, and the shear force was computed according to Mohr–Coulomb [23] friction law. The normal force between two ice cells was given by Equation (3). The tangential contact force between oceanic ices was expressed by Equation (4), and $K_{\mathrm{t}}$ was set to 60% of normal stiffness.

(3)$F_{\mathrm{n}}^n = K_{\mathrm{n}}{\delta }_{ij}-C_{\mathrm{n}}{V}_{ij}·n$

(4)$F_{\mathrm{t}}^n = \min \left(\left(F_{\mathrm{t}}^{n-1}-K_{\mathrm{t}}\Delta t\left(V_{ij}·t\right)\right),\mathsf{\mu }{F}_{\mathrm{n}}^n\right)$

2.3. Eulerian–Lagrangian Method

The diffusion and drift of oil film were computed by the conservation law of mass and momentum based on the convection–diffusion equation to simplify the coupling of the particle motion, chemical dynamic and hydrodynamic equations [24].

The ocean transport convection–diffusion equation was expressed in Lagrangian format as Equations (5)–(7). Differential dispersion was expressed by Equations (5)–(8) [25].

(5)$\frac{DC}{D\mathrm{t}} = \frac{\partial }{\partial z}\left(k\frac{\partial C}{\partial z}\right)+K_H{\nabla }^{2}C$

(6)${\nabla }^2 = \frac{\partial }{\partial {\mathrm{x}}^2}+\frac{\partial }{\partial {\mathrm{y}}^2}$

(7)$\frac{D}{Dt} = \frac{\partial }{\partial t}+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z}$

(8)$\frac{DC}{D\mathrm{t}} = \frac{C_Q^{K+1}-C_P^K}{\Delta t}$

After the time $\Delta t$, this micelle arrived at position $Q$, and position $P$ could be calculated by the velocity field along the trajectory of the substance micelle [26].

The amount of oil spill was expressed as a release of particles. Thus, when the thickness of an oil film on the water surface reached a terminal value (the mass balance of oil spill), the oil layer could be regarded as a spatial mapping distribution of particles. Slick thickness $h$ in grid cells was calculated by Equation (9).

(9)$h =\sum \frac{V_i}{\Delta A}$

The Euler–Lagrangian method was then used to establish a drift–diffusion model of oil spill on the seawater surface. The equation of motion was expressed as Equation (10). The diffusion coefficient of the oil film was expressed as Equation (11).

(10)$\frac{\partial h}{\partial t}+\frac{\partial h\left(u+{\tau }_x/f\right)}{\partial x}+\frac{\partial h\left(v+{\tau }_y/f\right)}{\partial y} = \frac{{\partial }^2{D}_{s}h}{\partial x^2}+\frac{{\partial }^2{D}_{s}h}{\partial y^2}+R_h$

(11)$D_{\mathrm{s}} = g{h}^2\left({\rho }_w-{\rho }_o\right){\rho }_o/{\rho }_{w}f$

3. Numerical Simulation of Oil Spill

3.1. Simulation of Oil-Spill Drift

The wind and current were assumed to affect the centroid of the oil film. The Eulerian–Lagrangian tracking method was adopted to simulate the drift trajectory of the oil film centroid on the water surface as Equation (12) which explained the initial position $S_0$ of the oil film centroid drifting to the new position $S$ after a time $\Delta t$. $U_L$ was calculated by Equation (13) and $a_w$ was commonly ranging from 0.03 to 0.04. $S$ was expressed as Equations (14) and (15). The velocity vector of oil-spill drift is shown in Figure 1.

(12)$S =S_0+{{\displaystyle \int }}_0^{0+\Delta t}{U}_{L}dt$

(13)$U_L=a_f{U}_f+a_w{U}_w$

(14)$x^{,}=x_0+u\Delta t+i_w·U_w·sin\theta$

(15)$y^{\prime }=y_0+v\Delta t+i_w·U_w·cos\theta$

3.2. Simulation of Ice-Floes Drift

The short-time prediction model of ice-floes movement was based on momentum conservation and mass continuity equations, therein, the critical step was to determine the constitutive relation of the ice internal force. The constitutive relation of the ice floes could be divided into two categories, based on the continuum and discrete medium. According to different ice densities, the viscoplastic constitution was adopted for continuous ice floes, and the collision-type constitutive relation was employed for the edge ice belts with discrete medium properties. Two types of equations were combined for the hydrodynamic model of ice floes, as shown in Equations (16)–(18). Equation (16) presented the momentum balance equation. Equations (17) and (18) presented mass conservation; and $diffusion$ represented the term of energy diffusion. In the prediction model, the Coriolis force, the ocean surface slope and thermodynamic terms were assumed negligible.

(16)$M\frac{D\overrightarrow{U}}{Dt} = -Mf\overrightarrow{k}\times \overrightarrow{U}+A\overrightarrow{{\tau }_a}+A\overrightarrow{{\tau }_w}-Mg·gradH+\overrightarrow{F}$

(17)$\frac{\partial h_a}{\partial t}+\overrightarrow{U}·\nabla h_a = {\Phi }_h+diffusion$

(18)$\frac{\partial A}{\partial t}+\overrightarrow{U}·\nabla A = {\Phi }_A+diffusion$

3.3. Numerical Model

A continuous oil-spill accident, which happened in the geographical range of 66°24′ N to 66°27′ N and 15°42′ W to 15°54′ W, was simulated here. It was assumed as a static point source in the environmental conditions of waters near Iceland. As shown in Figure 2, the location of the oil spill was identified as Point A (66.43° N, 15.79° W). The parameters used in this model are shown in Table 1. Some parameters on the behavior of oil refer to the obtained marine characteristics [27] and the sensitivity studies by Ross and Dickins [28].

The computational domain (Figure 3) had a dimension of 2000 × 1000 × 500 m, with the air domain on the top, water domain on the bottom, and the free surface at the height of 400 m. The smallest grid cell was 0.005 m, and there were 1.35 million cells in total. The oil-spill spout was set as 2 m in diameter. The oil spurted vertically from 0.5 m underwater, and was located at 100 m to the left end of the domain. The left side of the domain was the speed inlet; the right side was the pressure outlet; the upper and lower boundaries were wave-eliminating.

The grids were locally refined around the free surface and the oil spout. The girds refinement aimed to capture the flowing details preferably without significantly increasing the amount of calculation, especially in the small domain near the free surface and the oil spout (to guarantee grids with more than 20 layers in the wave height direction, 80 grids at one wavelength, surface grids refinement near the oil spout of 25%), as shown in Figure 4, and a three-dimensional implicit unsteady model of the marine environment was established. The dominating equations involved in the numerical calculation were the fluid continuity equation, mass conservation equation, and momentum equation (Navier–Stokes equation) [29]. The VOF method was adopted to describe the three-phase oil–water–air flow. The waves were fifth-order VOF waves with a wavelength of 60 m and a wave height of 1.5 m, which were initialized and shown in Figure 5. The drift of ice floes was simulated by DEM. An ice ejector was constructed using Lagrangian discrete particles, and each ice floe was simulated as a round, flat particle with a thickness of 0.2 m [30]. It was presumed that there were mutual collisions between ice floes without fragmentation and the particles were randomly sprayed in accordance with the derivative components of grid arrangement. After the ejector generated ice floe particles for 100 s, the oil-spill accident was initiated. This time point could be approximated as the completion of the ice-field initialization. The $k-\epsilon$ model was selected to present the turbulence mode and the computational step was set to 1 s.

A series of meshes ranging from 10.2 to 14.6 million elements were made as an example case, where 10.2, 12.2 and 14.6 million elements are designated as the coarse mesh, the medium mesh and the fine mesh, respectively. The computational time of the medium mesh was approximately 80 h, which was 3 times faster than that of the fine mesh using the 64-bit operating computational system, 3.8 GHz, 32 Core Processors with 8 GB RAM. The computational time of the coarse mesh was approximately 50 h with the same computational parameters. Thus, the medium mesh parameters were considered appropriate for this study.

3.4. Characteristics of Oil Spill in Icy and Ice-Free Waters

The oil spill behavior was simulated for both open (ice-free) and icy water conditions under the same sea state. The oil diffusion trajectories for open waters and icy waters are shown in Figure 6 and Figure 7, respectively. The blue color represents the oil film and the other colors represent the ice floes on the sea surface. In icy waters, the extent of the oil spill was obviously suppressed and its distribution was relatively narrow. Its maximum drift distance and maximum diffusion area radius were much smaller than those in open waters. At the time of 500 s, in open waters, the maximum drift distance of the oil spill was 1.2 km in the horizontal direction and 0.6 km in the vertical direction; in icy waters, the maximum oil-drift distance was only 0.7 km in the horizontal direction and 0.1 km in the vertical direction. The main direction of oil-spill diffusion was affected by wind and waves under both conditions. In ice-free waters, the oil film was widely spread on the sea surface, whereas in icy waters, a band-like distribution was observed.

Although the presence of ice was observed slowing down the oil diffusion rate, the blending of ice and oil would increase the difficulty of oil recovery. The response measures for the oil spill should, therefore, be adopted in a timely manner according to the oil-spill trajectory. In the forecast area of the oil spill, the enclosure and control facilities should be deployed to prevent the wide-area oil-spill accidents.

In icy waters, the oil diffusion was affected by wave action on the free surface. These diffusion trajectories are shown in Figure 8 and Figure 9, where the red portion represents the seawater, the blue portion represents air and the interface between the two is the free surface. The black parts indicated the ice floe particles that rose and fell with the waves. The oil diffusion was represented by contours and the oil spout was located at y = 0 m. Oil drift could be observed along the horizontal direction with time and the thickness of the oil film changed continuously, representing a simultaneous process of drift and diffusion.

The width and thickness of the oil film are shown in Table 2. Note that the oil particles presented not only on the free surface but also below-the-water surface. This portion of oil gradually rose within the water until contacting the ice, and then spread to the surrounding water. The drift of the oil spill was evident in the simulation and the drift distance changed significantly. However, the change in the thickness of the oil film was not prominent. Therefore, it was concluded that the expansion mode of an oil spill in medium icy waters was dominated by drift motion.

3.5. Interaction of Ice and Oil

Considering current, wind and waves, the movement of ice floes with 6DOF included rolling, pitching and heaving motions on the water surface. The overset grid [31] technique was used to simulate the flow field of moving ice with large displacement. This technique simplified the construction of computational grids about complex geometries and allowed relative motion between embedded grids. Based on interpolating connections in every step in the time direction, at the state of ice floes flow associated with the non-linear wave, the background mesh established in the static domain was adopted to simulate the far field in water and the valid configurations of overset mesh were refined around boundaries consisting of ice particles.

In the oil-spill model, oil density and ice diameter were set at 840 kg/m³ and 1.5 m, respectively. As shown in Figure 10 and Figure 11, there were two main movements: (a) oil trapped under the ice floe while it rose to the surface and encountered the ice; (b) oil trapped by slush or brash (broken) ice and diffused to the surface of ice, which was consistent with the features of oil spreading in medium ice-density waters [9,11]. At the plane section of y = 0 m and y = 0.25 m, the oil drifted to 4.06 m at 3 s, and the horizontal and vertical displacement of the ice flow in the vicinity of the water surface were 1.61 m and 0.19 m, respectively. The pitch amplitude only increased from 0 to 6.0 deg. This is a small increase mainly affected by currents and waves. Due to the fact that the velocity in the y direction was limited to 10% of the wind velocity, the roll motion was completely negligible.

4. Prediction of Oil-Spill Pollution Area

Once an oil-spill accident has occurred in icy waters, it is nearly impossible to clean up the oil by mechanical technology. Thereby, a low-temperature oil-dispersant could be adopted as the treatment method. The diffusion process of oil dispersants could be negligible because they attach themselves onto the surface of oil film [32]. However, the presence of fast-moving or drifting ice in Arctic waters will affect both oil and oil-dispersants irregularly during the spreading process. Therefore, it will be difficult to control the dosage of oil dispersants and might lead to severe damage to the Arctic Ocean environment. Accordingly, the changing extent of the oil-contaminated area in icy waters should be predicted as soon as possible in order to introduce effective emergency treatment.

For this purpose, approximate modeling technology utilizing a regression analysis to fit factors and responses was adopted. The RSM was used to establish an approximate model to depict the relationship between the responses for its benefits of simple mathematical expression, small calculation amount and fast convergence speed. The response surface in the form of polynomial function was used. The polynomial function response surface was verified with good fitting accuracy and calculation efficiency for the input–output relationship of most data structures. The relationship between factors and responses in the approximate modeling could be described by Equation (19).

(19)$\widehat{y}\left(x\right) = a+{\displaystyle \sum }_{i = 1}^{n^{\prime }}{b}_i{x}_i+{\displaystyle \sum }_{i = 1}^{n^{\prime }}{c}_i{x}_i^2+{\displaystyle \sum }_{i = 1}^{n^{\prime }-1}{\displaystyle \sum }_{j>i}^{n^{\prime }}{d}_{ij}{x}_i{x}_j+\epsilon$

The polynomial function was adopted by RSM to fit the approximate model of the design space. The procedure for constructing the approximate model of the oil spill area used a third-order RSM, which was listed as follows:

• 1.. The sample point $x_1$ was determined to construct the approximate model in the design space using the experimental design methodology, where ${\mathrm{x}}_i = \left(x_1,x_2,x_3,\dots ,x_m\right),$ and m was the number of design variables. The data of sample points were the centroid coordinates (X, Y), area S, perimeter C’, horizontal maximum distance L and vertical maximum distance D of the contaminated area, which were extracted from the numerical results of oil spill diffusion at different times as presented in Section 3.4 (i.e., 40 s, 80 s, 120 s, 160 s, and 200 s).

• 2.. A series of sample pairs of designated variables and responses $\left\{\left(x_i,y_i\right)\right\}$ were computed by the model, where ${\mathrm{y}}_i = \left(y_1,y_2,y_3,\dots ,y_n\right)$ was the response of the oil-spill area, and n was the number of sample points designed for the test.

• 3.. The regression analysis using RSM was conducted to obtain an approximate model of the oil-spill area. The multiple correlation coefficient was adopted to evaluate the credibility of the response surface model. If the requirements were met, this model could be applied to the subsequent optimization computation. Otherwise, the credibility of the model could be improved by changing the model type and increasing the sample points.

The oil-spill area at a certain future time was predicted by the function that described the relationship between variables and responses. The appropriate equivalent weight of oil dispersants could be determined based on the oil-spill area to implement the emergency treatment. In order to shorten the calculation time, the prediction of the oil-spill pollution area was used as the continuation and supplement of the oil trajectory simulation. Accordingly, the response surface model could be established by using source data of the oil spill, which were attained from the numerical calculation results of Section 3 shown in Table 3. In order to improve the accuracy of fitting function, the functional relationship among t, $\left(\mathrm{X},\mathrm{Y}\right)$, C’, L, and D was correlated by the fourth-order and sixth-order models. The design variable was t and the responses were $\mathrm{X}$, $\mathrm{Y}$, C’, L, and D. The fourth-order and sixth-order models are shown in Equation (20).

(20)$\begin{array}{c}\mathrm{x}=0.00000003t^4+0.0000006t^3-0.0027t^2+0.8342t+94.618\\ \mathrm{y}=-0.00000000005t^6+0.00000003t^5-0.000008t^4+0.0009t^3-0.0498t^2\\ +1.1752\mathrm{t}-4.6404\\ \mathrm{L}=0.0000003{\mathrm{t}}^4-0.0001t^3+0.0077t^2+1.2593t+97.143\\ \mathrm{D}=-0.0000001{\mathrm{t}}^4+0.00002t^3-0.0291\mathrm{t}+9.6553\\ {\mathrm{C}}^{\prime }=-0.000000002t^6+0.000001t^5-0.0003t^4+0.0243t^3-1.0821t^2+21.627\mathrm{t}\\ -67.344\end{array}$

The minimum number of sample points depended on the model order and the number of input variables. A third-order response surface model could be constructed for the oil-spill area S, for which the input variables were the perimeter C’, centroid position (X, Y), horizontal maximum distance L and maximum vertical distance D. Thus, the number of variables M was 5. The minimum number of sample points must be equal to $\left(M+1\right)\left(M+2\right)/2+M$, which was 26 in this research. Therefore, the fitting equation of the oil-spill area S is shown in Equation (21).

(21)$\begin{array}{cc}\mathrm{S}=165789.84& +109.3{\mathrm{C}}^{\prime }-4672.69\mathrm{X}-6503.4\mathrm{Y}+1138.7\mathrm{L}+973.86\mathrm{D}+0.18{C^{\prime }}^2+36.87X^2\\ & -349.34Y^2+1.26L^2-5.86D^2+0.1C^{\prime }\ast X-14.07C^{\prime }\ast Y-0.7C^{\prime }\ast L\\ & +0.72C^{\prime }\ast D+72.93X\ast Y-15.32X\ast L -10.34X\ast D-0.16Y\ast L+37.99Y\\ & \ast D+2.34L\ast D-4.12{C^{\prime }}^3-0.05X^3+26.73Y^3+0.01L^3-0.02D^3\end{array}$

The space model of factors and responses is shown in Figure 12 and Figure 13. The local effects and global effects for responses were shown in Figure 14 and Figure 15, respectively.

Obviously, factor $\mathrm{Y}$ played a dominant role in both local effects and global effects for $\mathrm{S}$. Therefore, it was essential to predict the oil-spill area by designated centroid coordinate $\mathrm{Y}$ in icy waters. Emergency treatment measures for the oil spill could then be organized without considering non-sensitivity factors $L$ and $C^{\prime }$.

The reliability of this approximate model was evaluated by the multiple correlation coefficient $R^2$ as given in Equation (22):

(22)$R^2=1-{\displaystyle \sum }_{i =1}^{n^{″}}{\left(y_i-\widehat{y_i}\right)}^2/{\displaystyle \sum }_{i =1}^{n^{″}}{\left(y_i-\overline{y_i}\right)}^2$

The value of $R^2$ was used to measure the degree of conformity between the approximate model and the sample points. $R^2$ being close to 1.00 indicated that the approximate model had a relatively high degree of credibility. In this case, $R^2$ was 0.9 and achieved the established threshold, which meant the approximate model could provide a relatively credible prediction of oil-spill diffusion.

5. Conclusions

In this study, a three-dimensional DEM-based numerical model of ship oil-spill movement in Arctic icy waters was established. The trajectory model of oil film coupled with marine loading including wind, current and waves utilized the multiple tracking method of oil drift and ice floe drift. The accident scenario was simulated with reference to a real experiment. The interaction details of oil and ice floes were expanded by a 6DOF model to describe the features of oil spreading. Subsequently, the correlation of parameters on the behaviour of oil in cold water were examined by sensitivity analysis, and the area of the oil spill was predicted using an RSM model in order to conduct rapid emergency treatment. The RSM model was evaluated by reliability analysis to indicate a high degree of credibility.

This model was presumably based on the unique characteristics of the Arctic Ocean and combined effects of wind, wave, current and ice. The wind was regarded as a steady wind of winter; the wave was simulated as a fifth-order Stokes wave in non-linear motion over time; the current was the relatively stable polar snap coupled with the wave; the ice floe was sprayed for 100 s before the oil’s spraying and stopped immediately when the oil started to spray. The oil spill was assumed as a continuous spill in a static point and the amount of oil spill was assumed to be 3 tons. The following conclusions were drawn:

• 1.. The presence of ice floes had a significant suppression effect on the movement of the oil spill. The maximum drift distance of oil film and the maximum radius of oil diffusion in icy waters were found to be much smaller than those in open waters. With the combined action of wind, wave, current and ice, the oil film was observed to expand to a larger area in ice-free waters but form a long and narrow belt in icy waters.

• 2.. The oil film diffusion and drift processes occurred simultaneously. The oil film in waters with medium-density ice was observed to be distributed on both the water surface and underwater. The underwater oil gradually rose until it came into contact with the ice, and then spread to the surroundings. The changes in the horizontal drift distance of the oil spill were more evident than those in oil-film thickness. Thus, oil-film drift constituted the primary behavior of oil-spill expansion.

• 3.. The regression function of time t, the centroid coordinates (X, Y), perimeter C’, horizontal maximum distance L and vertical maximum distance D were generated on the base of the approximate model method. This function exhibited a relatively high degree of credibility. It could predict the extent of oil spill in icy waters scientifically and reasonably, which could provide support for establishing a rapid emergency oil-spill cleanup system in icy waters, such as those found in the increasingly accessible Arctic Ocean.

Author Contributions

Formal analysis, W.L., X.L. and Z.D.; Methodology, W.L., X.L. and J.L. (Jianguo Lin); Supervision, J.L. (Jianguo Lin); Writing—original draft, W.L., P.G. and Q.M.; Writing—review and editing, J.L. (Jiamin Liu), H.W. and Z.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Engineering Laboratory for Mobile Source Emission Control Technology (NELMS2018A06) and Fundamental Research Funds for the Central Universities of China (3132017097, 3132018177, 3132019304).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

DPM	Discrete Phase Model
DEM	Discrete Element Method
RSM	Response Surface Method
VOF	Volume of Fluid Method
6DOF	Six Degrees of Freedom

Notations

$$u$$	Velocity component of fluid in x direction
$$v$$	Velocity component of fluid in y direction
$$\mathrm{t}$$	Flow time
$$F_q$$	Ratio of the volume of q-th phase fluid in the cell to the total cell volume
$${\phi }_q$$	Fluid parameter
$$F_{\mathrm{n}}^n$$	Normal contact force at time $$n$$
$$K_{\mathrm{n}}$$	Normal contact stiffness
$$C_{\mathrm{n}}$$	Contact viscous coefficient
$${\delta }_{ij}$$	Overlap value between the two disks $$i$$ and $$j$$
$$V_{ij}$$	Relative velocity vector of the contact faces
$$n$$	Normal unit vector of the contact cells
$$F_{\mathrm{t}}^n$$	Tangential forces at time $$n$$
$$F_{\mathrm{t}}^{n-1}$$	Tangential forces at time $$n-1$$
$$K_{\mathrm{t}}$$	Tangential stiffness between particles
$$\Delta t$$	Time step
$$t$$	Tangential unit vector between contact cells
$$\mu$$	Coefficient of sliding friction
$$C$$	Concentration
$$k$$	Vertical diffusion coefficient
$$K_H$$	Horizontal diffusion coefficient
$$P$$	Position of the marked substance micelle at time $$\Delta t$$
$$K$$	Normal coordinates
$$\Delta A$$	Number of cells in the computational domain
$$V_i$$	Volume occupied by the particles on the ocean surface in grid cell
$${\tau }_x/f$$	Shear stresses of the wind in the $$x$$ direction
$${\tau }_y/f$$	Shear stresses of the wind in the $$y$$ direction
$$D_{\mathrm{s}}$$	Diffusion coefficient of the oil film
$$f$$	Friction coefficient of the oil–water interface
$$g$$	Acceleration due to gravity
$$R_h$$	Physical–chemical kinetic terms
$$S_0$$	Initial position of the oil-film centroid
$$S$$	New position of the oil-film centroid after a time $$\Delta t$$
$$U_L$$	Velocity of oil-film drift
$$U_f$$	Velocity on water surface
$$U_w$$	Velocity of wind
$$a_f$$	Horizontal flow coefficient
$$a_w$$	Wind frag coefficient
$$x_0$$,$$y_0$$	Initial centroid position
$$x^{,}$$,$$y^{,}$$	Centroid position after a time $$\Delta t$$
$$u$$	Current velocity in the $$x$$ direction
$$v$$	Current velocity in the $$y$$ direction
$$i_w$$	Correction factor
$$\theta$$	Wind angle
$$Mf\overrightarrow{k}\times \overrightarrow{U}$$	Coriolis force
$$A\overrightarrow{{\tau }_a}$$	Force of the wind above ice floe
$$A\overrightarrow{{\tau }_w}$$	Force of the current above ice floe
$$Mg·gradH$$	Effects of the ocean surface slope
$$\overrightarrow{F}$$	Non-linear ice internal force
$$H$$	Gravitational potential on the ocean surface
$$f$$	Coriolis force coefficient
$$\overrightarrow{k}$$	Unit vector perpendicular to the ocean surface
$$\overrightarrow{U}$$	Velocity of the ice floe
$$A$$	Ice density
$$h_a$$	Average ice thickness in the grid
$$M$$	Unit area mass of ice floe
$${\Phi }_h$$	Thermodynamic term determined by ice thickness
$${\Phi }_A$$	Thermodynamic term determined by ice density
$$\mathrm{y}\left(x\right)$$	Response of a sample point
$$\tilde{y}\left(x\right)$$	Response of an approximate model
$$x_i,x_j$$	Design parameter value
$$n^{\prime }$$	Number of designated parameters;
$$\epsilon$$	Error between the response value and true value
S	Area of oil film
C’	Perimeter of oil film
L	Horizontal maximum distance of oil film
D	Vertical maximum distance of oil film
$$n^{″}$$	Sample numbers
$$y_i$$	Real value of simulation program
$$\widehat{y_i}$$	Estimated value of response surface model
$$\overline{y_i}$$	Mean of true response

Figures and Tables

Figure 1. Velocity vector synthesis of oil-spill drift.

Figure 2. Computational domain of Icelandic waters.

Figure 3. Computational domain.

Figure 4. Grids’ refinement around the oil spout.

Figure 5. Initialization of free surface wave.

Figure 6. Oil diffusion trajectories for open waters.

Figure 6. Oil diffusion trajectories for open waters.

Figure 7. Oil diffusion trajectories for icy waters.

Figure 8. Three-phase distribution of oil–water–air for icy waters at the oil emerging time.

Figure 9. Three phase distribution of oil–water–air for icy waters at 10.14 s.

Figure 9. Three phase distribution of oil–water–air for icy waters at 10.14 s.

Figure 10. Oil-spreading characteristics of the ice particle with six-degrees-of-freedom (6DOF) movement at y = 0 m plane section.

Figure 11. Oil spreading characteristics of the ice particle with 6DOF movement at y = 0.25 m plane section.

Figure 11. Oil spreading characteristics of the ice particle with 6DOF movement at y = 0.25 m plane section.

Figure 12. Three-dimensional space model of factors (X, Y) and responses (S).

Figure 13. Three-dimensional space model of factors (L, D) and responses (S).

Figure 14. Local effects (gradients) for S.

Figure 15. Global effects (average gradients) for S.

Table 1

The parameters used in the oil-spill model.

Properties	Parameters
Location	66.43° N, 15.79° W
Release situation	Accidental ship leakage
Spill rate	15 kg/s
Spill duration	200 s
Oil density	840 kg/m3
Oil viscosity	500 mPa·s
seawater density	1025 kg/m3
ice density	910 kg/m3
Ice diameter	5 m
Ice concentration	Medium (30–80%)

Table 2

Numerical results of drift and diffusion of oil film.

Position of y/m	Time(t)/s	Width of the Oil Film/m	Thickness of the Oil Film/m
y = 5	2.34–10.14	2.67–6.46	1.82–1.99
y = 10	5.64–10.14	2.96–9.51	1.44–1.94
y = 15	9.48–10.14	2.72–4.10	1.23–1.34

Table 3

Oil-spill parameters extracted from oil-spill trajectory simulation.

Time(t) /s	$$\mathbf{Area}\left(\mathbf{S}\right)/{\mathbf{m}}^2$$	$$\mathbf{Perimeter}\left(\mathbf{C}\text{’}\right)/\mathbf{m}$$	$$\mathbf{Centroid}\mathbf{Centroid}\left(\mathbf{X}\right)/\mathbf{m}$$	$$\mathbf{Coordinate}\mathbf{Coordinate}\left(\mathbf{Y}\right)/\mathbf{m}$$	$$\mathbf{Horizontal}\mathbf{Maximum}\mathbf{Distance}\left(\mathbf{L}\right)/\mathbf{m}$$	$$\mathbf{Vertical}\mathbf{Maximum}\mathbf{Distance}\left(\mathbf{D}\right)/\mathbf{m}$$
5	65.33	40.88	103.99	3.09	108.06	7.28
10	100.23	41.6	104.01	3.14	111.91	7.5
20	214.91	63.75	109.16	4.5	121.7	12.05
25	331.59	86.64	113.91	4.93	130.03	13.23
30	354.71	95.52	116.56	4.46	135.42	11.73
40	476.29	121.62	123	4.31	149.24	11.76
50	500.27	148.83	128.49	2.75	162.32	8.27
60	524.14	186.14	133.57	2.56	181.61	9.08
70	772.1	236.14	143.35	2.71	207.39	10.28
75	775.34	254.93	145.48	3.61	208.27	11.01
80	1138.79	302.18	154.09	4.98	224.39	15.93
90	1173.98	309.93	145.95	5.33	220.1	17.39
95	1262.44	346.96	147.97	6.16	220.89	21.58
100	1535.68	351.8	152.47	7.63	222.3	28.44
110	1618.11	408.45	151.36	7.93	222.72	30.97
120	2236.59	484.88	162.31	9.24	242.96	35.7
125	2291.43	454.33	168.08	8.13	251.08	32.12
130	2369.06	474.85	168.84	7.27	256.49	29.33
140	2423.48	432.39	176.93	6.63	266.22	24.62
150	2819.96	489.75	181.67	6.85	285.7	29.66
160	3344.54	690.65	173.01	11.14	270.35	55.5
170	3415.69	703.94	189.05	10.32	288	38.41
175	3674.69	883.42	190.71	11.51	306.49	45.38
180	3943.54	987.46	194.54	12.25	308.54	44.42
190	4225.53	942.32	202.43	10.35	338.16	41.04
195	4361.73	727.68	200.25	11.74	315.29	43.76
200	4792.95	707.69	207.89	10.7	340.31	41.7

References

1. Pietri, D.; Soule, A.B.; Kershner, J.; Soles, P.; Sullivan, M. The arctic shipping and environmental management agreement: A regime for marine pollution. Coast. Manag.; 2008; 36, pp. 508-523. [DOI: https://dx.doi.org/10.1080/08920750802398735]

2. Eide, M.S.; Endresen, O.; Breivik, O.; Brude, O.W.; Ellingsen, I.H.; Røang, K.; Hauge, J.; Brett, P.O. Prevention of oil spill from shipping by modelling of dynamic risk. Mar. Pollut. Bull.; 2007; 54, pp. 1619-1633. [DOI: https://dx.doi.org/10.1016/j.marpolbul.2007.06.013]

3. Shiau, B.S.; Tsai, R.S. A numerical simulation of oil spill spreading on the coastal waters. J. Chin. Inst. Eng.; 1994; 17, pp. 473-484. [DOI: https://dx.doi.org/10.1080/02533839.1994.9677614]

4. Johansen, Ø. DeepBlow—A Lagrangian Plume Model for Deep Water Blowouts. Spill Sci. Technol. Bull.; 2000; 6, pp. 103-111. [DOI: https://dx.doi.org/10.1016/S1353-2561(00)00042-6]

5. Yapa, P.D.; Li, Z. Simulation of oil spills from underwater accidents I: Model development. J. Hydraul. Res.; 1998; 35, pp. 673-688. [DOI: https://dx.doi.org/10.1080/00221689709498401]

6. Wang, P.; Li, Z.J. A Primary Numerical Model of Oil Spilling in Icing Sea. J. Glaciol. Geocryol.; 2003; 25, pp. 334-337.

7. Chen, X.; Yu, J.X.; Li, Z.G.; Wu, C.H.; Jiang, M.R.; Yang, Z.L. Application of Discrete Phase Model and Empirical Model in Oil Spill Simulation of Sea Surface. China Offshore Platf.; 2019; 34, pp. 47-53.

8. Ambjorn, C. Seatrack Web, Forecasts of Oil Spills, a New Version. Environ. Res. Eng. Manag.; 2007; 41, pp. 60-66.

9. Venkatesh, S.; El-Tahan, H.; Comfort, G.; Abdelnour, R. Modelling the behaviour of oil spills in ice-infested waters. Atmos. Ocean; 1990; 28, pp. 303-329. [DOI: https://dx.doi.org/10.1080/07055900.1990.9649380]

10. Wang, R.S.; Liu, Q.Z.; Chen, W.B.; Jin, H.T.; Liu, C.H. Numerical simulation and tests of sea ice drift process in the bohai sea. Oceanol. Limnol. Sin.; 1994; 3, 10.

11. Xia, D.W.; Xu, J.Z. Investigation on numerical modelling of spill oil movement in sea ice infested waters. Acta Oceanol. Sin.; 1998; 1, pp. 113-122.

12. Huang, J.; Cao, Y.J.; Gao, S. Numerical simulation of oil spill drift-diffusion in the Bohai sea. Mar. Sci.; 2014; 38, pp. 100-107.

13. Huang, Y.; Song, M.R.; Guan, P. Study on the Varying Regularity of Oil Spill Properties During the Ice Period of the Bohai Sea. J. Ocean Technol.; 2016; 35, pp. 1-5.

14. Yu, J.A.; Zhang, B.; Liu, Q.Z.; Chen, W.B.; Wang, R.S. Numerical experiment on the behaviour of oil spills in ice-infested waters in the bohai sea. Oceanol. Limnol. Sin.; 1999; 30, pp. 557-562.

15. Wilkinson, J.P.; Wadhams, P.; Hughes, N.E. Modelling the spread of oil under fast sea ice using three-dimensional multibeam sonar data. Geophys. Res. Lett.; 2007; 34, L22506. [DOI: https://dx.doi.org/10.1029/2007GL031754]

16. Xiao, M.; Gao, Q.J.; Lin, J.G.; Li, W.; Liang, X. Simulation of Submarine Pipeline Oil Spill Based on Wave Motion. Proceedings of the International Conference on Computer Modeling and Simulation; Hainan, China, 22–24 January 2010; Volume 1, pp. 433-437.

17. Liu, X.F.; Meng, L.K.; Zhao, C.Y. DEM-based Oil and Gas Migration Pathway Simulation of Oil-and Gas-Bearing Basin. Editor. Board Geomat. Inf. Sci. Wuhan Univ.; 2004; 29, pp. 371-375.

18. Ma, E.; Xu, Z.M. Technological process and optimum design of organic materials vacuum pyrolysis and indium chlorinated separation from waste liquid crystal display panels. J. Hazard. Mater.; 2013; 263, pp. 610-617. [DOI: https://dx.doi.org/10.1016/j.jhazmat.2013.10.020]

19. Lv, G.L.; Liu, Z.M. Simulation and Analysis of Diffusion of Oil and Gas Leakage from Submarine Pipeline Based on CFD. J. Jimei Univ.; 2017; 22, pp. 51-57.

20. Ren, B.; Li, X.L.; Wang, Y.X. An Irregular Wave Maker of Active Absorption with VOF Method. China Ocean Eng.; 2008; 4, pp. 94-105.

21. Li, Z.G.; Jiang, M.R.; Yu, J.X. Numerical simulation on the oil spill for the submarine pipeline based on VOF method. Ocean Eng.; 2016; 34, pp. 100-110.

22. Li, Z.L.; Liu, Y.; Sun, S.S.; Lu, Y.L.; Ji, S.Y. Analysis of ship maneuvering performances and ice loads on ship hull with discrete element model in broken-ice fields. Chin. J. Theor. Appl. Mech.; 2013; 45, pp. 868-877.

23. Mohr, O. Welche Umstande bedingen die Elastizitatsgrenze undden bruch eines materials. Z. Des. Ver. Dtsch. Ing.; 1900; 44, pp. 1524-1530.

24. Egberongbe, F.O.; Nwilo, P.C.; Badejo, O.T. Oil Spill Disaster Monitoring along Nigerian Coastline. Proceedings of the 5th FIG Regional Conference: Promoting Land Administration and Good Governance; Accra, Ghana, 8–11 March 2006.

25. Cheng, R.T.; Casulli, V.; Milford, S.N. Eulerian-Lagrangian Solution of the Convection-Dispersion Equation in Natural Coordinates. Water Resour. Res.; 1984; 20, pp. 944-952. [DOI: https://dx.doi.org/10.1029/WR020i007p00944]

26. Yang, Z.Y. B-Spline interpolation in Eular-Lagrangian method for convection-dispersion equation. J. Ocean Univ. Qingdao; 1994; 2, pp. 143-151.

27. Rasmussen, T.L.; Thomsen, E.; Kuijpers, A.; Stefan, W. Late warming and early cooling of the sea surface in the Nordic seas during MIS 5e (Eemian Interglacial). Quat. Sci. Rev.; 2003; 22, pp. 809-821. [DOI: https://dx.doi.org/10.1016/S0277-3791(02)00254-8]

28. DF Dickins Associates. Field Research Spills to Investigate the Physical and Chemical Fate of Oil in Pack Ice; Environmental Studies Revolving Funds: Canada, 1987; 118p.

29. Sun, H.; Lu, P.; Li, Z.J. Numerical Simulation on Dynamic Characteristics of Flow Field under Ice. Math. Pract. Theory; 2015; 45, pp. 69-77.

30. Tu, X.C. Ice Resistance Prediction and Parameters Sensitivity Study for Polar Geophysical Prospecting Vessel. Master Thesis; Jiangsu University of Science and Technology: Zhenjiang, China, 2019.

31. Ma, S.; Ge, W.P.; Duan, W.Y.; Liu, H.X. Simulation of free decay roll for C11 container ship based on overset gird. J. Huazhong Univ. Sci. Technol.; 2017; 45, pp. 34-39.

32. Li, Y.; Liu, B.X.; Lan, G.X.; Ma, L.; Gao, C. Study on Spectrum of Oil Film in Ice-Infested Waters. Guang Pu Xue Yu Guang Pu Fen Xi; 2010; 30, pp. 1018-1021.