## 3 Numerical Simulation of Oil Spill
### Numerical Model  
A continuous oil-spill accident, which happened in the geographical range of 66\({}^{\circ}\)24\({}^{\prime}\) N to 66\({}^{\circ}\)27\({}^{\prime}\) N and 15\({}^{\circ}\)42\({}^{\prime}\) W to 15\({}^{\circ}\)54\({}^{\prime}\) 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\({}^{\circ}\) N, 15.79\({}^{\circ}\) 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 [PERSON] and [PERSON] [28]. The computational domain (Figure 3) had a dimension of 2000 \(\times\) 1000 \(\times\) 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. \begin{table}
\begin{tabular}{l l} \hline \hline
**Properties** & **Parameters** \\ \hline Location & 66.43\({}^{\circ}\) N, 15.79\({}^{\circ}\) W \\ \hline Release situation & Accidental ship leakage \\ \hline Spill rate & 15 kg/s \\ \hline Spill duration & 200 s \\ \hline Oil density & 840 kg/m\({}^{3}\) \\ \hline Oil viscosity & 500 mPa-s \\ \hline seawater density & 1025 kg/m\({}^{3}\) \\ \hline ice density & 910 kg/m\({}^{3}\) \\ \hline Ice diameter & 5 m \\ \hline Ice concentration & Medium (30–80\%) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The parameters used in the oil-spill model. Figure 2: Computational domain of Icelandic waters. 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-\varepsilon\) model was selected to present the turbulence mode and the computational step was set to 1 s.  
Figure 4: Grids’ refinement around the oil spout. Figure 5: Initialization of free surface wave. Figure 3: Computational domain. 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.