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

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

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

[[partial derivative]z/[partial derivative]t] + [[partial derivative](hu)/[partial derivative]x] + [[partial derivative](hv)/[partial derivative]y] = 0. (1)

The equations for conservation of momentum are given by

[mathematical expression not reproducible], (2)

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, [A.sub.h] is the horizontal eddy viscosity coefficient, p is the density of water, and [[tau].sub.xs] and [[tau].sub.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 [[tau].sub.xs] and [[tau].sub.ys], which can be computed by the following empirical formula:

[[tau].sub.xs] = 0.125[C.sub.D][W.sub.x][absolute value of [??]] [[tau].sub.ys] = 0.125[C.sub.D][W.sub.y][absolute value of [??]], (3)

where [W.sub.x] and [W.sub.y] are the component velocities of the wind on the sea surface in the x and y directions, respectively, [??] is the wind velocity vector, and the wind drag coefficient [C.sub.D] can be obtained from the following Heaps empirical formula:

[mathematical expression not reproducible]. (4)

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:

[L.sub.1] = [K.sub.1][(gV[DELTA]).sup.1/4][t.sup.1/2] (5)

[L.sub.2] = [K.sub.2][[g(1 - [DELTA])[DELTA]].sup.1/6] [v.sup.-1/12.sub.w][t.sup.1/4][V.sup.1/4] (6)

[L.sub.3] = [K.sub.3][g.sup.-1/4][[delta].sup.1/2][[rho].sup.-3/4.sub.w][t.sup.3/4], (7)

where [DELTA] = 1 - [[rho].sub.0]/[[rho].sub.w], where [[rho].sub.0] and [[rho].sub.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, [delta] = [[delta].sub.wa] - [[delta].sub.ao] - [[delta].sub.ow], where [[delta].sub.wa], [[delta].sub.ao], and [[delta].sub.ow] are the water-air, oil slick-air, and oil slick-water surface tension coefficients, respectively, [v.sub.w] is the kinematic viscosity of water, [K.sub.1] = 1.35, [K.sub.2] = 1.60, and [K.sub.3] = 0.48 are experiment constants, t is the duration of the oil spill, V = [[summation].sup.n.sub.i=1][Q.sub.i](t)[DELTA]t[1 - k(t - i[DELTA]t)] is the spill volume, Q is the spill discharge, k is the synthetic attenuation coefficient, and n = t/[DELTA]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 Figure 2. The drift velocity [[??].sub.r] can be written as

[mathematical expression not reproducible], (8)

where [[??].sub.c] is the surface current velocity, [[??].sub.w] is the wind velocity at 10 m above the water surface, [K.sub.c] = 1 is the current factor, and [K.sub.w] = 0.035 is the wind drift factor [35].

The spreading thickness h can be determined from the mass conservation equation, as follows:

[[partial derivative](Ch)/[partial derivative]t] + [[partial derivative](uCh)/[partial derivative]x] + [[partial derivative](vCh)/[partial derivative]y] = ([[PHI].sub.s] + [[PHI].sub.b] + R), (9)

where C is the oil slick concentration, [[PHI].sub.s] and [[PHI].sub.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 thickness h can be computed from Fick's law:

h = [[Ve.sup.-kt]/2[pi][[sigma].sub.s][[sigma].sub.n]] exp(-[[s.sup.2]/2[[sigma].sub.s.sup.2]] - [[n.sup.2]/2[[sigma].sub.n.sup.2]]), (10)

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 [[sigma].sub.s] = [a.sub.s][t.sup.1.17] and [[sigma].sub.n] = (1/[square root of 10])[[sigma].sub.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 [DELTA]t > 0 and [X.sub.n] = X(n[DELTA]t), we have

[X.sub.n] = [X.sub.n-1] + [sigma][square root of [DELTA]t][W.sub.n] ([sigma] > 0), (11)

where {[W.sub.n]} are independent random numbers on N(0, 1) and the increment [X.sub.n] - [X.sub.n-k] depends only on k variables ([W.sub.n-k+1] ... [W.sub.n]) (k < n) corresponding to (n - k, n), so [X.sub.n] - [X.sub.n-k] follows the normal distribution N(0, [sigma][square root of k[DELTA]t]).

Specifically, supposing that position coordinates of an oil particle are R([t.sub.i]) and R([t.sub.i+1]) at times [t.sub.i] and [t.sub.i+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

[mathematical expression not reproducible], (12)

where [gamma] is the random number ranging from 0 to 1 and [[??].sub.l] is the drift vector in the period of [DELTA]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:

[mathematical expression not reproducible]. (13)

[[??].sub.s] is the spreading vector in the period of [DELTA]t, which is given by

[mathematical expression not reproducible]. (14)

The discrete form of transport distance of the labeled oil particle can be obtained by

[mathematical expression not reproducible]. (15)

In addition, [[??].sub.r] can be obtained from (8), and [??]([t.sub.i+1]) and [??]([t.sub.i]) 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]:

[F.sub.v] = ln [1 + (B[T.sub.G]/T)[theta] exp (A - B[T.sub.0]/T)]/(T/B[T.sub.G]), (16)

where [F.sub.v] is the volume fraction evaporated, A and B are the constants, usually selected as 6.3 and 10.3 for crude oils, respectively, [T.sub.G] is the slope of distillation curve, T is the surface temperature of the oil slick, [T.sub.0] = 542.6 - 30.275API + 1.565[API.sup.2] - 0.03439[API.sup.3] + 0.0002604[API.sup.4] is the initial boiling point [38], and API is the density of spilt oil following the classification of the American Petroleum Institute; [theta] = 0.0025[([u.sub.w] + 1).sup.0.78] x 2.437/h is the exposure coefficient of oil slick, and S(i) 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]

d[F.sub.w]/dt = [C.sub.1][([u.sub.w] + [u.sub.0]).sup.2] (1 - [C.sub.2][F.sub.w])Q(t), (17)

where [F.sub.w] is the emulsification fraction, [C.sub.1] is the absorption rate, usually selected as 2 x [10.sup.-6], [C.sub.2] is the water content, usually adopted as 1.33, [u.sub.0] 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 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)).

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.

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

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

3.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 x [10.sup.4] [km.sup.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)).

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, [K.sub.1], [M.sub.2], [O.sub.1], and [S.sub.2], 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 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]. The equations for these three criteria are shown as follows:

[mathematical expression not reproducible], (18)

where [[eta].sup.m.sub.i] are the modeled results and [[eta].sup.o.sub.i] 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 [+ or -]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 [+ or -]0.10 m/s at three stations (U1, U4, and U7), and, for the flow direction, the maximum RSME is 17.63[degrees] at Station U1 and the BIAS is below [+ or -]2[degrees] 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.

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

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

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

4.2. 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).

4.3. Area of Oil Slicks versus Time. Figures 20-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.

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.385 [km.sup.2] after 48 h. The minimum area occurred in the western and middle locations and reached 0.823 [km.sup.2], which was continuously spilling under the influence of northwest and northeast winds, respectively.

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

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

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.

https://doi.org/10.1155/2018/2736102

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was financially supported by the Opening Foundation of Key Laboratory of Marine Spill Oil Identification and Damage Assessment Technology, State Oceanic Administration (SOA). The authors greatly appreciate the assistance from Dr. Yangyang Li for subject research.

References

[1] T. M. Alves, E. Kokinou, and G. Zodiatis, "A three-step model to assess shoreline and offshore susceptibility to oil spills: the south aegean (crete) as an analogue for confined marine basins," Marine Pollution Bulletin, vol. 86, no. 1-2, pp. 443-457, 2014.

[2] T. M. Alves, E. Kokinou, G. Zodiatis, R. Lardner, C. Panagiotakis, and H. Radhakrishnan, "Modelling of oil spills in confined maritime basins: the case for early response in the Eastern Mediterranean Sea," Environmental Pollution, vol. 206, article no. 8069, pp. 390-399, 2015.

[3] T. M. Alves, E. Kokinou, G. Zodiatis, H. Radhakrishnan, C. Panagiotakis, and R. Lardner, "Multidisciplinary oil spill modeling to protect coastal communities and the environment of the Eastern Mediterranean Sea," Scientific Reports, vol. 6, Article ID 36882, 2016.

[4] H. A. Espedal and T. Wahl, "Satellite SAR oil spill detection using wind history information," International Journal of Remote Sensing, vol. 20, no. 1, pp. 49-65, 1999.

[5] C. Brekke and A. H. S. Solberg, "Oil spill detection by satellite remote sensing," Remote Sensing of Environment, vol. 95, no. 1, pp. 1-13, 2005.

[6] J. C. Dietrich, C. J. Trahan, M. T. Howard et al., "Surface trajectories of oil transport along the Northern Coastline of the Gulf of Mexico," Continental Shelf Research, vol. 41, pp. 17-47, 2012.

[7] H. Yang, B. Hong, and S. Chen, "Research and application process of marine oil spill models," Transactions of Oceanology and Limnology, vol. 2, pp. 156-163, 2007 (Chinese).

[8] X. Lou and S. G. Liu, "Review in theory and study of oil spill models," Environmental Science and Management, vol. 33, no. 10, article 61, pp. 33-37, 2008 (Chinese).

[9] G. Coppini, M. De Dominicis, G. Zodiatis et al., "Hindcast of oil-spill pollution during the Lebanon crisis in the Eastern Mediterranean, July-August 2006," Marine Pollution Bulletin, vol. 62, no. 1, pp. 140-153, 2011.

[10] G. Zodiatis, M. De Dominicis, L. Perivoliotis et al., "The mediterranean decision support system for marine safety dedicated to oil slicks predictions," Deep-Sea Research Part II-Topical Studies in Oceanography, vol. 133, pp. 4-20, 2016.

[11] W. J. Guo, Numerical simulation of oil spill based on POM, Dalian University of Technology, 2007 (Chinese).

[12] American Society of Civil Engineers, "State-of-the-art review of modelling transport and fate of oil spills," Journal of Hydraulic Engineering, vol. 122, no. 11, pp. 594-609, 1996.

[13] J. A. Galt, G. Y. Watabayashi, D. L. Payton, and J. C. Petersen, "Trajectory analysis for the Exxon Valdez: hindcast study," in Proceedings of the 1991 Oil Spill Conference, vol. 1991, pp. 629-634, Washington DC., Wash, USA.

[14] E. Howlett, K. Jayko, and M. L. Spaulding, "Interfacing real-time information with OILMAP," in Proceeding of the 16th Arctic and Marine Oil Spill Program Technical Seminar, pp. 517-527, Ottawa, Canada, 1993.

[15] M. Leech, M. Walker, M. Wiltshire et al., "OSIS--a windows-3 oil spill information-system," in Proceedings of the 16th Arctic and Marine Oil Spill Program (AMOP) Technical Seminar, pp. 549-572, Calgary, Canada.

[16] O. M. Aamo, M. Reed, and K. Downing, "Oil spill contingency and response (oscar) model system: sensitivity studies," in Proceedings of the 1997 International Oil Spill Conference--Improving Environmental Protection, vol. 1997, pp. 429-438, FT Lauderdale, FL, USA.

[17] J. K. Jolliff, S. Ladner, R. Crout et al., "Forecasting the ocean's optical environment using the BioCast system," Oceanography, vol. 27, no. 3, pp. 68-79, 2014.

[18] M. Skedsmo, R. Ayasse, N. Soleng, and M. Indregard, "Oil spill detection and response using satellite imagery, insight to technology and regulatory context," in Proceedings of the SPE International Conference and Exhibition on Health, Safety, Security, Environment, and Social Responsibility 2016, April 2016.

[19] M. Marghany, "Automatic Detection of Oil Spill Disasters Along Gulf of Mexico Using RADARSAT-2 SAR Data," Journal of the Indian Society of Remote Sensing, vol. 45, no. 3, pp. 503-511, 2017.

[20] J. K. O. Gjosteen, "Oil spreading in cold waters--A model suitable for broken ice," in Proceedings of the 11th International Offshore and Polar Engineering Conference (ISOPE '01), Stavanger, Norway, 2001.

[21] J. H. Wang and Y. M. Shen, "Development of an integrated model system to simulate transport and fate of oil spills in seas," Science China Technological Sciences, vol. 53, no. 9, pp. 2423-2434, 2010.

[22] J. H. Wang and Y. M. Shen, "Oil spill simulation system for complex terrain," Scientia Sinica (Technologica), vol. 40, no. 11, pp. 1367-1377, 2010 (Chinese).

[23] J. Wang and Y. Shen, "Modeling oil spills transportation in seas based on unstructured grid, finite-volume, wave-ocean model," Ocean Modelling, vol. 35, no. 4, pp. 332-344, 2010.

[24] J.-H. Wang and J.-S. Zhang, "Specification of turbulent diffusion by random walk method for oil dispersion modeling," Applied Mechanics and Materials, vol. 212-213, pp. 1161-1167, 2012.

[25] M. De Dominicis, N. Pinardi, G. Zodiatis, and R. Archetti, "MEDSLIK-II, a Lagrangian marine surface oil spill model for short-term forecasting-Part 2: numerical simulations and validations," Geoscientific Model Development, vol. 6, no. 6, pp. 1871-1888, 2013.

[26] Z. Deng, T. Yu, X. Jiang et al., "Bohai Sea oil spill model: a numerical case study," Marine Geophysical Research, vol. 34, no. 2, pp. 115-125, 2013.

[27] Y. Lu, X. Li, Q. Tian et al., "Progress in marine oil spill optical remote sensing: detected targets, spectral response characteristics, and theories," Marine Geodesy, vol. 36, no. 3, pp. 334-346, 2013.

[28] M. De Dominicis, S. Falchetti, F. Trotta et al., "A relocatable ocean model in support of environmental emergencies," Ocean Dynamics, vol. 64, no. 5, pp. 667-688, 2014.

[29] Y. C. Zeng, J. P. Yang, and C. W. Yu, "Mixed Euler-Lagrange approach to modeling fiber motion in high speed air flow," Applied Mathematical Modelling, vol. 29, no. 3, pp. 253-261, 2005.

[30] E. Capo, A. Orfila, J. M. Sayol et al., "Assessment of operational models in the Balearic Sea during a MEDESS-4MS experiment," Deep-Sea Research Part II: Topical Studies in Oceanography, vol. 133, pp. 118-131, 2016.

[31] W. Y. Tan, Computational Shallow Water Dynamics: Application of Finite Volume Method, Tsinghua University Press, Beijing, China, 1998.

[32] Y. F. Xu, Numerical Simulation of Wave and Analysis of Its Flow Field Structure [Master Thesis], Harbin Institute of Technology, 2013.

[33] J. A. Fay, The Spread of Oil Slicks on a Calm Sea/Oil on the Sea, Springer, 1969.

[34] H. M. Li, Numerical Simulation of the Spread-Diffusion Process of Oil Released from Seabed in Penglai 19-3 Oilfield Area [PhD Thesis], Ocean University of China, 2013 (Chinese).

[35] L. X. Huang, G. X. Zhang, and Z. Z. Wan, "The spread of oil in the sea," Chinese Journal of Environmental Engineering, vol. 3, no. 1, pp. 7-11, 1982.

[36] F. Yu, J. Li, S. Cui, Y. Zhao, Q. Feng, and G. Chen, "A hindcast method to simulate oil spill trajectories for the Bohai Sea, Northeast China," Ocean Engineering, vol. 124, pp. 363-370, 2016.

[37] W. Stiver and D. MacKay, "Evaporation rate of spills of hydrocarbons and petroleum mixtures," Environmental Science & Technology, vol. 18, no. 11, pp. 834-840, 1984.

[38] H. T. Shen and P. D. Yapa, "Oil slick transport in eivers," Journal of Hydraulic Engineering, vol. 114, no. 5, pp. 529-543, 1988.

[39] D. A. Mackay, A Mathematical Model of Oil Spill Behaviour, Ottawa, ontario, Canada, 1980.

[40] D. A. Mackay and I. Buist A, Mascarenhas R, Patersons. Oil Spill Processed and Models, Ottawa, Ontario, Canada, 1980.

[41] W. Q. Zhao and Z. H. Wu, "Determination of the dimension of an oil film by instantaneous oil slick on the sea surface," Journal of Chengdu University of Science and Technology, vol. 41, no. 5, pp. 63-72,1988 (Chinese).

[42] R. D. Ray, "A global ocean tide model from TOPEX/POSEIDON altimetry: GOT99. 2," Tech. Rep. 209478, NASA Technical Memorandum, 1999.

[43] TSDIWTE, Hydrometry Test Analysis Report of Tourism Construction Project in the Western Penglai Coast, Tianjin Research Institute for Water Transport Engineering, Ministry of Transportation, 2011.

[44] State Standard of the People's Republic of China, "Specifications for identification system of spilled oils on the sea (GB/T 21247-2007)," Tech. Rep. 21247, Standards Press of China, Beijing, China, 2007 (Chinese).

Darning Li (iD), Xingchen Tang (iD), Yanqing Li, Xiao Wang, and Hongqiang Zhang (iD)

State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, 135 Yaguan Rd. N., Jinnan District, Tianjin 300350, China

Correspondence should be addressed to Xingchen Tang; xingchentang1108@yeah.net

Received 5 September 2017; Revised 21 November 2017; Accepted 19 December 2017; Published 17 January 2018

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

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

Caption: 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).

Caption: 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).

Caption: 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).

Caption: Figure 6: The major (a) and minor (b) axes of the oil slick versus time.

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

Caption: Figure 8: (a) Bathymetry and (b) unstructured grids for the model domain.

Caption: Figure 9: Comparison of the tidal level between the modeled (solid line) and the observed (dots) results at three stations (H1, H2, and H3).

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

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

Caption: Figure 12: Distributions of the flow field at the times of ebb (a) and flood (b).

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

Caption: Figure 14: Wind rose diagram for Luanjiakou District in 2002-2006.

Caption: 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).

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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.

Caption: 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).

Caption: 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).
```Table 1: Comparison of the major axes scales of the oil slick.

Spill volume ([m.sup.3])     [10.sup.2]    [10.sup.3]    [10.sup.4]

Simulated values of this        15.33         28.14         45.1
paper (km)
Simulated values of [41]        12.67         25.91         43.54
(km)

Spill volume ([m.sup.3])     [10.sup.5]

Simulated values of this        64.83
paper (km)
Simulated values of [41]        65.49
(km)

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

Spill volume ([m.sup.3])           [10.sup.2]    [10.sup.3]

Simulated values of this              5.99          11.69
paper (km)
Simulated values of [41] (km)         5.18          11.74

Spill volume ([m.sup.3])           [10.sup.4]    [10.sup.5]

Simulated values of this              20.69         34.15
paper (km)
Simulated values of [41] (km)         21.10         34.04

Table 3: Comparison of the simulated and experimental results.

Item                      Initial size (cm)     Final size (cm)

Simulated results                15                   21
Experimental results             15                   22

Item                      Movement distance (m)     Movement time (s)

Simulated results                 1.17                     30
Experimental results               1.2                     30

Table 4: Statistical errors at tidal survey stations for model
verification.

Tidal level
Station    MAE (cm)     RSME (cm)     BIAS (cm)    Station

H1           9.18         11.04         -8.11         U1
H2           8.29         10.32         -6.83         U4
H3           10.02        12.10         -9.13         U7

Flow velocity
Station    MAE (m/s)      RSME (m/s)     BIAS (m/s)    MAE (deg)

H1            0.09           0.11           0.06         12.83
H2            0.06           0.08           -0.02        10.55
H3            0.07           0.09           -0.03        11.72

Flow direction
Station      RSME (deg)      BIAS (deg)

H1              17.63           1.63
H2              14.98           -1.98
H3              15.18           1.06

Table 5: Properties of the oil.

Density         Water content
Name                    (kg/[m.sup.3])     of emulsion (%)     API

Condensate oil               830.5               74           38.874
Low sulfur fuel oil           972                80           14.08

Table 6: Wind conditions of the model.

Wind direction       No wind      Southwest wind (SW)

Wind speed (m/s)        0                 4.9
Note                             Maximum wind direction

Wind direction         South wind (S)      Northwest wind (NW)

Wind speed (m/s)            2.0                    3.4
Note                  Wind Direction 1      Wind Direction 2

Wind direction        Northeast wind (NE)

Wind speed (m/s)              2.7
Note                   Wind Direction 3
```