# Mathematical Modeling of Marine Oil Spills in the Luanjiakou District, near the Port of Yantai

Authors: Daming Li; Xingchen Tang; Yanqing Li; Xiao Wang; Hongqiang Zhang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

This paper presents a simulation method for oil spills in a multi-island area. The simulation considers three parts, which consist of(1) the spreading of an oil slick on its edge as well as the diffusion and drift under dynamic actions, (2) the evaporation and spreading thickness of an oil slick in its interior, and (3) the adsorption and emulsification near shorelines and islands. The Euler-Lagrange method is adopted to track the spill location and particles positions on the edge of oil slicks. A mathematical model of marine oil spills is established for the Luanjiakou District of the Port of Yantai. The flow field verification shows that the BIAS of tidal level, flow velocity, and flow direction is below ±10 cm, 0.11 m/s, and ±2°, respectively, and the oil spill verification captures satisfactory results. Hence, the proposed model could reproduce the oil spill process in this region. Then, we simulate oil spills under various operating conditions. It is concluded that the transport of oil slicks is mainly influenced by flood/ebb currents, whereas the wind plays a major role in the drift and thickness of oil slicks. The study provides an important reference to controlling and handling of accidental oil spills.


---

## Body

## 1. Introduction

With the rapid development of modern industry, human demand for energy has grown. The increasing energy demands have resulted in the increased exploitation of fossil energy, such as crude oil and gas [1–3]. Consequently, the risk of oil spill accidents has also increased. Coupled with the booming development of maritime transportation, ship grounding accidents, collisions, and capsized tankers have also contributed to the frequent occurrence of oil spill pollution. These accidents can result in severe damage to the marine environment [4–6]. Therefore, in recent years, governments worldwide have strengthened the overall management of the water transportation of oil and have fostered scientific research on the movement properties of oil in water [1–3, 7].Spreading is the horizontal expansion of an oil slick due to gravity, inertia, viscosity, and surface tension forces, which plays an important role in the fate process of surface spilled oil. The transport and fate of spilled oil in water are processes affected by dynamic factors, nondynamic factors, and variable oil properties [8]. Dynamic factors include the gravity, inertia, viscosity, and surface tension forces, wind, tidal currents, and randomness of the motion process [9, 10]. Nondynamic factors mainly consist of weathering processes such as dissolution, evaporation, photooxidation, emulsification, biodegradation, and other chemical changes. Thus, to accurately simulate the transport process of an oil slick under the influence of these dynamic factors is a significant component of oil spill modeling [11]. At present, oil spill models chiefly consider the following four aspects: (1) the spreading and drift of spilled oil on the sea, (2) the evaporation and adsorption of spilled oil, (3) the dispersion of spilled oil in water, and (4) the trajectories and fate of spilled oil [3].Early researchers adopted various numerical methods to simulate the movement of an oil slick based on the advection-diffusion equation. In the mid-1990s, an oil-particles model was developed by Johansen, Elliot, and Reed [12]. Other commercial models, such as NOAA [13], OLIMAP [14], OSIS [15], and OSCAR [16], have been subsequently adopted to simulate oil movement and distribution in water. In addition, the bio-optical forecasting (BioCast) system can be used to forecast the ocean’s color optical environment and the nowcast-forecast studies of ocean color data streams in physical circulation models [17], which have promoted the further development of oil spill models. What is more, Radarsat SAR technology is an effective way for marine oil spill detection, which can identify and extract oil spill information on the sea surface quickly and accurately and provide remote sensing image to calibrate the modeling result and improve the forecast precision. It is of great significance for the protection of the marine environment and nearshore ecological environment [18, 19].In recent years, there have been breakthroughs in terms of understanding the complex geometry for oil spill model research. Gjosteen [20] developed an oil-spreading model based on the conservation of volume and momentum, which is suitable for coupling with the discrete-element ice model and other complex boundaries for oil spreading. Dietrich et al. [6] applied a coupled SWAN + ADCIRC model on a high-resolution computational mesh and the highly efficient Lagrangian particle transport model to simulate the surface trajectories of the oil spill from the Macondo well, northern Gulf of Mexico, during the spring of 2010. Wang and Shen [21, 22] developed a three-dimensional integrated model that provides great flexibility for modeling oil spill accidents in complex geometries such as tidal creeks, barriers, and islands. Currently, an unstructured grid is primarily used to simulate the solid boundary for the study of oil spills in complex terrains [6, 21–28]. Alves et al. [1] proposed a three-step model to predict and assess shoreline and offshore susceptibility to oil spills for confined basins with islands, which comprises the leading edge in the study of oil spills in complex terrains. Then, Alves et al. [2] simulated a series of oil spill accidents in the Eastern Mediterranean Sea, which included confined maritime basins. Nevertheless, further studies need to be done in the problem of oil particles penetrating the solid boundary such as a complex multi-island terrain or hydraulic structures.In this study, an oil spill simulation method in a multi-island area is presented, and the processing modes for oil slick longshore transport and its penetration-resistant boundaries are developed. In addition, a local search method that can specify the search radius is proposed and adopted. What is more, the Euler-Lagrange method is adopted to track the spill location and the position of particles on the edge of oil slicks, which can easily calculate the oil slick area. Based on the Monte Carlo method, a mathematical model for marine oil spills is established to simulate the movement of oil spills in the Luanjiakou District of the Port of Yantai (Figure1).Figure 1
Map of the Luanjiakou District near the Port of Yantai and survey stations.


## 2. Theories and Methods

The Euler-Lagrange systems are divided into two parts, that is, Euler and Lagrange approaches. The Euler approach describes the distributions of the variables in flow field at any time, whose computational mesh is fixed in space, while the Lagrange approach traces each particle from a certain time and describes its trajectory, whose computational mesh is fixed on the centroid of the research object and can be used to simulate the trajectory of an oil slick [29, 30]. Fundamentally, the accuracy of the velocity information for the hydrodynamic field established by the Eulerian approach has a crucial role in the successful prediction of the oil spill trajectory by adopting the Lagrange approach [25]. Thus, the coupling approach can give full play to the advantages of both approaches and avoid their defects, which is the reason why it is widely used in solving two-dimensional hydrodynamic problems by using finite element analysis method.


### 2.1. Hydrodynamic Module

In the coastal area, the horizontal movement scale of the tidal current is much larger than the vertical movement scale and the hydraulic parameters are unconspicuous in the vertical direction, so the flow field can be expressed by the average flow quantities along the direction of the water depth [31]. If the three-dimensional simulation is adopted, the calculation would be greatly increased and the simulation conditions would be more complex due to the impact of the vertical stratification on open boundary conditions, wind stress conditions, bottom friction forms, and so on. As said above, the governing equations of hydrodynamic model are the two-dimensional depth-averaged shallow water circulation equations, which are discretized by the finite element method. The continuity equation is given by(1)∂z∂t+∂hu∂x+∂hv∂y=0.The equations for conservation of momentum are given by(2)∂u∂t+u∂u∂x+v∂u∂y+g∂z∂x=-gu2+v2c2hu+fv+τxsρh+Ah∂2u∂x2+∂2u∂y2∂v∂t+u∂v∂x+v∂v∂y+g∂z∂y=-gu2+v2c2hv-fu+τysρh+Ah∂2v∂x2+∂2v∂y2,where u and v are the component velocities of current in the x and y directions, respectively, z is the water level, g is the gravitational acceleration, c is Chezy’s friction coefficient, h is the water depth, f is the Coriolis coefficient, Ah is the horizontal eddy viscosity coefficient, ρ is the density of water, and τxs and τys are the wind stress components on the sea surface in the x and y directions, respectively.Waves are generally induced by the wind on the sea surface. According to the temporal and spatial variation, the waves can be divided into regular and irregular waves. Regular waves have constant amplitudes and wavelengths, whose waveforms do not vary with time and space. It is possible to form such waves only when the problem is two-dimensional, the water depth is constant, and the disturbance source generating wave periodically varies with time. Irregular wave patterns are transient, whose elements vary with time and space. The waves generated when the wind blows over the sea surface are a common type of irregular wave [32]. Thus, the wind wave is just taken into account in the wave calculation, that is, wind-induced shearing stresses on water surfaces τxs and τys, which can be computed by the following empirical formula:(3)τxs=0.125CDWxW→τys=0.125CDWyW→,where Wx and Wy are the component velocities of the wind on the sea surface in the x and y directions, respectively, W→ is the wind velocity vector, and the wind drag coefficient CD can be obtained from the following Heaps empirical formula:(4)CD=0.564×10-3,W→≤4.917-0.12+0.13W→×10-3,4.917<W→<19.2212.513×10-3,W→≥19.221.


### 2.2. Oil Spill Module

#### 2.2.1. Spreading

Fay’s [33] three-phase spreading theory is adopted to study the spreading of the oil slick on the still water surface, which is based on laboratory hydrostatic experiments [34]. The spreading diameter of the oil slick in each phase can be expressed by Fay’s empirical formulas:Gravity-inertia spreading phase:(5)L1=K1gVΔ1/4t1/2Gravity-viscous spreading phase:(6)L2=K2g1-ΔΔ1/6υw-1/12t1/4V1/4Surface tension viscous spreading phase:(7)L3=K3g-1/4δ1/2ρw-3/4t3/4,where Δ=1-ρ0/ρw, where ρ0 and ρw are the density of oil particles and water, respectively, L is the spreading diameter of the oil slick, and the subscripts denote different phases, δ=δwa-δao-δow, where δwa, δao, and δow are the water-air, oil slick-air, and oil slick-water surface tension coefficients, respectively, υw is the kinematic viscosity of water, K1=1.35, K2=1.60, and K3=0.48 are experiment constants, t is the duration of the oil spill, V=∑i=1nQi(t)Δt[1-k(t-iΔt)] is the spill volume, Q is the spill discharge, k is the synthetic attenuation coefficient, and n=t/Δt is the number of oil spill times.


#### 2.2.2. Drift and Thickness

The drift of the oil slick is a vector sum of surface current and wind, of which the velocity vector is shown in Figure2. The drift velocity u⃑r can be written as(8)u⃑r=Kcu→c+Kwu→w,where u→c is the surface current velocity, u→w is the wind velocity at 10 m above the water surface, Kc=1 is the current factor, and Kw=0.035 is the wind drift factor [35].Figure 2
Velocity vector of oil slick drift.The spreading thicknessh can be determined from the mass conservation equation, as follows:(9)∂Ch∂t+∂uCh∂x+∂vCh∂y=-Φs+Φb+R,where C is the oil slick concentration, Φs and Φb are the spill flux of the upper and lower surface of the oil slick, respectively, and R is the loss of the oil slick in the chemical and physical processes.The spill flux is very similar to the mass transfer flux in molecular diffusion, so the spreading thicknessh can be computed from Fick’s law:(10)h=Ve-kt2πσsσnexp⁡-s22σs2-n22σn2,where s and n are the natural coordinates in the oil slick drift direction and the direction perpendicular to s, respectively, k is the attenuation coefficient of oil, and σs=ast1.17 and σn=1/10σn are the standard deviation of the oil slick thickness in the s and n directions, respectively.


#### 2.2.3. Particle Motion Model of Oil Slick

Due to the influence of various dynamic factors, such as wind, wave, and current, the diffusion of spilled oil on the sea surface has certain randomness at any time, which can be properly described by the Monte Carlo method [36]. It captures the multiple sampling data of the function based on the sampling of each of the random variables and then calculates the function value of each group from the independent sampling data, so as to determine the probability distribution of the function. It is applied to the problem of oil spill diffusion, which is to obtain the movement direction and displacement of oil particles by giving each of the tracked particles a set of random numbers under the premise of determining the disturbance intensity and time scale. Namely, the trajectories of oil particles are captured by adding a random term to the result obtained by the Lagrange method. The essence is to help supplement and revise the Euler-Lagrange systems.The Monte Carlo method is adopted to calculate the oil movement in the present study. First, the spill location and the position of particles on the edge of oil slicks are tracked and recorded by using the Euler-Lagrange method. Next, the diffusion random number is added to the module. As a result, the action of the wave-guide and wind-induced currents on dispersion and fragmentation of oil slicks is taken into account to describe the trajectory and irregular shape of these same spills.Assuming the sampling stepΔt>0 and Xn=X(nΔt), we have(11)Xn=Xn-1+σΔtWnσ>0,where Wn are independent random numbers on N0,1 and the increment Xn-Xn-k depends only on k variables Wn-k+1⋯Wnk<n corresponding to n-k,n, so Xn-Xn-k follows the normal distribution N0,σkΔt.Specifically, supposing that position coordinates of an oil particle areRti and Rti+1 at times ti and ti+1, respectively, and R is the movement distance of an oil particle under the actions of spreading, drift, and so on, we will then have(12)R→ti+1=R→ti+R→R→=γR→s+R→l,where γ is the random number ranging from 0 to 1 and R→l is the drift vector in the period of Δt, which can be obtained by integrating the Lagrange velocity as follows. The Lagrange velocity can be approximately represented by the Euler velocity in the calculation:(13)R→l=∫titi+1u→rdt.R→s is the spreading vector in the period of Δt, which is given by(14)R→s=l→ti+1-l→ti.The discrete form of transport distance of the labeled oil particle can be obtained by (15)R→ti+1=R→ti+Δt·u→r+γl→ti+1-l→ti.In addition, u→r can be obtained from (8), and l→ti+1 and l→ti can be obtained from (5)–(7). As mentioned above, we can calculate the transport position of each of the oil particles. A large number of oil particles can reflect the behavior processes of marine oil spills.


#### 2.2.4. Evaporation and Emulsification

The evaporation rate is influenced by the temperature, waves, wind speed, and oil slick areas, among other factors. Hence, the evaporation amount of surface oil slick can be calculated by the following distillation formula [37]:(16)Fv=ln⁡1+BTG/Tθexp⁡A-BT0/TT/BTG,where Fv is the volume fraction evaporated, A and B are the constants, usually selected as 6.3 and 10.3 for crude oils, respectively, TG is the slope of distillation curve, T is the surface temperature of the oil slick, T0=542.6-30.275API+1.565API2-0.03439API3+0.0002604API4 is the initial boiling point [38], and API is the density of spilt oil following the classification of the American Petroleum Institute; θ=0.0025(uw+1)0.78×2.437/h is the exposure coefficient of oil slick, and S(t) is the area of the oil slick.When drifting on the sea surface under the influence of wind and waves, oil particles disperse to the aqueous phase and water particles also disperse to the oil phase continuously. Subsequently, an oily emulsion is generated. The emulsification equation is given by [39, 40](17)dFwdt=C1uw+u021-C2FwQt,where Fw is the emulsification fraction, C1 is the absorption rate, usually selected as 2×10-6, C2 is the water content, usually adopted as 1.33, u0 is the emulsification correction factor in the ocean environment, and Q(t) is the emulsifying amount of the oil slick.


#### 2.2.5. Computation Mode

In this study, a simulation method for oil spills in a multi-island area is presented to simultaneously observe and study the edge and centroid motion of an oil slick (see Figure3). It is suggested that a number of discrete nodes are distributed along the edge of the oil slick, and there is a line along the edge of the oil slick between the nodes, which is called the “oil surface.” The number of nodes can be increased or decreased appropriately depending on the degree of density so that the edge interface can be expressed by a continuous and smooth edge line. The interface is referred to as “combination of particles and oil surface.” This way, the motion quantities of the discrete nodes can be calculated. Therefore, the model can entirely simulate the motion process of the oil slick, including the spreading of the oil slick on its edge, the diffusion and drift under the dynamic actions of wind, waves, and currents, the evaporation and thickness of the oil slick in its interior, and the adsorption and emulsification of the oil slick near shorelines and islands.Figure 3
Description of the computing mode of the oil slick (the area surrounded by the solid line represents the oil slick, black points represent discrete nodes along the edge of the oil slick, line between discrete nodes represents oil surface, and the area surrounded by the dashed line represents the combination of particles and oil surface).Marine oil spill models usually cover large areas using many grids. Furthermore, in most calculations one does not only need to determine the scope of the search unit, but also need to ascertain whether or not the search node is in this unit. In addition, the centroid and edge of an oil slick are not necessarily near the previous location because the transport of the oil slick with water movement may be very large over a short period. However, using the global search method (i.e., searching the entire study area) would lead to the huge calculation. Therefore, the local search method is proposed in this paper, which specifies the search radius, thereby reducing the amount of computation (see Figure4). As shown in the figure, the position of the node in the previous moment is taken as the circle center and the search radius is provided. In addition, the unit number is arbitrary and its centroid coordinate is provided. This way, we can determine whether or not the unit is within the search range.Figure 4
Schematic diagram of the local search method (red circular area for the search range, pink point for the circle center, yellow point for the search node, yellow arrow for the search radius, blue solid line for the contour line of an oil slick in the previous moment, and blue dashed line for the contour line of an oil slick in the present moment).During oil spills around multi-island areas, coastal structures such as breakwaters, quays, jetties, wharfs, and docks are likely obstacles to the spreading and transport of oil slicks [3, 10]. When transporting along these obstacles, a portion of the oil slick would be adsorbed in the structures. Note that the permanent absorption is taken into account in this study. Hence, the mode of the penetration-resistant boundary that is the case where oil particles are transported along the coast and adsorbed on it and do not penetrate the solid boundary, is developed (see Figure 5(a)). In addition, the mode can be used for real-time detection of the solid boundary. Then, the adsorption unit and location of oil particles can be ascertained using the unit information recorded by the local search method (see Figure 5(a)). This strategy is a good way to avoid the unlikely case of oil particles penetrating the solid boundary when the current velocity is relatively large (see Figure 5(b)).Figure 5
Comparison of different movement conditions of oil particles (black point) when arriving at the solid boundary (solid line) ((a) represents the modes of the penetration-resistant boundary as well as the longshore transport and adsorption of the oil slick and (b) represents the unlikely case of oil particles penetrating the solid boundary).
(a)
(b)


#### 2.2.6. Oil Spill Verification

(
1) Oil Spill on a Still Water Surface. The spreading and extension of an oil slick are some of the main differences between oil spill diffusion and concentration diffusion, which is reflected by the major and minor axes of the oil slick changing with time. Thus, the scales of the major and minor axes of the oil slick after an instantaneous oil spill are simulated under different oil volumes (Figure 6). A comparison of the numerical results with the results obtained by Zhao and Wu [41] shows good conformity in the major/minor axes scales (see Tables 1 and 2). Moreover, the numerical results of the two studies converge with increasing oil volume. There is slightly larger discrepancy between simulated and reference results for major axes as compared to minor axes. The reason for this is that the major axes of the oil slick are deeply influenced by many factors, such as wind, waves, and currents.Table 1
Comparison of the major axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
15.33
28.14
45.1
64.83
Simulated values of [41] (km)
12.67
25.91
43.54
65.49Table 2
Comparison of the minor axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
5.99
11.69
20.69
34.15
Simulated values of [41] (km)
5.18
11.74
21.10
34.04Figure 6
The major (a) and minor (b) axes of the oil slick versus time.
(a)
(b)(
2) Oil Spill on a Flowing Water Surface. The oil slick diffusion and drift experiment were carried out in a flume 25 cm long and 60 cm wide. The flow section for experimental observation is 1.17 m, in which the flow is uniform, and the mean flow velocity is approximately 0.04 m/s. The flume experiment and simulated results are shown in Figure 7, in which (a) and (c) are the oil slick diffusion and drift at different times in the flume and (b) is the simulated result. A comparison of the simulated and experimental results is shown in Table 3, which shows that the results are in good agreement with each other.Table 3
Comparison of the simulated and experimental results.
Item
Initial size (cm)
Final size (cm)
Movement distance (m)
Movement time (s)
Simulated results
15
21
1.17
30
Experimental results
15
22
1.2
30Figure 7
Comparison of the flume experiment (a, c) and the simulated result (b) of the spreading and drift of the oil slick.
(a)
(b)
(c)


## 2.1. Hydrodynamic Module

In the coastal area, the horizontal movement scale of the tidal current is much larger than the vertical movement scale and the hydraulic parameters are unconspicuous in the vertical direction, so the flow field can be expressed by the average flow quantities along the direction of the water depth [31]. If the three-dimensional simulation is adopted, the calculation would be greatly increased and the simulation conditions would be more complex due to the impact of the vertical stratification on open boundary conditions, wind stress conditions, bottom friction forms, and so on. As said above, the governing equations of hydrodynamic model are the two-dimensional depth-averaged shallow water circulation equations, which are discretized by the finite element method. The continuity equation is given by(1)∂z∂t+∂hu∂x+∂hv∂y=0.The equations for conservation of momentum are given by(2)∂u∂t+u∂u∂x+v∂u∂y+g∂z∂x=-gu2+v2c2hu+fv+τxsρh+Ah∂2u∂x2+∂2u∂y2∂v∂t+u∂v∂x+v∂v∂y+g∂z∂y=-gu2+v2c2hv-fu+τysρh+Ah∂2v∂x2+∂2v∂y2,where u and v are the component velocities of current in the x and y directions, respectively, z is the water level, g is the gravitational acceleration, c is Chezy’s friction coefficient, h is the water depth, f is the Coriolis coefficient, Ah is the horizontal eddy viscosity coefficient, ρ is the density of water, and τxs and τys are the wind stress components on the sea surface in the x and y directions, respectively.Waves are generally induced by the wind on the sea surface. According to the temporal and spatial variation, the waves can be divided into regular and irregular waves. Regular waves have constant amplitudes and wavelengths, whose waveforms do not vary with time and space. It is possible to form such waves only when the problem is two-dimensional, the water depth is constant, and the disturbance source generating wave periodically varies with time. Irregular wave patterns are transient, whose elements vary with time and space. The waves generated when the wind blows over the sea surface are a common type of irregular wave [32]. Thus, the wind wave is just taken into account in the wave calculation, that is, wind-induced shearing stresses on water surfaces τxs and τys, which can be computed by the following empirical formula:(3)τxs=0.125CDWxW→τys=0.125CDWyW→,where Wx and Wy are the component velocities of the wind on the sea surface in the x and y directions, respectively, W→ is the wind velocity vector, and the wind drag coefficient CD can be obtained from the following Heaps empirical formula:(4)CD=0.564×10-3,W→≤4.917-0.12+0.13W→×10-3,4.917<W→<19.2212.513×10-3,W→≥19.221.


## 2.2. Oil Spill Module

### 2.2.1. Spreading

Fay’s [33] three-phase spreading theory is adopted to study the spreading of the oil slick on the still water surface, which is based on laboratory hydrostatic experiments [34]. The spreading diameter of the oil slick in each phase can be expressed by Fay’s empirical formulas:Gravity-inertia spreading phase:(5)L1=K1gVΔ1/4t1/2Gravity-viscous spreading phase:(6)L2=K2g1-ΔΔ1/6υw-1/12t1/4V1/4Surface tension viscous spreading phase:(7)L3=K3g-1/4δ1/2ρw-3/4t3/4,where Δ=1-ρ0/ρw, where ρ0 and ρw are the density of oil particles and water, respectively, L is the spreading diameter of the oil slick, and the subscripts denote different phases, δ=δwa-δao-δow, where δwa, δao, and δow are the water-air, oil slick-air, and oil slick-water surface tension coefficients, respectively, υw is the kinematic viscosity of water, K1=1.35, K2=1.60, and K3=0.48 are experiment constants, t is the duration of the oil spill, V=∑i=1nQi(t)Δt[1-k(t-iΔt)] is the spill volume, Q is the spill discharge, k is the synthetic attenuation coefficient, and n=t/Δt is the number of oil spill times.


### 2.2.2. Drift and Thickness

The drift of the oil slick is a vector sum of surface current and wind, of which the velocity vector is shown in Figure2. The drift velocity u⃑r can be written as(8)u⃑r=Kcu→c+Kwu→w,where u→c is the surface current velocity, u→w is the wind velocity at 10 m above the water surface, Kc=1 is the current factor, and Kw=0.035 is the wind drift factor [35].Figure 2
Velocity vector of oil slick drift.The spreading thicknessh can be determined from the mass conservation equation, as follows:(9)∂Ch∂t+∂uCh∂x+∂vCh∂y=-Φs+Φb+R,where C is the oil slick concentration, Φs and Φb are the spill flux of the upper and lower surface of the oil slick, respectively, and R is the loss of the oil slick in the chemical and physical processes.The spill flux is very similar to the mass transfer flux in molecular diffusion, so the spreading thicknessh can be computed from Fick’s law:(10)h=Ve-kt2πσsσnexp⁡-s22σs2-n22σn2,where s and n are the natural coordinates in the oil slick drift direction and the direction perpendicular to s, respectively, k is the attenuation coefficient of oil, and σs=ast1.17 and σn=1/10σn are the standard deviation of the oil slick thickness in the s and n directions, respectively.


### 2.2.3. Particle Motion Model of Oil Slick

Due to the influence of various dynamic factors, such as wind, wave, and current, the diffusion of spilled oil on the sea surface has certain randomness at any time, which can be properly described by the Monte Carlo method [36]. It captures the multiple sampling data of the function based on the sampling of each of the random variables and then calculates the function value of each group from the independent sampling data, so as to determine the probability distribution of the function. It is applied to the problem of oil spill diffusion, which is to obtain the movement direction and displacement of oil particles by giving each of the tracked particles a set of random numbers under the premise of determining the disturbance intensity and time scale. Namely, the trajectories of oil particles are captured by adding a random term to the result obtained by the Lagrange method. The essence is to help supplement and revise the Euler-Lagrange systems.The Monte Carlo method is adopted to calculate the oil movement in the present study. First, the spill location and the position of particles on the edge of oil slicks are tracked and recorded by using the Euler-Lagrange method. Next, the diffusion random number is added to the module. As a result, the action of the wave-guide and wind-induced currents on dispersion and fragmentation of oil slicks is taken into account to describe the trajectory and irregular shape of these same spills.Assuming the sampling stepΔt>0 and Xn=X(nΔt), we have(11)Xn=Xn-1+σΔtWnσ>0,where Wn are independent random numbers on N0,1 and the increment Xn-Xn-k depends only on k variables Wn-k+1⋯Wnk<n corresponding to n-k,n, so Xn-Xn-k follows the normal distribution N0,σkΔt.Specifically, supposing that position coordinates of an oil particle areRti and Rti+1 at times ti and ti+1, respectively, and R is the movement distance of an oil particle under the actions of spreading, drift, and so on, we will then have(12)R→ti+1=R→ti+R→R→=γR→s+R→l,where γ is the random number ranging from 0 to 1 and R→l is the drift vector in the period of Δt, which can be obtained by integrating the Lagrange velocity as follows. The Lagrange velocity can be approximately represented by the Euler velocity in the calculation:(13)R→l=∫titi+1u→rdt.R→s is the spreading vector in the period of Δt, which is given by(14)R→s=l→ti+1-l→ti.The discrete form of transport distance of the labeled oil particle can be obtained by (15)R→ti+1=R→ti+Δt·u→r+γl→ti+1-l→ti.In addition, u→r can be obtained from (8), and l→ti+1 and l→ti can be obtained from (5)–(7). As mentioned above, we can calculate the transport position of each of the oil particles. A large number of oil particles can reflect the behavior processes of marine oil spills.


### 2.2.4. Evaporation and Emulsification

The evaporation rate is influenced by the temperature, waves, wind speed, and oil slick areas, among other factors. Hence, the evaporation amount of surface oil slick can be calculated by the following distillation formula [37]:(16)Fv=ln⁡1+BTG/Tθexp⁡A-BT0/TT/BTG,where Fv is the volume fraction evaporated, A and B are the constants, usually selected as 6.3 and 10.3 for crude oils, respectively, TG is the slope of distillation curve, T is the surface temperature of the oil slick, T0=542.6-30.275API+1.565API2-0.03439API3+0.0002604API4 is the initial boiling point [38], and API is the density of spilt oil following the classification of the American Petroleum Institute; θ=0.0025(uw+1)0.78×2.437/h is the exposure coefficient of oil slick, and S(t) is the area of the oil slick.When drifting on the sea surface under the influence of wind and waves, oil particles disperse to the aqueous phase and water particles also disperse to the oil phase continuously. Subsequently, an oily emulsion is generated. The emulsification equation is given by [39, 40](17)dFwdt=C1uw+u021-C2FwQt,where Fw is the emulsification fraction, C1 is the absorption rate, usually selected as 2×10-6, C2 is the water content, usually adopted as 1.33, u0 is the emulsification correction factor in the ocean environment, and Q(t) is the emulsifying amount of the oil slick.


### 2.2.5. Computation Mode

In this study, a simulation method for oil spills in a multi-island area is presented to simultaneously observe and study the edge and centroid motion of an oil slick (see Figure3). It is suggested that a number of discrete nodes are distributed along the edge of the oil slick, and there is a line along the edge of the oil slick between the nodes, which is called the “oil surface.” The number of nodes can be increased or decreased appropriately depending on the degree of density so that the edge interface can be expressed by a continuous and smooth edge line. The interface is referred to as “combination of particles and oil surface.” This way, the motion quantities of the discrete nodes can be calculated. Therefore, the model can entirely simulate the motion process of the oil slick, including the spreading of the oil slick on its edge, the diffusion and drift under the dynamic actions of wind, waves, and currents, the evaporation and thickness of the oil slick in its interior, and the adsorption and emulsification of the oil slick near shorelines and islands.Figure 3
Description of the computing mode of the oil slick (the area surrounded by the solid line represents the oil slick, black points represent discrete nodes along the edge of the oil slick, line between discrete nodes represents oil surface, and the area surrounded by the dashed line represents the combination of particles and oil surface).Marine oil spill models usually cover large areas using many grids. Furthermore, in most calculations one does not only need to determine the scope of the search unit, but also need to ascertain whether or not the search node is in this unit. In addition, the centroid and edge of an oil slick are not necessarily near the previous location because the transport of the oil slick with water movement may be very large over a short period. However, using the global search method (i.e., searching the entire study area) would lead to the huge calculation. Therefore, the local search method is proposed in this paper, which specifies the search radius, thereby reducing the amount of computation (see Figure4). As shown in the figure, the position of the node in the previous moment is taken as the circle center and the search radius is provided. In addition, the unit number is arbitrary and its centroid coordinate is provided. This way, we can determine whether or not the unit is within the search range.Figure 4
Schematic diagram of the local search method (red circular area for the search range, pink point for the circle center, yellow point for the search node, yellow arrow for the search radius, blue solid line for the contour line of an oil slick in the previous moment, and blue dashed line for the contour line of an oil slick in the present moment).During oil spills around multi-island areas, coastal structures such as breakwaters, quays, jetties, wharfs, and docks are likely obstacles to the spreading and transport of oil slicks [3, 10]. When transporting along these obstacles, a portion of the oil slick would be adsorbed in the structures. Note that the permanent absorption is taken into account in this study. Hence, the mode of the penetration-resistant boundary that is the case where oil particles are transported along the coast and adsorbed on it and do not penetrate the solid boundary, is developed (see Figure 5(a)). In addition, the mode can be used for real-time detection of the solid boundary. Then, the adsorption unit and location of oil particles can be ascertained using the unit information recorded by the local search method (see Figure 5(a)). This strategy is a good way to avoid the unlikely case of oil particles penetrating the solid boundary when the current velocity is relatively large (see Figure 5(b)).Figure 5
Comparison of different movement conditions of oil particles (black point) when arriving at the solid boundary (solid line) ((a) represents the modes of the penetration-resistant boundary as well as the longshore transport and adsorption of the oil slick and (b) represents the unlikely case of oil particles penetrating the solid boundary).
(a)
(b)


### 2.2.6. Oil Spill Verification

(
1) Oil Spill on a Still Water Surface. The spreading and extension of an oil slick are some of the main differences between oil spill diffusion and concentration diffusion, which is reflected by the major and minor axes of the oil slick changing with time. Thus, the scales of the major and minor axes of the oil slick after an instantaneous oil spill are simulated under different oil volumes (Figure 6). A comparison of the numerical results with the results obtained by Zhao and Wu [41] shows good conformity in the major/minor axes scales (see Tables 1 and 2). Moreover, the numerical results of the two studies converge with increasing oil volume. There is slightly larger discrepancy between simulated and reference results for major axes as compared to minor axes. The reason for this is that the major axes of the oil slick are deeply influenced by many factors, such as wind, waves, and currents.Table 1
Comparison of the major axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
15.33
28.14
45.1
64.83
Simulated values of [41] (km)
12.67
25.91
43.54
65.49Table 2
Comparison of the minor axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
5.99
11.69
20.69
34.15
Simulated values of [41] (km)
5.18
11.74
21.10
34.04Figure 6
The major (a) and minor (b) axes of the oil slick versus time.
(a)
(b)(
2) Oil Spill on a Flowing Water Surface. The oil slick diffusion and drift experiment were carried out in a flume 25 cm long and 60 cm wide. The flow section for experimental observation is 1.17 m, in which the flow is uniform, and the mean flow velocity is approximately 0.04 m/s. The flume experiment and simulated results are shown in Figure 7, in which (a) and (c) are the oil slick diffusion and drift at different times in the flume and (b) is the simulated result. A comparison of the simulated and experimental results is shown in Table 3, which shows that the results are in good agreement with each other.Table 3
Comparison of the simulated and experimental results.
Item
Initial size (cm)
Final size (cm)
Movement distance (m)
Movement time (s)
Simulated results
15
21
1.17
30
Experimental results
15
22
1.2
30Figure 7
Comparison of the flume experiment (a, c) and the simulated result (b) of the spreading and drift of the oil slick.
(a)
(b)
(c)


## 2.2.1. Spreading

Fay’s [33] three-phase spreading theory is adopted to study the spreading of the oil slick on the still water surface, which is based on laboratory hydrostatic experiments [34]. The spreading diameter of the oil slick in each phase can be expressed by Fay’s empirical formulas:Gravity-inertia spreading phase:(5)L1=K1gVΔ1/4t1/2Gravity-viscous spreading phase:(6)L2=K2g1-ΔΔ1/6υw-1/12t1/4V1/4Surface tension viscous spreading phase:(7)L3=K3g-1/4δ1/2ρw-3/4t3/4,where Δ=1-ρ0/ρw, where ρ0 and ρw are the density of oil particles and water, respectively, L is the spreading diameter of the oil slick, and the subscripts denote different phases, δ=δwa-δao-δow, where δwa, δao, and δow are the water-air, oil slick-air, and oil slick-water surface tension coefficients, respectively, υw is the kinematic viscosity of water, K1=1.35, K2=1.60, and K3=0.48 are experiment constants, t is the duration of the oil spill, V=∑i=1nQi(t)Δt[1-k(t-iΔt)] is the spill volume, Q is the spill discharge, k is the synthetic attenuation coefficient, and n=t/Δt is the number of oil spill times.


## 2.2.2. Drift and Thickness

The drift of the oil slick is a vector sum of surface current and wind, of which the velocity vector is shown in Figure2. The drift velocity u⃑r can be written as(8)u⃑r=Kcu→c+Kwu→w,where u→c is the surface current velocity, u→w is the wind velocity at 10 m above the water surface, Kc=1 is the current factor, and Kw=0.035 is the wind drift factor [35].Figure 2
Velocity vector of oil slick drift.The spreading thicknessh can be determined from the mass conservation equation, as follows:(9)∂Ch∂t+∂uCh∂x+∂vCh∂y=-Φs+Φb+R,where C is the oil slick concentration, Φs and Φb are the spill flux of the upper and lower surface of the oil slick, respectively, and R is the loss of the oil slick in the chemical and physical processes.The spill flux is very similar to the mass transfer flux in molecular diffusion, so the spreading thicknessh can be computed from Fick’s law:(10)h=Ve-kt2πσsσnexp⁡-s22σs2-n22σn2,where s and n are the natural coordinates in the oil slick drift direction and the direction perpendicular to s, respectively, k is the attenuation coefficient of oil, and σs=ast1.17 and σn=1/10σn are the standard deviation of the oil slick thickness in the s and n directions, respectively.


## 2.2.3. Particle Motion Model of Oil Slick

Due to the influence of various dynamic factors, such as wind, wave, and current, the diffusion of spilled oil on the sea surface has certain randomness at any time, which can be properly described by the Monte Carlo method [36]. It captures the multiple sampling data of the function based on the sampling of each of the random variables and then calculates the function value of each group from the independent sampling data, so as to determine the probability distribution of the function. It is applied to the problem of oil spill diffusion, which is to obtain the movement direction and displacement of oil particles by giving each of the tracked particles a set of random numbers under the premise of determining the disturbance intensity and time scale. Namely, the trajectories of oil particles are captured by adding a random term to the result obtained by the Lagrange method. The essence is to help supplement and revise the Euler-Lagrange systems.The Monte Carlo method is adopted to calculate the oil movement in the present study. First, the spill location and the position of particles on the edge of oil slicks are tracked and recorded by using the Euler-Lagrange method. Next, the diffusion random number is added to the module. As a result, the action of the wave-guide and wind-induced currents on dispersion and fragmentation of oil slicks is taken into account to describe the trajectory and irregular shape of these same spills.Assuming the sampling stepΔt>0 and Xn=X(nΔt), we have(11)Xn=Xn-1+σΔtWnσ>0,where Wn are independent random numbers on N0,1 and the increment Xn-Xn-k depends only on k variables Wn-k+1⋯Wnk<n corresponding to n-k,n, so Xn-Xn-k follows the normal distribution N0,σkΔt.Specifically, supposing that position coordinates of an oil particle areRti and Rti+1 at times ti and ti+1, respectively, and R is the movement distance of an oil particle under the actions of spreading, drift, and so on, we will then have(12)R→ti+1=R→ti+R→R→=γR→s+R→l,where γ is the random number ranging from 0 to 1 and R→l is the drift vector in the period of Δt, which can be obtained by integrating the Lagrange velocity as follows. The Lagrange velocity can be approximately represented by the Euler velocity in the calculation:(13)R→l=∫titi+1u→rdt.R→s is the spreading vector in the period of Δt, which is given by(14)R→s=l→ti+1-l→ti.The discrete form of transport distance of the labeled oil particle can be obtained by (15)R→ti+1=R→ti+Δt·u→r+γl→ti+1-l→ti.In addition, u→r can be obtained from (8), and l→ti+1 and l→ti can be obtained from (5)–(7). As mentioned above, we can calculate the transport position of each of the oil particles. A large number of oil particles can reflect the behavior processes of marine oil spills.


## 2.2.4. Evaporation and Emulsification

The evaporation rate is influenced by the temperature, waves, wind speed, and oil slick areas, among other factors. Hence, the evaporation amount of surface oil slick can be calculated by the following distillation formula [37]:(16)Fv=ln⁡1+BTG/Tθexp⁡A-BT0/TT/BTG,where Fv is the volume fraction evaporated, A and B are the constants, usually selected as 6.3 and 10.3 for crude oils, respectively, TG is the slope of distillation curve, T is the surface temperature of the oil slick, T0=542.6-30.275API+1.565API2-0.03439API3+0.0002604API4 is the initial boiling point [38], and API is the density of spilt oil following the classification of the American Petroleum Institute; θ=0.0025(uw+1)0.78×2.437/h is the exposure coefficient of oil slick, and S(t) is the area of the oil slick.When drifting on the sea surface under the influence of wind and waves, oil particles disperse to the aqueous phase and water particles also disperse to the oil phase continuously. Subsequently, an oily emulsion is generated. The emulsification equation is given by [39, 40](17)dFwdt=C1uw+u021-C2FwQt,where Fw is the emulsification fraction, C1 is the absorption rate, usually selected as 2×10-6, C2 is the water content, usually adopted as 1.33, u0 is the emulsification correction factor in the ocean environment, and Q(t) is the emulsifying amount of the oil slick.


## 2.2.5. Computation Mode

In this study, a simulation method for oil spills in a multi-island area is presented to simultaneously observe and study the edge and centroid motion of an oil slick (see Figure3). It is suggested that a number of discrete nodes are distributed along the edge of the oil slick, and there is a line along the edge of the oil slick between the nodes, which is called the “oil surface.” The number of nodes can be increased or decreased appropriately depending on the degree of density so that the edge interface can be expressed by a continuous and smooth edge line. The interface is referred to as “combination of particles and oil surface.” This way, the motion quantities of the discrete nodes can be calculated. Therefore, the model can entirely simulate the motion process of the oil slick, including the spreading of the oil slick on its edge, the diffusion and drift under the dynamic actions of wind, waves, and currents, the evaporation and thickness of the oil slick in its interior, and the adsorption and emulsification of the oil slick near shorelines and islands.Figure 3
Description of the computing mode of the oil slick (the area surrounded by the solid line represents the oil slick, black points represent discrete nodes along the edge of the oil slick, line between discrete nodes represents oil surface, and the area surrounded by the dashed line represents the combination of particles and oil surface).Marine oil spill models usually cover large areas using many grids. Furthermore, in most calculations one does not only need to determine the scope of the search unit, but also need to ascertain whether or not the search node is in this unit. In addition, the centroid and edge of an oil slick are not necessarily near the previous location because the transport of the oil slick with water movement may be very large over a short period. However, using the global search method (i.e., searching the entire study area) would lead to the huge calculation. Therefore, the local search method is proposed in this paper, which specifies the search radius, thereby reducing the amount of computation (see Figure4). As shown in the figure, the position of the node in the previous moment is taken as the circle center and the search radius is provided. In addition, the unit number is arbitrary and its centroid coordinate is provided. This way, we can determine whether or not the unit is within the search range.Figure 4
Schematic diagram of the local search method (red circular area for the search range, pink point for the circle center, yellow point for the search node, yellow arrow for the search radius, blue solid line for the contour line of an oil slick in the previous moment, and blue dashed line for the contour line of an oil slick in the present moment).During oil spills around multi-island areas, coastal structures such as breakwaters, quays, jetties, wharfs, and docks are likely obstacles to the spreading and transport of oil slicks [3, 10]. When transporting along these obstacles, a portion of the oil slick would be adsorbed in the structures. Note that the permanent absorption is taken into account in this study. Hence, the mode of the penetration-resistant boundary that is the case where oil particles are transported along the coast and adsorbed on it and do not penetrate the solid boundary, is developed (see Figure 5(a)). In addition, the mode can be used for real-time detection of the solid boundary. Then, the adsorption unit and location of oil particles can be ascertained using the unit information recorded by the local search method (see Figure 5(a)). This strategy is a good way to avoid the unlikely case of oil particles penetrating the solid boundary when the current velocity is relatively large (see Figure 5(b)).Figure 5
Comparison of different movement conditions of oil particles (black point) when arriving at the solid boundary (solid line) ((a) represents the modes of the penetration-resistant boundary as well as the longshore transport and adsorption of the oil slick and (b) represents the unlikely case of oil particles penetrating the solid boundary).
(a)
(b)


## 2.2.6. Oil Spill Verification

(
1) Oil Spill on a Still Water Surface. The spreading and extension of an oil slick are some of the main differences between oil spill diffusion and concentration diffusion, which is reflected by the major and minor axes of the oil slick changing with time. Thus, the scales of the major and minor axes of the oil slick after an instantaneous oil spill are simulated under different oil volumes (Figure 6). A comparison of the numerical results with the results obtained by Zhao and Wu [41] shows good conformity in the major/minor axes scales (see Tables 1 and 2). Moreover, the numerical results of the two studies converge with increasing oil volume. There is slightly larger discrepancy between simulated and reference results for major axes as compared to minor axes. The reason for this is that the major axes of the oil slick are deeply influenced by many factors, such as wind, waves, and currents.Table 1
Comparison of the major axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
15.33
28.14
45.1
64.83
Simulated values of [41] (km)
12.67
25.91
43.54
65.49Table 2
Comparison of the minor axes scales of the oil slick.
Spill volume (m3)
102
103
104
105
Simulated values of this paper (km)
5.99
11.69
20.69
34.15
Simulated values of [41] (km)
5.18
11.74
21.10
34.04Figure 6
The major (a) and minor (b) axes of the oil slick versus time.
(a)
(b)(
2) Oil Spill on a Flowing Water Surface. The oil slick diffusion and drift experiment were carried out in a flume 25 cm long and 60 cm wide. The flow section for experimental observation is 1.17 m, in which the flow is uniform, and the mean flow velocity is approximately 0.04 m/s. The flume experiment and simulated results are shown in Figure 7, in which (a) and (c) are the oil slick diffusion and drift at different times in the flume and (b) is the simulated result. A comparison of the simulated and experimental results is shown in Table 3, which shows that the results are in good agreement with each other.Table 3
Comparison of the simulated and experimental results.
Item
Initial size (cm)
Final size (cm)
Movement distance (m)
Movement time (s)
Simulated results
15
21
1.17
30
Experimental results
15
22
1.2
30Figure 7
Comparison of the flume experiment (a, c) and the simulated result (b) of the spreading and drift of the oil slick.
(a)
(b)
(c)


## 3. Model Setup and Verification

The Luanjiakou District is located in the western portion of Penglai-Yantai City, Shandong Peninsula. The district faces the Miaodao Islands, whose eastern coastline extends in the direction of Penglai City and the Yellow Sea and the western coastline extends in the direction of the Laizhou Gulf (see Figure1).


### 3.1. Study Area

The model domain and its bathymetry are shown in Figure8(a). The length of the domain is approximately 100 km, and its width is approximately 40 km extending to deep water, covering a sea area of approximately 4.6 × 104 km2. There are three open sea boundaries around, that is, the left, right, and upper straight boundaries. Triangular grids covering this domain were generated by the finite element method with a high grid resolution in the harbor, channel, and artificial island regions, with the following total number of grids and nodes: 47 and 244 and 24 and 350. The maximum grid spacing is approximately 2 km, and the minimum is approximately 0.025 km (see Figure 8(b)).Figure 8
(a) Bathymetry and (b) unstructured grids for the model domain.
(a)
(b)


### 3.2. Boundary Condition

To account for the lack of observational data, the astronomical tide, we induced the tidal level condition at the three open boundaries. Four main constituents in this domain are considered, that is, K1, M2, O1, and S2, whose harmonic constants can be derived from the global ocean tide model from the United States Department of the Navy [42], so that the tidal levels processes can be obtained at the open sea boundaries. Moreover, observational data are used for the landward boundaries.


### 3.3. Flow Field Verification

According to historical data [43], the survey stations are shown in Figure 1. The data from three survey stations (H1, H2, and H3), from 0:00 on July 4 to 18:00 on July 7, 2011, are adopted to validate tidal levels. The data from nine survey stations (U1, U2, U3, U4, U5, U6, U7, U8, and U9) of the diurnal tide, from 10:00 on July 5 to 14:00 on July 6, 2011, are used to validate flow velocity and directions.The validation results of the tidal level are shown in Figure9, which indicates that variations between the observed and the modeled results are in good agreement with each other. However, the tidal range is slightly different between the two. At high tide, the modeled values are smaller than the observed values, while at low tide the modeled values are larger than the observed values. This result could be related to datum selection prior to the modeling.Figure 9
Comparison of the tidal level between the modeled (solid line) and the observed (dots) results at three stations (H1, H2, and H3).There are many diurnal tide survey stations (see Figure1). Stations U1, U4, and U7 are used to illustrate our verifications of the flow velocity and direction (see Figures 10 and 11). In Figures 10 and 11, the variations of the flow velocity and direction between the observed and the modeled results are consistent at the three stations considered (U1, U4, and U7) except that there are deviations at individual times. The reason for this discrepancy may be associated with the accuracy of the observed data.Figure 10
Comparison of flow velocity between the modeled (solid line) and the observed (dots) results at three stations (U1, U4, and U7).Figure 11
Comparison of flow direction between the modeled (solid line) and the observed (dots) results at three stations (U1, U4, and U7).In particular, three criteria are adopted to assess the model performance for tidal level, flow velocity, and flow direction simulation, including the mean absolute error (MAE), the root mean square error (RMSE), and bias (BIAS) [19]. The equations for these three criteria are shown as follows:(18)MAE=1N∑i=1Nηim-ηioRMSE=1N∑i=1Nηim-ηio2BIAS=1N∑i=1nηim-ηio,where ηim are the modeled results and ηio are the observed results. The statistical errors for the differences between the simulated and observed results can be found in Table 4, from which it can be seen that, for the tidal level, the maximum RSME is 12.10 cm at Station H3 and the BIAS is below ±10 cm at three stations (H1, H2, and H3), for the flow velocity, the maximum RSME is 0.11 m/s at Station U1 and the BIAS is below ±0.10 m/s at three stations (U1, U4, and U7), and, for the flow direction, the maximum RSME is 17.63° at Station U1 and the BIAS is below ±2° at three stations (U1, U4, and U7).Table 4
Statistical errors at tidal survey stations for model verification.
Station
Tidal level
Station
Flow velocity
Flow direction
MAE (cm)
RSME (cm)
BIAS (cm)
MAE (m/s)
RSME (m/s)
BIAS (m/s)
MAE (deg)
RSME (deg)
BIAS (deg)
H1
9.18
11.04
-
8.11
U1
0.09
0.11
0.06
12.83
17.63
1.63
H2
8.29
10.32
-
6.83
U4
0.06
0.08
-
0.02
10.55
14.98
-
1.98
H3
10.02
12.10
-
9.13
U7
0.07
0.09
-
0.03
11.72
15.18
1.06The distributions of the flow field at ebb and flood periods are shown in Figure12. The results indicate that, during the ebb period, the velocities along the shoreline are much larger than those near the islands, because the water converges into the deep areas. During the flood period, velocity differences between the shoreline and the islands are less obvious. At both times, the tendencies of the flow field were well reflected by the model.Figure 12
Distributions of the flow field at the times of ebb (a) and flood (b).
(a)
(b)In summary, the hydrodynamic field can serve as the basis for studying marine oil spills in our study area.


### 3.4. Concentration Diffusion Verification

In the concentration diffusion verification of an oil slick, the results of a dyestuff tracing experiment carried out by South China Sea Institute of Oceanology, Academia Sinica, from 2:30 to 5:30 on January 29, 2002, were compared with the modeled results, as shown in Figure13. The figure shows that the diffusion tendency and range of the oil slick are relatively consistent, which provides the basis for the selection of the diffusion coefficient. It is indicated that the model can be adopted to reflect the actual oil slick movement in the region.Figure 13
Comparison between the experimental result (a) and the modeled result (b) of the concentration diffusion of the oil slick.
(a)
(b)


## 3.1. Study Area

The model domain and its bathymetry are shown in Figure8(a). The length of the domain is approximately 100 km, and its width is approximately 40 km extending to deep water, covering a sea area of approximately 4.6 × 104 km2. There are three open sea boundaries around, that is, the left, right, and upper straight boundaries. Triangular grids covering this domain were generated by the finite element method with a high grid resolution in the harbor, channel, and artificial island regions, with the following total number of grids and nodes: 47 and 244 and 24 and 350. The maximum grid spacing is approximately 2 km, and the minimum is approximately 0.025 km (see Figure 8(b)).Figure 8
(a) Bathymetry and (b) unstructured grids for the model domain.
(a)
(b)


## 3.2. Boundary Condition

To account for the lack of observational data, the astronomical tide, we induced the tidal level condition at the three open boundaries. Four main constituents in this domain are considered, that is, K1, M2, O1, and S2, whose harmonic constants can be derived from the global ocean tide model from the United States Department of the Navy [42], so that the tidal levels processes can be obtained at the open sea boundaries. Moreover, observational data are used for the landward boundaries.


## 3.3. Flow Field Verification

According to historical data [43], the survey stations are shown in Figure 1. The data from three survey stations (H1, H2, and H3), from 0:00 on July 4 to 18:00 on July 7, 2011, are adopted to validate tidal levels. The data from nine survey stations (U1, U2, U3, U4, U5, U6, U7, U8, and U9) of the diurnal tide, from 10:00 on July 5 to 14:00 on July 6, 2011, are used to validate flow velocity and directions.The validation results of the tidal level are shown in Figure9, which indicates that variations between the observed and the modeled results are in good agreement with each other. However, the tidal range is slightly different between the two. At high tide, the modeled values are smaller than the observed values, while at low tide the modeled values are larger than the observed values. This result could be related to datum selection prior to the modeling.Figure 9
Comparison of the tidal level between the modeled (solid line) and the observed (dots) results at three stations (H1, H2, and H3).There are many diurnal tide survey stations (see Figure1). Stations U1, U4, and U7 are used to illustrate our verifications of the flow velocity and direction (see Figures 10 and 11). In Figures 10 and 11, the variations of the flow velocity and direction between the observed and the modeled results are consistent at the three stations considered (U1, U4, and U7) except that there are deviations at individual times. The reason for this discrepancy may be associated with the accuracy of the observed data.Figure 10
Comparison of flow velocity between the modeled (solid line) and the observed (dots) results at three stations (U1, U4, and U7).Figure 11
Comparison of flow direction between the modeled (solid line) and the observed (dots) results at three stations (U1, U4, and U7).In particular, three criteria are adopted to assess the model performance for tidal level, flow velocity, and flow direction simulation, including the mean absolute error (MAE), the root mean square error (RMSE), and bias (BIAS) [19]. The equations for these three criteria are shown as follows:(18)MAE=1N∑i=1Nηim-ηioRMSE=1N∑i=1Nηim-ηio2BIAS=1N∑i=1nηim-ηio,where ηim are the modeled results and ηio are the observed results. The statistical errors for the differences between the simulated and observed results can be found in Table 4, from which it can be seen that, for the tidal level, the maximum RSME is 12.10 cm at Station H3 and the BIAS is below ±10 cm at three stations (H1, H2, and H3), for the flow velocity, the maximum RSME is 0.11 m/s at Station U1 and the BIAS is below ±0.10 m/s at three stations (U1, U4, and U7), and, for the flow direction, the maximum RSME is 17.63° at Station U1 and the BIAS is below ±2° at three stations (U1, U4, and U7).Table 4
Statistical errors at tidal survey stations for model verification.
Station
Tidal level
Station
Flow velocity
Flow direction
MAE (cm)
RSME (cm)
BIAS (cm)
MAE (m/s)
RSME (m/s)
BIAS (m/s)
MAE (deg)
RSME (deg)
BIAS (deg)
H1
9.18
11.04
-
8.11
U1
0.09
0.11
0.06
12.83
17.63
1.63
H2
8.29
10.32
-
6.83
U4
0.06
0.08
-
0.02
10.55
14.98
-
1.98
H3
10.02
12.10
-
9.13
U7
0.07
0.09
-
0.03
11.72
15.18
1.06The distributions of the flow field at ebb and flood periods are shown in Figure12. The results indicate that, during the ebb period, the velocities along the shoreline are much larger than those near the islands, because the water converges into the deep areas. During the flood period, velocity differences between the shoreline and the islands are less obvious. At both times, the tendencies of the flow field were well reflected by the model.Figure 12
Distributions of the flow field at the times of ebb (a) and flood (b).
(a)
(b)In summary, the hydrodynamic field can serve as the basis for studying marine oil spills in our study area.


## 3.4. Concentration Diffusion Verification

In the concentration diffusion verification of an oil slick, the results of a dyestuff tracing experiment carried out by South China Sea Institute of Oceanology, Academia Sinica, from 2:30 to 5:30 on January 29, 2002, were compared with the modeled results, as shown in Figure13. The figure shows that the diffusion tendency and range of the oil slick are relatively consistent, which provides the basis for the selection of the diffusion coefficient. It is indicated that the model can be adopted to reflect the actual oil slick movement in the region.Figure 13
Comparison between the experimental result (a) and the modeled result (b) of the concentration diffusion of the oil slick.
(a)
(b)


## 4. Results and Discussion

The port has 10,000-tonne tanker berths and the channel is an important shipping route for oil tankers and ships. Hence, the simulation assumes that spill locations are evenly distributed in the western, middle, and eastern portions of the port covering the entire channel, which are all the high-risk oil spill areas.According to the relevant specifications [44], the scenario simulations of marine oil spills are assumed and carried out in two ways: instantaneous and continuous. The condensate oil is used for the instantaneous oil spill scenario, and the spill volume is approximately 8000 t. For convenience of comparison, the low sulfur fuel oil is utilized for the continuous oil spill scenario, whose spill volume is constant and the duration is 10 h. The properties of the spilt oil are shown in Table 5.Table 5
Properties of the oil.
Name
Density (kg/m3)
Water content of emulsion (%)
API
Condensate oil
830.5
74
38.874
Low sulfur fuel oil
972
80
14.08In this region, the prevalent wind directions are SSW and S and the frequency is 15.14%. The static wind frequency is 0.47%. The strong wind directions are N, NW, and NNE, and the instantaneous maximum wind speed is 28 m/s [43]. The wind rose diagram for Luanjiakou District in 2002–2006 is shown in Figure 14. Together with live telecast data, the wind conditions in the model were set as shown in Table 6, in which Wind Direction 1 predominates in the sea area and the islands near the Miaodao Strait, Wind Direction 2 blows against the shoreline around the artificial islands, and Wind Direction 3 is unfavorable to the dock and harbor. The simulation time step was 60 s, and the time length was 48 h. To control the time, the initial minimum distinguishable spacing was 15 m, and the maximum distinguishable spacing was set as 100 m.Table 6
Wind conditions of the model.
Wind direction
No wind
Southwest wind (SW)
South wind (S)
Northwest wind (NW)
Northeast wind (NE)
Wind speed (m/s)
0
4.9
2.0
3.4
2.7
Note
Maximum wind direction
Wind Direction 1
Wind Direction 2
Wind Direction 3Figure 14
Wind rose diagram for Luanjiakou District in 2002–2006.


### 4.1. Spill Trajectories

The trajectories of instantaneous oil spills from the western portion of the channel under five wind conditions are shown in Figure15. In the figure, it can be seen that, in the case of no wind (Figure 15(a)), the oil slick migrated with flood/ebb currents and the area trajectory radiated towards the surrounding areas from the spill location because the ebb and flood velocities were roughly the same. When the oil spread to the narrow waterway of the Miaodao Strait, the ebb velocity increased and an oil slick zone protruding into the open sea appeared. Under the influence of southwest winds (Figure 15(b)), the oil slick after spill migrated towards the ebb, because the breakwater had little effect on the migration of the oil slick along the wind and flood/ebb directions. When removing the preventive area of the breakwater, the oil slick quickly spread to the Miaodao Islands, and the scope swept by the area trajectories was relatively large. Under the influence of south winds (Figure 15(c)), the oil slick approached the breakwater and then migrated towards the ebb due to the resistance of the breakwater. When removing the preventive area of the breakwater, the oil slick insufficiently spread, so the scope swept by the area trajectories was relatively small. Under the influence of northwest winds (Figure 15(d)), most of the oil slick after spill entered the Luanjiakou Port because the tidal current velocity was relatively small. Under the influence of northeast winds (Figure 15(e)), after drifting some distance with the ebb current, the oil slick moved to the southwest through passenger ferry berths and the port due to the combined effect of the wind and the flood current. Finally, part of the oil slick reached the western shoreline.Figure 15
Trajectories of instantaneous oil spills (red line) from the western portion of the channel (black star symbol for the western spill location) under five wind conditions ((a) represents oil spill trajectory in the case of no wind, (b) represents oil spill trajectory under the influence of southwest winds, (c) represents oil spill trajectory under the influence of south winds, (d) represents oil spill trajectory under the influence of northwest winds, and (e) represents oil spill trajectory under the influence of northeast winds).
(a)
(b)
(c)
(d)
(e)


### 4.2. Movement Process of Oil Slicks

Figures16 and 17 show the transport processes of instantaneous oil spills that occurred in the western portion of the channel in the case of no wind and the eastern portion of the channel under the influence of south winds, respectively. The figures show that oil slicks after spill migrated with the tidal current and wind, and they spread by themselves.Figure 16
Transport processes of instantaneous oil spills (red area) from the western portion of the channel (black star symbol for the western spill location) in the case of no wind.Figure 17
Transport processes of instantaneous oil spills (red area) from the eastern portion of the channel (red star symbol for the eastern spill location) under the influence of south winds.Figures18 and 19 show the transport processes of continuous oil spills that appeared in the western portion of the channel in the case of no wind and the eastern portion of the channel under the influence of south winds, respectively. The figures indicate that oil slicks after spill mixed with each other and that a narrow oil slick was formed. Then, oil slicks migrated with tidal current and wind, and they spread by themselves.Figure 18
Transport processes of continuous oil spills (red area) from the western portion of the channel (black star symbol for the western spill location) in the case of no wind.Figure 19
Transport processes of continuous oil spills (red area) from the eastern portion of the channel (red star symbol for the eastern spill location) under the influence of south winds.From Section2.2.2, it can be seen that the transport velocity of oil slicks is related to the local current velocity and the wind speed and that the spreading velocity is influenced by the spill volume, the density of the oil, and the surrounding terrain. Therefore, the instantaneously spilled oil drifted in the shape of the approximate ellipse. After bursting, an irregular multilayer ring was formed (see Figures 16 and 17). Conversely, the continuously spilled oil drifted in the shape of a narrow strip, and an irregular single-layer ring was finally formed (see Figures 18 and 19).


### 4.3. Area of Oil Slicks versus Time

Figures20–24 show the relationship of the slick area of instantaneous and continuous oil spills versus time. The results show that, in the case of no wind (Figure 20), the spreading area of instantaneous and continuous oil spills reached the maximums within 48 h. Under the influence of southwest winds (Figure 21), the maximum spreading area appeared in the eastern spill location. Under the influence of south winds (Figure 22), the maximum spreading area appeared in the middle spill location. Under the influence of northwest winds (Figure 23), the maximum spreading area of an instantaneous oil spill appeared in the western spill location, and the maximum spreading area of a continuous oil spill appeared in the middle spill location. Under the influence of northeast winds (Figure 24), the maximum spreading area of the instantaneous oil spill appeared in the western spill location and the maximum spreading area of the continuous oil spill appeared in the eastern spill location.Figure 20
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) in the case of no wind.
(a)
(b)Figure 21
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of southwest winds.
(a)
(b)Figure 22
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of south winds.
(a)
(b)Figure 23
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northwest winds.
(a)
(b)Figure 24
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northeast winds.
(a)
(b)From Figures20–24, it can be concluded that the maximum spreading area of oil slicks occurred in the eastern location, which spilled quickly under the influence of southwest winds and reached 109.385 km2 after 48 h. The minimum area occurred in the western and middle locations and reached 0.823 km2, which was continuously spilling under the influence of northwest and northeast winds, respectively.


### 4.4. Thickness of Oil Slicks versus Time

Figures25 and 26 show the relationship of the slick thickness of instantaneous and continuous oil spills versus time under different conditions. It can be observed that the thickness of oil slicks was relatively large in the beginning and gradually decreased with spreading and drift. When obstructed by the shoreline, oil slicks accumulated and the thickness suddenly increased or remained constant. After spilling for 48 h, the maximum thickness of oil slicks was approximately 9.998 mm, which mainly occurred under the influence of northwest and northeast winds. Due to the small current velocity near the shoreline, harbors, and islands, the wind squeezed oil slicks and limited the spreading and drift of them, forming a thicker oil slick area in the vicinity.Figure 25
Slick thickness of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) in the case of no wind.
(a)
(b)Figure 26
Slick thickness of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northeast winds.
(a)
(b)


### 4.5. Fate Process of Oil Volume

In the present study, the oil fate mainly includes the oil on the sea surface, evaporated, emulsified, and adsorbed near the shoreline after coming ashore. Figure27 shows the fate processes of the instantaneous oil spills, where the following can be observed: the initial oil volume on the sea surface is relatively large and then decreased slowly after the 48 hours due to evaporation, emulsification, and adsorption; evaporated and emulsified oil volume relate to the wind speed on the sea surface, whose tendencies are gradually increasing and then tend to be stable; the oil slick would be adsorbed when coming ashore, so the corresponding oil volume is also increasing.Figure 27
Fate processes of the instantaneous oil spill that occurred in the west of the channel in the case without wind (a) and in the east of the channel under the action of northwest wind (b).
(a)
(b)Figure28 shows the fate processes of the continuous oil spills, where it can be observed that the oil volume on the sea surface gradually increases during the initial 10 h and then the tendency is basically consistent with the instantaneous oil spill. And the other fate processes are in agreement with the instantaneous oil spill.Figure 28
Fate processes of the continuous oil spill that occurred in the west of the channel in the case without wind (a) and in the east of the channel under the action of northwest wind (b).
(a)
(b)


### 4.6. Future Work

The scenario simulations of marine oil spills in this study were preliminary using a two-dimensional oil spill model, which is actually a large-scale simulation in large areas. Further work remains to be done to improve the model performance, such as the multiscale simulation. For instance, the vertical diffusion of spilled oil in the water column can be carried out by the advanced SPH (Smoothed Particle Hydrodynamics) method, that is, the mesh-free particle method, which describes the transport of an oil slick with a series of particles and is more in coincidence with the idea of “oil-particles” model. In addition, the acquisition and usage of remote sensing information are essential to simulate and predict the marine oil spills in the near future due to its wide area coverage and the all-weather and all-day capabilities.


## 4.1. Spill Trajectories

The trajectories of instantaneous oil spills from the western portion of the channel under five wind conditions are shown in Figure15. In the figure, it can be seen that, in the case of no wind (Figure 15(a)), the oil slick migrated with flood/ebb currents and the area trajectory radiated towards the surrounding areas from the spill location because the ebb and flood velocities were roughly the same. When the oil spread to the narrow waterway of the Miaodao Strait, the ebb velocity increased and an oil slick zone protruding into the open sea appeared. Under the influence of southwest winds (Figure 15(b)), the oil slick after spill migrated towards the ebb, because the breakwater had little effect on the migration of the oil slick along the wind and flood/ebb directions. When removing the preventive area of the breakwater, the oil slick quickly spread to the Miaodao Islands, and the scope swept by the area trajectories was relatively large. Under the influence of south winds (Figure 15(c)), the oil slick approached the breakwater and then migrated towards the ebb due to the resistance of the breakwater. When removing the preventive area of the breakwater, the oil slick insufficiently spread, so the scope swept by the area trajectories was relatively small. Under the influence of northwest winds (Figure 15(d)), most of the oil slick after spill entered the Luanjiakou Port because the tidal current velocity was relatively small. Under the influence of northeast winds (Figure 15(e)), after drifting some distance with the ebb current, the oil slick moved to the southwest through passenger ferry berths and the port due to the combined effect of the wind and the flood current. Finally, part of the oil slick reached the western shoreline.Figure 15
Trajectories of instantaneous oil spills (red line) from the western portion of the channel (black star symbol for the western spill location) under five wind conditions ((a) represents oil spill trajectory in the case of no wind, (b) represents oil spill trajectory under the influence of southwest winds, (c) represents oil spill trajectory under the influence of south winds, (d) represents oil spill trajectory under the influence of northwest winds, and (e) represents oil spill trajectory under the influence of northeast winds).
(a)
(b)
(c)
(d)
(e)


## 4.2. Movement Process of Oil Slicks

Figures16 and 17 show the transport processes of instantaneous oil spills that occurred in the western portion of the channel in the case of no wind and the eastern portion of the channel under the influence of south winds, respectively. The figures show that oil slicks after spill migrated with the tidal current and wind, and they spread by themselves.Figure 16
Transport processes of instantaneous oil spills (red area) from the western portion of the channel (black star symbol for the western spill location) in the case of no wind.Figure 17
Transport processes of instantaneous oil spills (red area) from the eastern portion of the channel (red star symbol for the eastern spill location) under the influence of south winds.Figures18 and 19 show the transport processes of continuous oil spills that appeared in the western portion of the channel in the case of no wind and the eastern portion of the channel under the influence of south winds, respectively. The figures indicate that oil slicks after spill mixed with each other and that a narrow oil slick was formed. Then, oil slicks migrated with tidal current and wind, and they spread by themselves.Figure 18
Transport processes of continuous oil spills (red area) from the western portion of the channel (black star symbol for the western spill location) in the case of no wind.Figure 19
Transport processes of continuous oil spills (red area) from the eastern portion of the channel (red star symbol for the eastern spill location) under the influence of south winds.From Section2.2.2, it can be seen that the transport velocity of oil slicks is related to the local current velocity and the wind speed and that the spreading velocity is influenced by the spill volume, the density of the oil, and the surrounding terrain. Therefore, the instantaneously spilled oil drifted in the shape of the approximate ellipse. After bursting, an irregular multilayer ring was formed (see Figures 16 and 17). Conversely, the continuously spilled oil drifted in the shape of a narrow strip, and an irregular single-layer ring was finally formed (see Figures 18 and 19).


## 4.3. Area of Oil Slicks versus Time

Figures20–24 show the relationship of the slick area of instantaneous and continuous oil spills versus time. The results show that, in the case of no wind (Figure 20), the spreading area of instantaneous and continuous oil spills reached the maximums within 48 h. Under the influence of southwest winds (Figure 21), the maximum spreading area appeared in the eastern spill location. Under the influence of south winds (Figure 22), the maximum spreading area appeared in the middle spill location. Under the influence of northwest winds (Figure 23), the maximum spreading area of an instantaneous oil spill appeared in the western spill location, and the maximum spreading area of a continuous oil spill appeared in the middle spill location. Under the influence of northeast winds (Figure 24), the maximum spreading area of the instantaneous oil spill appeared in the western spill location and the maximum spreading area of the continuous oil spill appeared in the eastern spill location.Figure 20
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) in the case of no wind.
(a)
(b)Figure 21
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of southwest winds.
(a)
(b)Figure 22
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of south winds.
(a)
(b)Figure 23
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northwest winds.
(a)
(b)Figure 24
Area of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northeast winds.
(a)
(b)From Figures20–24, it can be concluded that the maximum spreading area of oil slicks occurred in the eastern location, which spilled quickly under the influence of southwest winds and reached 109.385 km2 after 48 h. The minimum area occurred in the western and middle locations and reached 0.823 km2, which was continuously spilling under the influence of northwest and northeast winds, respectively.


## 4.4. Thickness of Oil Slicks versus Time

Figures25 and 26 show the relationship of the slick thickness of instantaneous and continuous oil spills versus time under different conditions. It can be observed that the thickness of oil slicks was relatively large in the beginning and gradually decreased with spreading and drift. When obstructed by the shoreline, oil slicks accumulated and the thickness suddenly increased or remained constant. After spilling for 48 h, the maximum thickness of oil slicks was approximately 9.998 mm, which mainly occurred under the influence of northwest and northeast winds. Due to the small current velocity near the shoreline, harbors, and islands, the wind squeezed oil slicks and limited the spreading and drift of them, forming a thicker oil slick area in the vicinity.Figure 25
Slick thickness of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) in the case of no wind.
(a)
(b)Figure 26
Slick thickness of instantaneous (a) and continuous (b) oil spills versus time in different spill locations (black line for the western spill location, blue line for the middle spill location, and red line for the eastern spill location) under the influence of northeast winds.
(a)
(b)


## 4.5. Fate Process of Oil Volume

In the present study, the oil fate mainly includes the oil on the sea surface, evaporated, emulsified, and adsorbed near the shoreline after coming ashore. Figure27 shows the fate processes of the instantaneous oil spills, where the following can be observed: the initial oil volume on the sea surface is relatively large and then decreased slowly after the 48 hours due to evaporation, emulsification, and adsorption; evaporated and emulsified oil volume relate to the wind speed on the sea surface, whose tendencies are gradually increasing and then tend to be stable; the oil slick would be adsorbed when coming ashore, so the corresponding oil volume is also increasing.Figure 27
Fate processes of the instantaneous oil spill that occurred in the west of the channel in the case without wind (a) and in the east of the channel under the action of northwest wind (b).
(a)
(b)Figure28 shows the fate processes of the continuous oil spills, where it can be observed that the oil volume on the sea surface gradually increases during the initial 10 h and then the tendency is basically consistent with the instantaneous oil spill. And the other fate processes are in agreement with the instantaneous oil spill.Figure 28
Fate processes of the continuous oil spill that occurred in the west of the channel in the case without wind (a) and in the east of the channel under the action of northwest wind (b).
(a)
(b)


## 4.6. Future Work

The scenario simulations of marine oil spills in this study were preliminary using a two-dimensional oil spill model, which is actually a large-scale simulation in large areas. Further work remains to be done to improve the model performance, such as the multiscale simulation. For instance, the vertical diffusion of spilled oil in the water column can be carried out by the advanced SPH (Smoothed Particle Hydrodynamics) method, that is, the mesh-free particle method, which describes the transport of an oil slick with a series of particles and is more in coincidence with the idea of “oil-particles” model. In addition, the acquisition and usage of remote sensing information are essential to simulate and predict the marine oil spills in the near future due to its wide area coverage and the all-weather and all-day capabilities.


## 5. Conclusions

In this paper, a simulation method for the spreading and drift of an oil slick in a multi-island area and the mode of the penetration-resistant solid boundary are presented. To improve the computation efficiency, a local search method that can specify the search radius is adopted. The Euler-Lagrange method is adopted to track the spill location and the position of particles on the edge of oil slicks in order to calculate the slick area easily. Based on the Monte Carlo method, a mathematical model for marine oil spills was established for the Luanjiakou District, near the Port of Yantai. A series of verifications of the tidal current field and the movement of an oil slick show that the model can reflect the actual oil slick movement.The model has been applied to simulate the movement of oil slicks, including the trajectory, transport, area, thickness, and fate processes. It was concluded that the scope of spill trajectories was the largest under the influence of southwest winds, and it was the smallest under the influence of northwest winds; the transport of oil slicks was mainly affected by flood/ebb currents and oil slicks could reciprocate with flood/ebb currents; the spreading area of instantaneously spilled oil reached the maximum in the eastern spill location, under southwest winds, after spilling for 48 h. The minimum oil area appeared in the western and middle spill locations, which continuously spilled oil under the influence of northwest and northeast winds, respectively; the wind had a significant influence on drift and thickness of oil slicks, especially when the flow velocity was relatively small. The fate processes of oil volume on the sea surface gradually increase during the initial 10 h and subsequently the variation tendency is basically consistent with the instantaneous oil spill. The fate processes of evaporated, emulsified, and adsorbed oil volume of two types of oil spills are basically the same.Overall, the proposed model provides a reasonable method for the study of marine oil spills. Moreover, the simulation results will be helpful for controlling and handling of accidental oil spills efficiently.



---

*Source: 2736102-2018-01-17.xml*