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

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.11m/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.


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][2][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][5][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][2][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 advectiondiffusion 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][22][23][24][25][26][27][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 multiisland terrain or hydraulic structures.
In this study, an oil spill simulation method in a multiisland 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 (Figure 1).

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.

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 The equations for conservation of momentum are given by where  and V are the component velocities of current in the  and  directions, respectively,  is the water level,  is the gravitational acceleration,  is Chezy's friction coefficient, ℎ is the water depth,  is the Coriolis coefficient,  ℎ is the horizontal eddy viscosity coefficient,  is the density of water, and   and   are the wind stress components on the sea surface in the  and  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 twodimensional, 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   and   , which can be computed by the following empirical formula: where   and   are the component velocities of the wind on the sea surface in the  and  directions, respectively,  →  is the wind velocity vector, and the wind drag coefficient   can be obtained from the following Heaps empirical formula:  [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: Gravity-viscous spreading phase: Surface tension viscous spreading phase: where Δ = 1 −  0 /  , where  0 and   are the density of oil particles and water, respectively,  is the spreading diameter of the oil slick, and the subscripts denote different phases,  =   −   −   , where   ,   , and   are the water-air, oil slick-air, and oil slick-water surface tension coefficients, respectively,   is the kinematic viscosity of water,  1 = 1.35,  2 = 1.60, and  3 = 0.48 are experiment constants,  is the duration of the oil spill,  = ∑  =1   ()Δ[1 − ( − Δ)] is the spill volume,  is the spill discharge,  is the synthetic attenuation coefficient, and  = /Δ is the number of oil spill times.

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 Figure 2. The drift velocity  ⇀   can be written as where  →   is the surface current velocity,  →   is the wind velocity at 10 m above the water surface,   = 1 is the current factor, and   = 0.035 is the wind drift factor [35].
The spreading thickness ℎ can be determined from the mass conservation equation, as follows: where  is the oil slick concentration, Φ  and Φ  are the spill flux of the upper and lower surface of the oil slick, respectively, and  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 thickness ℎ can be computed from Fick's law: where  and  are the natural coordinates in the oil slick drift direction and the direction perpendicular to , respectively,  is the attenuation coefficient of oil, and   =    1.17 and   = (1/ √ 10)  are the standard deviation of the oil slick thickness in the  and  directions, respectively.

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 Δ > 0 and   = (Δ), we have where {  } are independent random numbers on (0, 1) and the increment   −  − depends only on  variables ( −+1 ⋅ ⋅ ⋅   ) ( < ) corresponding to ( − , ), so   −  − follows the normal distribution (0,  √ Δ).Specifically, supposing that position coordinates of an oil particle are (  ) and ( +1 ) at times   and  +1 , respectively, and  is the movement distance of an oil particle under the actions of spreading, drift, and so on, we will then have where  is the random number ranging from 0 to 1 and  →   is the drift vector in the period of Δ, 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: →   is the spreading vector in the period of Δ, which is given by The discrete form of transport distance of the labeled oil particle can be obtained by In addition,  →   can be obtained from (8), and  →  ( +1 ) and  →  (  ) 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.

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]: where  V is the volume fraction evaporated,  and  are the constants, usually selected as 6.3 and 10.3 for crude oils, respectively,   is the slope of distillation curve,  is the surface temperature of the oil slick,  0 = 542.6 − 30.275API + 1.565API 2 − 0.03439API 3 + 0.0002604API 4 is the initial boiling point [38], and API is the density of spilt oil following the classification of the American Petroleum Institute;  = 0.0025(  + 1) 0.78 × 2.437/ℎ is the exposure coefficient of oil slick, and () 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] where   is the emulsification fraction,  1 is the absorption rate, usually selected as 2 × 10 −6 ,  2 is the water content, usually adopted as 1.33,  0 is the emulsification correction factor in the ocean environment, and () is the emulsifying amount of the oil slick.

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 Figure 3).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.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 Figure 4).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.
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)).

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.
(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.

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 Figure 1).

Study Area.
The model domain and its bathymetry are shown in Figure 8(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 × 10 4 km 2 .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)).

Boundary Condition.
To account for the lack of observational data, the astronomical tide, we induced the tidal   [42], so that the tidal levels processes can be obtained at the open sea boundaries.Moreover, observational data are used for the landward boundaries.

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 validation results of the tidal level are shown in Figure 9, 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.
There are many diurnal tide survey stations (see Figure 1).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.
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] where    are the modeled results and    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).
The distributions of the flow field at ebb and flood periods are shown in Figure 12.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.
In summary, the hydrodynamic field can serve as the basis for studying marine oil spills in our study area.

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 Figure 13.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.

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.
In 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.    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.

Movement Process of Oil Slicks. Figures 16 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.
Figures 18 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.
From Section 2.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    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.From Figures 20-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.385km 2 after 48 h.The minimum area occurred in the western and middle locations and reached 0.823 km 2 , which was continuously spilling under the influence of northwest and northeast winds, respectively.

Thickness of Oil Slicks versus Time. Figures 25 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.

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.Figure 27 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 28 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.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.

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

Figure 1 :
Figure 1: Map of the Luanjiakou District near the Port of Yantai and survey stations.

Figure 2 :
Figure 2: Velocity vector of oil slick drift.

Figure 3 :
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).

Figure 4 :
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).

Figure 5 :Figure 6 :
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).

Figure 7 :
Figure 7: Comparison of the flume experiment (a, c) and the simulated result (b) of the spreading and drift of the oil slick.

Figure 11 :Figure 12 :
Figure 11: Comparison of flow direction between the modeled (solid line) and the observed (dots) results at three stations (U1, U4, and U7).

Table 4 :Figure 13 :
Figure 13: Comparison between the experimental result (a) and the modeled result (b) of the concentration diffusion of the oil slick.

4. 1 .
Spill Trajectories.The trajectories of instantaneous oil spills from the western portion of the channel under five wind conditions are shown in Figure 15.In the figure, it

Figure 15 :
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).

Figure 16 :
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. ) .

4. 3 .Figure 17 :Figure 18 :
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.

Figure 19 :Figure 20 :
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.

Figure 21 :Figure 22 :
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.

Figure 23 :Figure 24 :
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.

Figure 25 :
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.

Figure 26 :
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.

Table 1 :
Comparison of the major axes scales of the oil slick.

Table 2 :
Comparison of the minor axes scales of the oil slick.

Table 3 :
Comparison of the simulated and experimental results.the three open boundaries.Four main constituents in this domain are considered, that is, K 1 , M 2 , O 1 , and S 2 , whose harmonic constants can be derived from the global ocean tide model from the United States Department of the Navy

Table 5 :
Properties of the oil.

Table 6 :
Wind conditions of the model.