Numerical Simulation of Oil Spill in Ocean

Abstract

The spreading of oil in an open ocean may cause serious damage to a marine environmental system. Thus, an accurate prediction of oil spill is very important to minimize coastal damage due to unexpected oil spill accident. The movement of oil may be represented with a numerical model that solves an advection-diffusion-reaction equation with a proper numerical scheme. In this study, the spilled oil dispersion model has been established in consideration of tide and tidal currents simultaneously. The velocity components in the advection-diffusion-reaction equation are obtained from the shallow-water equations. The accuracy of the model is verified by applying it to a simple but significant problem. The results produced by the model agree with corresponding analytical solutions and field-observed data. The model is then applied to predict the spreading of an oil spill in a real coastal environment.

1. Introduction

In an ocean environment, a spilled oil due to an unexpected accident may cause a serious damage. Thus, an accurate prediction of behaviors of split oil is very significant to keep coastal environmental system. The area of oil spreading can be predicted numerically by solving the proper equations governing the flow field and the associated mass transport phenomenon. The most reasonable choice is probably an advection-diffusion equation consisted of both advection and diffusion. The equation is difficult to solve numerically because it is represented by both hyperbolic- and parabolic-type partial differential equations. A standard split operator approach is a plausible and practical choice for solving the advection-diffusion-type partial differential equation [1]. The hyperbolic (advection) and parabolic (diffusion) components of the equation are solved separately using methods that properly and simply simulate the physical behaviors. Numerical solutions obtained from the advection and diffusion equations are then combined.

An oil spill accident is regarded as a kind of a disaster because it causes not only fatal destruction of the marine environment but also enormous cost of the disaster prevention and the damage compensation. Thus, it is important to accurately predict the spread range of the spilled oil as an early stage countermeasure against a disaster. In the last three decades, many investigators have studied the transport processes of oil spills based on the trajectory method [2–4]. Those methods have been applied in river-lake system [5–7] and seas [8–11]. Some commercial oil spill models, such as, COZOIL [12], NOAA [13], OILMAP [14], WOSM [15], have been used to determine the oil movement and distribution in the water body. However, transport of oil spills has been conducted considering tidal currents simultaneously in only few researches. Furthermore, oil spills processes in a field of strong tide and tidal currents have not been investigated.

In this study, the spilled oil dispersion model has been established in consideration of tide and tidal currents simultaneously. The Hebei Spirit oil spill that occurred on December 7, 2007 is the largest oil spill accident occurred in the Yellow Sea. In the Yellow Sea, tidal currents are very strong and should be considered to investigate related coastal processes. Therefore, the accuracy of the model for predicting tide and tidal currents is very important to investigate oil spreading in this area. Subsequently, by computing the diffusion distribution on the Hebei Spirit oil spill considering tidal currents simultaneously, verification has been made through comparison of the diffusion distribution and the field-observed data between at 8 p.m. on December 7, 2007 when the spilled oil flowed into the whole Manripo and Sindu-ri seashore and at 11 a.m. on December 11, 2007 when the satellite photographs exist. Also, this study aims to prepare for a possible accident in the future and provide the basic materials for establishing the disaster prevention of the oil spill pollution.

In the following section, the governing equations are described first. The numerical model is then presented and detailed description is followed. Numerical simulations of tide and tidal currents are conducted to test the applicability of the spilled oil spreading model in the Yellow Sea, and the verification has been made comparing numerical results with the field observed data. Also, the model is verified conducting the numerical calculations on the simple diffusion distribution problem. The model is then applied to predict the Hebei Spirit oil spill in a real coastal environment. Finally, concluding remarks are made.

2. Governing Equations

The behaviors of tide and tidal currents may be described by the following nonlinear shallow-water equations [16, 17]:𝜕𝜂𝜕𝑃𝜕𝑡+𝜕𝜕𝑃𝑥+𝜕𝜕𝑄𝑦=0,(2.1)𝜕𝜕𝑄𝑡+𝜕𝜕𝑥𝑃𝐻²+𝜕𝜕𝑦𝑃𝐻𝑄+𝑔𝐻𝜕𝜕𝑥𝜂+𝜏𝜌^𝑥=0,(2.2)𝜕𝑡+𝜕𝜕𝑥𝑃𝐻𝑄+𝜕𝜕𝑦𝑄𝐻²+𝑔𝐻𝜕𝜕𝜂𝑦+𝜏𝜌^𝑦=0.(2.3) In (2.1)–(2.3), 𝜂 represents the free surface displacement and 𝐻 is the total water depth defined as𝐻=𝜂+ℎ with ℎ being a local still water depth (Figure 1), 𝑃 and 𝑄 are volume flux components along the 𝑥- and 𝑦-axis directions defined as 𝑃=𝑢𝐻 and 𝑄=𝑣𝐻 with 𝑢 and 𝑣 being the depth-integrated velocity components in the 𝑥- and 𝑦-axis directions, and 𝜏_𝑥 and 𝜏_𝑦 are bottom frictional effects of the 𝑥- and 𝑦-axis directions, respectively.

Figure 1 

A definition sketch of the free surface displacement and local water depth.

The study is focused to investigate the oil spreading in a limited area in the Southern part of the Korean Peninsula and thus the Cartesian coordinate system is then employed. The boundary conditions along the offshore are provided by a large-scale numerical model [18]. The bottom frictional effects may play a significant role in the very shallow zone. Thus, the Manning’s empirical formula is employed in the study given as:𝜏_𝑥𝜌=𝐻𝑔𝑛^27/3𝜏_𝑦𝑃𝑃²+𝑄²,𝜌=𝐻𝑔𝑛^27/3𝑄𝑃²+𝑄².(2.4) In (2.4), 𝑛 is the Manning’s roughness coefficient. In this study, the Manning coefficient of 𝑛=0.013 has been used.

It is well known that the change of sea level is very complicated and affected by many physical phenomena such as waves, tides, and others. Among these, the tide may be the most dominant force for the oil spreading in the Yellow Sea. In this study, therefore, the tide is only investigated in a very limited area without considering any interaction with other physical phenomena.

In this study, a source or sink can be described to a following equation [19]:𝑆(𝐶)=𝑑𝑆𝑑𝑡=𝑑𝑑𝐶𝑆𝑑𝑑𝐶𝑡,(2.5) in which 𝐶 represents the concentration of the oil in the sea. The source or sink 𝑆(𝐶) is assumed to be nonlinear and is linearized using the Newton-Raphson iterative technique [20] as:𝑆=𝑆^∗+𝑑𝑆𝑑𝐶|||^∗𝐶−𝐶^∗,(2.6) in which the superscript (∗) represents a previous step and the sedimentation rate, 𝑘 can be derived as the following relation [21]:𝑘=𝑑𝑆𝑑𝐶|||^∗.(2.7)

By applying an operator splitting approach to (2.1), the following advection equation can first be obtained by:𝜕𝐶𝜕𝑡+𝜕(𝜕𝐶𝑥𝑈)+𝜕(𝜕𝐶𝑦𝑉)=0,(2.8) and the diffusion-reaction equation can then be derived by𝜕𝐶𝜕𝑡=𝐻1𝜕𝜕𝑥𝐻𝐷^𝑥𝑥𝜕𝐶𝜕𝑥+𝐻𝐷^𝑥𝑦𝜕𝐶𝜕𝑦+𝐻1𝜕𝜕𝑦𝐻𝐷^𝑦𝑥𝜕𝐶𝜕𝑥+𝐻𝐷^𝑦𝑦𝜕𝐶𝜕𝑦+𝑆,(2.9) in which 𝐷_𝑥𝑥, 𝐷_𝑥𝑦, 𝐷_𝑦𝑥, 𝐷_𝑦𝑦 represent the diffusion coefficients.

3. Numerical Model

The continuity and momentum equations given by (2.1)–(2.3) are solved by the finite difference approximations. The detailed description of the finite difference method can be found in literatures [2, 3] and is not repeated here. Firstly, the advection equation can be discretized by using the SOWMAC (second-order wave equation method for advective calculation) scheme based on the characteristic method and is given as [22](1−+𝜒𝛼)₂𝐶^{𝑛+1𝑖−1}−2𝐶^𝑛𝑖−1+𝐶^𝑛𝑖−1+𝜒𝐶^𝑖𝑛+1−2𝐶^𝑛𝑖+𝐶^𝑖𝑛−1𝜃𝐶^𝑛+1𝑖+1−2𝐶^𝑖𝑛+1+𝐶^{𝑛+1𝑖−1}+(1−𝜃)𝐶^𝑛𝑖+1−2𝐶^𝑛𝑖+𝐶^𝑛𝑖−1=0,(3.1) in which 𝛼 is the Courant number defined as 𝛼=𝑈Δ𝑡/Δ𝑥, and superscript 𝑛 and subscript 𝑖 are time level and computational point, respectively. 𝜒 and 𝜃 are weighting factors of the finite difference scheme.

An unnecessary wave propagates upstream if the initial concentration distributions at both 𝑡=0 and 𝑡=−Δ𝑡 are not consistent with the downstream advection progress. To overcome this difficulty, the values of the concentration at time step (𝑛−1) can be employed as shown in Figure 2 instead of time step (𝑛−1). Drawing the characteristic curves which go downstream through points (𝑖,𝑛) and (𝑖+1,𝑛), time step (𝑛−1) can be introduced corresponding to time (𝑛Δ𝑡−Δ𝑥/𝑈) when each characteristic curve intersects the (𝑖−1-)axis or 𝑖-axis. According to the concept of characteristics for a pure advection, the concentrations are conserved between two time steps and are given as follows:𝐶𝐶^{(𝑛−1)𝑖−1}_{𝑖(𝑛−1)}=𝐶^𝑛𝑖,=𝐶^𝑛𝑖+1.(3.2)

Figure 2 

Estimation of concentration at a previous time step.

Using (3.2), (3.1) can be rewritten as2(1−=𝜒𝛼)(𝛼𝐶+1^{𝑛+1𝑖−1})−𝜃2𝐶𝐶^𝑛𝑖−1_𝑛+1𝑖+1+𝛼𝐶^𝑛𝑖+2𝜒𝐶^𝑖𝑛+1−(𝛼+1)𝐶^𝑛𝑖+𝛼𝐶^𝑛𝑖+1−2𝐶^𝑖𝑛+1+𝐶^{𝑛+1𝑖−1}+𝛼(𝛼+1)(1−𝜃)𝐶^𝑛𝑖+1−2𝐶^𝑛𝑖+𝐶^𝑛𝑖−1.(3.3) The optimal values of 𝜒 and 𝜃 in (3.3) do not seem to be constant, but instead are dependent on the Courant number 𝛼. To obtain the functional relationship of 𝜒 and 𝜃 to 𝛼, the Taylor series expansion on (3.3) may be employed [22].

The final expression of the advection equation can be written as follows:𝑝₁𝐶^{𝑛+1𝑖−1}+𝑝₂𝐶^𝑖𝑛+1+𝑝₃𝐶^𝑛+1𝑖+1=𝑝₄𝐶^𝑛𝑖−1+𝑝₅𝐶^𝑛𝑖+𝑝₆𝐶^𝑛𝑖+1,(3.4) in which the coefficients are given by𝑝𝑝₁₂=0.3776𝛼₀₊+0.3152𝛼₀₋−0.5467𝛼₁₊+0.4843𝛼₁₋+0.1691𝛼²,𝑝₃=1.3072+0.0624|𝛼|−0.3382𝛼²,𝑝₄=0.3152𝛼₀₊+0.3776𝛼₀₋+0.4843𝛼₁₊−0.5467𝛼₁₋+0.1691𝛼²,𝑝₅=0.3776𝛼₀₊+0.3152𝛼₀₋+0.5157𝛼₁₊−0.4533𝛼₁₋+0.1381𝛼²,𝑝₆=1.3072−0.0624|𝛼|−0.2762𝛼²,=0.3152𝛼₀₊+0.3776𝛼₀₋−0.4553𝛼₁₊+0.5157𝛼₁₋+0.1383𝛼².(3.5) The values of 𝛼₀₊, 𝛼₀₋, 𝛼₁₊ and 𝛼₁₋ are defined as:𝛼₀₊=AINT𝛼+1𝛼|𝛼|+1₀₋,=AINT1−𝛼𝛼1+|𝛼|₁₊,=|𝛼|+𝛼𝛼₁₋2,=|||𝛼−|𝛼|2|||,(3.6) in which AINT is one of the intrinsic functions in FORTRAN that carries out the function of truncating decimals.

Secondly, the diffusion-reaction equation can be written as:𝜕𝐶𝜕𝑡=𝜕𝜕𝑥𝐷^𝑥𝑥𝜕𝐶𝜕𝑥+𝐷^𝑥𝑦𝜕𝐶𝜕𝑦+𝜕𝜕𝑦𝐷^𝑦𝑥𝜕𝐶𝜕𝑥+𝐷^𝑦𝑦𝜕𝐶𝜕𝑦+𝑘𝐶.(3.7) The diffusion-reaction equation is discretized by a three-level locally implicit scheme with reference to the computational grid, as shown in Figure 3 [23].

Figure 3 

Computational mesh for the diffusion-reaction term.

Equation (3.7) contains two dimensions and is different from the form derived in Hobson et al. [23]. However, their numerical scheme can be easily extended to the 𝑦-direction in a similar manner and the final form of discretized equations using the finite difference approximation is given as:1+(4𝛼+4)𝐾+(16𝛼−4𝛽)𝐾²𝐶^{𝑛+1𝑖,𝑗}−𝐶^{𝑛−1𝑖,𝑗}=𝐷𝐾4𝐶_{𝑛−1𝑖+2,𝑗}2Δ𝑡+𝐶^{𝑛−1𝑖,𝑗+2}+𝐶^{𝑛−1𝑖−2,𝑗}+𝐶^{𝑛−1𝑖,𝑗−2}−4𝐶_𝑖,𝑗+2𝐶_{𝑛−1𝑖+1,𝑗+1}(2Δ𝑥)²+𝐶^{𝑛−1𝑖+1,𝑗−1}+𝐶^{𝑛−1𝑖−1,𝑗+1}+𝐶^{𝑛−1𝑖−1,𝑗−1}−4𝐶^{𝑛−1𝑖,𝑗}−4(1−𝛼)𝐶_{𝑛−1𝑖+1,𝑗+1}(2Δ𝑥)²+𝐶^{𝑛−1𝑖+1,𝑗−1}(Δ𝑥)²+𝐶^{𝑛−1𝑖−1,𝑗+1}+𝐶^{𝑛−1𝑖−1,𝑗−1}−4𝐶^{𝑛−1𝑖,𝑗}+𝐶(Δ𝑥)²_{𝑛−1𝑖+1,𝑗+1}+𝐶^{𝑛−1𝑖+1,𝑗−1}+𝐶^{𝑛−1𝑖−1,𝑗+1}+𝐶^{𝑛−1𝑖−1,𝑗−1}−4𝐶^{𝑛−1𝑖,𝑗}+2Δ𝑡(4𝐾+4𝛼𝐾+1)𝑘𝐶(Δ,𝑥)²(3.8) in which 𝐾=2𝐷Δ𝑡/(Δ𝑥)² is used and the diffusion coefficients are assumed as a constant value, 𝐷. In (3.8), 𝛼 and 𝛽 are the artificial values which are determined to give the correct amount of diffusion and avoid any time lag problems. In general, a whole family of diffusion schemes can be obtained by changing the values of 𝛼 and 𝛽. The optimum scheme of this type is the one which exhibits no time lag and is the most numerically stable.

In this study, (3.8) is solved by a two-dimensional three-level locally implicit scheme with second-order accuracy [23]. Equation (3.8) can be rewritten for 𝐶^{𝑛+1𝑖,𝑗} as follows:𝐶^{𝑛+1𝑖,𝑗}=𝐴+𝐵𝐶𝐶^{𝑛−1𝑖+2,𝑗}+𝐶^{𝑛−1𝑖,𝑗+2}_{𝑛−1𝑖+1,𝑗+1}+𝐶^{𝑛−1𝑖−2,𝑗}+𝐶^{𝑛−1𝑖,𝑗−2}+𝐸𝐶_{𝑛−1𝑖+1,𝑗}+𝐶^{𝑛−1𝑖+1,𝑗−1}+𝐶^{𝑛−1𝑖−1,𝑗+1}+𝐶^{𝑛−1𝑖−1,𝑗−1}+𝐶^{𝑛−1𝑖,𝑗+1}+𝐶^{𝑛−1𝑖−1,𝑗}+𝐶^{𝑛−1𝑖,𝑗−1}+𝐹𝐶^{𝑛−1𝑖,𝑗}+𝐺(𝑘𝐶),(3.9) in which 𝐴, 𝐵, 𝐸, 𝐹 and 𝐺 are given as:𝐴=𝐾²𝐵=1+(4𝛼+4)𝐾2𝐾+₂(16𝛼−4𝛽)𝐾²,𝐸=1+(4𝛼𝐾+[14−)𝐾4+(1(−16𝛼𝛼)−𝐾]4𝛽)𝐾²,𝐹=1+(14+𝛼+4𝛼4𝐾)𝐾++4((116−𝛼𝛽−)𝐾4𝛽₂)𝐾²,𝐺=1+(42𝛼Δ+𝑡(44)𝐾𝐾++4(𝛼1𝐾6𝛼+−14)𝛽)𝐾²,1+(4𝛼+4)𝐾+(16𝛼−4𝛽)𝐾².(3.10)

Since a first-order sink and zero-order reaction cannot generate numerical stability alone, this study is concerned only with the presence of a first-order term. Thus, the zero-order sedimentation rate, 𝑘₀, is assumed to be zero, and the first-order sedimentation rate of a suspended material, 𝑘₁, is described as:𝑘₁=−𝜆𝑊^𝑠𝐻,(3.11) and a first-order sink is defined as:sink=𝑘₁𝐶,(3.12) in which 𝜆 is a coefficient for calculation of the concentration at the bottom and, for convenience, is assumed to be 1, and 𝑊_𝑠 is the sedimentation rate, assumed to be 0.08 mm/sec.

4. Verification of the Numerical Model

The numerical model developed in this study is employed to simulate the change of tide and tidal currents near the accident area of the Hebei Spirit oil spill occurred on December 7, 2007. The Hebei Spirit oil spill is the largest oil spill accident occurred in the Yellow Sea located between Korea and China. In general, tide effects play a dominant role in coastal processes in the Yellow Sea, which is well known for the strongest tidal currents in the world. Therefore, the accuracy of the model for predicting tide and tidal currents is very important to investigate oil spreading in this area. To verify the model, free surface elevations and tidal ellipses at several points are computed and compared with the corresponding observed data.

Figure 4 shows the computational domain and sea level observation points. Sea bottom topography is also shown in the figure. Numerical simulations of tide and tidal currents were conducted during 16 days (December 7 ~ 22, 2007) and the harmonic analysis for computed results is presented. Table 1 displays a comparison of observed data and computed results using harmonic analysis. Computed results represent both amplitudes and phase differences very well. Figure 5 shows a comparison of numerically obtained free surface profiles and corresponding observed measurements at points of T-1, T-2, and T-3. The agreement between the numerical solutions and field measurements is quite reasonable.

Table 1 

Observed data and computed results using harmonic analysis.

Figure 4 

Computational domain and sea bottom topography.

(a) Free surface profile at Taean (T-1)

(b) Free surface profile at Sinjin-do (T-2)

(c) Free surface profile at Baeksajang (T-3)

Figure 5 

Comparison of observed and computed free surfaces of tides.

Tidal currents are very strong in the Yellow Sea and those may play an important role in oil spreading. In Figure 6, tidal ellipses at points of C-1 and C-2 are computed and compared with field-observed data to verify the accuracy of the numerical model. The numerical model represents well-strong tidal currents observed in this area.

(a) C-1 (observed data)

(b) C-1 (numerical model)

(c) C-2 (observed data)

(d) C-2 (numerical model)

Figure 6 

Comparison of observed and predicted tidal ellipses.

5. Application of the Model

The spilled oil dispersion model as established in previous sections is applied to predict the diffusion distribution on the Hebei Spirit oil spill. Subsequently, by conducting the numerical calculations on the diffusion distribution considering tide and tidal currents simultaneously, the prediction has been made through the comparison of the diffusion distribution observation between at 8 p.m. on December 7, 2007 when the spilled oil flowed into the whole Manripo and Sindu-ri seashore and at 11 a.m. on December 11, 2007 when the satellite photographs exist (Figure 7). In Figure 7, a target area of the numerical simulation is also described.

Figure 7 

Satellite image of the Hebei Spirit oil spill (December 11, 2007).

Figures 8(a) and 8(b) show distribution of simulated tidal currents field, 8(a) flood phase and 8(b) ebb phase, where the Hebei Spirit was stranded on duration of a spring tide. The results show that the tide current reaches the maximum speed in the northern area from the accident and mainly flows not perpendicular to the coastline but parallel to that. As a result, strong invasion of the oil spill to the coastline is not represented in the numerical simulation.

(a) Flood current at a spring tide

(b) Ebb current at a spring tide

Figure 8 

Distribution of tidal currents near the target area.

Figures 9(a) and 9(b) show numerically obtained diffusion distribution at 8 p.m. on December 7, 2007 and 11 a.m. on December 11, 2007. Figure 9(a) shows diffusion distribution after 13 hours since the Hebei Spirit had been stranded. The south-western currents on a flood phase and north-eastern currents on ebb phase are well reproduced and those are well reflected in the shape of oil distribution. Therefore, the numerical result shows similar distribution of the oil spreading with the satellite photograph (Figure 7). Since there is lack of quantitative observed data near the target area, numerical results are compared with a satellite image only in a qualitative viewpoint. On the other hand, by comparing the satellite photograph and Figure 9(b), numerical results do not represent strong invasion of the oil spill to the coastline. It may be occurred because effects of wind wave were neglected in the numerical simulation. Nevertheless, it is obvious that the model represents well the whole processes of the oil spreading, and the model can be applied to predict the oil spreading in a real coastal environment.

(a) 8 p.m. on December 7, 2007

(b) 11 a.m. on December 11, 2007

Figure 9 

Computed diffusion distribution of the Hebei Spirit oil spreading.

6. Concluding Remarks

In this study, a numerical model is employed to solve a two-dimensional advection-diffusion-reaction equation. The model is based on a standard split operator (fractional step) approach. Thus, the hyperbolic (advection) and parabolic (diffusion) portions of the equation are solved separately by using techniques describing properly the physical behaviors of each. In the model, the advection step is solved by using the SOWMAC scheme, while the diffusion-reaction step is done by a three-level locally implicit scheme.

The numerical model is first applied to an idealized problem to verify its accuracy. The model produces results agreeable with corresponding analytical solutions and field observed data. The model is then employed to investigate the behavior of the oil spreading in a coastal environment. The model yields reasonable results. The employed model could be used to forecast the behaviors of oil spreading in various practical situations.

Acknowledgment

This research was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education, Science and Technology (2011-0015386).