Chinese Journal of Oceanology and Limnology   2016, Vol. 34 Issue(1): 245-263     PDF       
http://dx.doi.org/10.1007/s00343-015-4230-7
Institute of Oceanology, Chinese Academy of Sciences
0

Article Information

CHEN Haibo (陈海波)1,2, AN Wei (安伟)2, YOU Yunxiang (尤云祥)1, LEI Fanghui (雷方辉)3, ZHAO Yupeng (赵宇鹏)2, LI Jianwei (李建伟)2
Modeling underwater transport of oil spilled from deepwater area in the South China Sea
Chinese Journal of Oceanology and Limnology, 2016, 34(1): 245-263
http://dx.doi.org/10.1007/s00343-015-4230-7

Article History

Received Aug. 24, 2014
accepted in principle Dec. 13, 2014;
accepted for publication Mar. 6, 2015
Modeling underwater transport of oil spilled from deepwater area in the South China Sea
CHEN Haibo (陈海波)1,2 , AN Wei (安伟)2, YOU Yunxiang (尤云祥)1, LEI Fanghui (雷方辉)3, ZHAO Yupeng (赵宇鹏)2, LI Jianwei (李建伟)2       
1 State Key Laboratory of Ocean Engineering, Shanghai Jiaotong University, Shanghai 200240, China;
2 China Offshore Environmental Service Co. Ltd ., Tianjin 300452, China;
3 China National Offshore Oil Corporation General Research Institute, Beijing 100027, China
ABSTRACT:Based on a Lagrangian integral technique and Lagrangian particle-tracking technique, a numerical model was developed to simulate the underwater transport of oil from a deepwater spill. This model comprises two submodels: a plume dynamics model and an advection-diff usion model. The former is used to simulate the stages dominated by the initial jet momentum and plume buoyancy of the spilled oil, while the latter is used to simulate the stage dominated by the ambient current and turbulence. The model validity was verifi ed through comparisons of the model predictions with experimental data from several laboratory fl ume experiments and a fi eld experiment. To demonstrate the capability of the model further, it was applied to the simulation of a hypothetical oil spill occurring at the seabed of a deepwater oil/gas fi eld in the South China Sea. The results of the simulation would be useful for contingency planning with regard to the emergency response to an underwater oil spill.
Keywordsunderwater oil spill     numerical simulation     contingency planning     deepwater oil/gas fi eld     South China Sea    
 1 INTRODUCTION

The worldwide increase in dem and for oil and dwindling onshore reserves have meant that off shore oil production has increased significantly since the 1990s. Particularly in China, rapid economic growth has caused a continuous increase in oil consumption. It has been reported that in 2013, the quantity of oil imported by China was 2.82×10 8 tons, representing year-on-year growth of 4.03%(http://www.chinairn.com/news/20140810/100701651.shtml). There are rich oil resources in the South China Sea(SCS)with some estimates placing the reserves at between 2.3×1010 and 3.0×1010 tons; however, about 70% lies in deepwater areas. Advances in modern technology mean that the production of oil from deepwater wells is commercially viable. In 2012, the first deepwater semi-submersible drilling vessel called “Off shore Oil 981”, owned by the China National Off shore Oil Corporation, was formally put into operation in the SCS. Its maximum operating depth is 3 000 m, indicating that China has the capability to conduct independent deepwater oil/gas exploration.

Because of the complex hydrodynamic environment of deepwater areas, the expansion of exploitation increases the possibility of accidental oil releases from well blowouts and pipeline or riser ruptures. There have been major oil spill accidents that have brought worldwide attention to the problems caused by such occurrences, e.g., the 2010 Deepwater Horizon oil spill in the Gulf of Mexico and the 2011 Penglai 19-3 oil field oil spill in the Bohai Sea of China. Such large oil spill accidents not only cause extensive damage to marine ecosystem and wildlife habitats, but they also endanger fishing and tourism industries. The complex hydrodynamic environment and adverse seasonal weather conditions in the SCS present the potential for accidental oil releases. As the SCS borders various countries and some heavily populated areas, the complex political environment means that an oil spill accident could easily become a serious international incident.

In recognition of the potential for oil spill accidents, many governmental agencies have to prepare emergency response contingency plans. Underst and ing the underwater behavior of oil through simulations is critical to any contingency plan, especially when the accident occurs in a deepwater area. Thus, a good oil spill model is very useful for planners, because model predictions may often be the only information available. A contingency plan combined with an oil spill model can help the operators better underst and how spilled oil will disperse as it moves up through the water column, how to track it, and how to clean it up once it reaches the surface. Therefore, since the mid-1970s, many oil spill models have been developed to predict the transport and fate of oil spilled in the marine environment. Earlier studies were concerned mainly with surface, near surface, and shallow water spills(Topham, 1975; McDougall, 1978; Fanneløp and Sjøen, 1980; Johansen, 1984; Elliott, 1986; Shen et al., 1986; Lee and Cheung, 1990; Rye et al., 1996; Yapa and Zheng, 1997; Lonin, 1999; Yapa et al., 1999). Subsequent eff orts were largely devoted to developing an enhanced comprehensive oil spill model that could describe deepwater oil spills. Two such well-known models are DeepBlow(Johansen, 2000) and ADMS/CDOG(Zheng et al., 2003). The most recent developments in underwater oil spill modeling have been summarized by Yapa et al.(2012).

In June 2000, the large-scale field experiment “DeepSpill” was conducted as a joint industry project in the deepwater(840 m)of the Norwegian Sea(Johansen et al., 2001), from which a considerable amount of experimental data were made available. Subsequently, based on these experimental data, the DeepBlow and ADMS/CDOG models were successfully verified(Chen and Yapa, 2003; Johansen et al., 2003). The Gulf of Mexico, as a typical deepwater region for oil/gas exploration(water depths >300 m), has also become an area of practical interest for researchers to further test their models using scenarios of hypothetical oil spills(Chen and Yapa, 2004; Dasanayaka and Yapa, 2009). More recently, after the 2010 Deepwater Horizon oil spill, which originated from a water depth of 1 500 m(Boufadel et al., 2014), the Gulf of Mexico has become a focus for deepwater oil spill research(Socolofsky et al., 2011; Yapa et al., 2012). In recent years, some researchers have begun to study numerical simulations of underwater oil spills in the off shore area of China. Many of these have focused on the shallow Bohai Sea(maximum depth 85 m)as their region of interest(Wang et al., 2008, 2013; Guo and Wang, 2009, Li et al., 2013), probably because current off shore oil exploitation and production by China is mainly concentrated within this area. However, little research has been published on how oil might behave if released accidentally in the SCS, especially in deepwater regions far from the mainl and. In this study, an off shore oil/gas field in the northern SCS was selected as the location for a hypothetical oil spill simulation.

According to Dasanayaka and Yapa(2009), for contingency planning and emergency response, an oil spill model must address two issues: the time it takes for the oil to appear at the surface and its approximate location. For this, the description of the underwater oil spill process must be as precise as possible. Based on the analysis of the underwater transport of oil released from deepwater and on previous works, we know that the underwater process of an oil spill is dominated principally by three factors: initial jet momentum, plume buoyancy, and ambient current and turbulence. Accordingly, in the present study, the entire underwater process of an oil spill was separated into three successive stages: the turbulent jet stage, buoyant plume stage, and advection-diff usion stage(Fig. 1). The turbulent jet stage exists when a mixture of oil and gas is released from a violent underwater blowout where the velocity at the orifice can reach 5-10 m/s. During this stage, the movement of the jet is dominated by the initial momentum, whereby the oil breaks into a large number of droplets of unequal sizes due to the high turbulence. Because ambient water is also entrained into the jet, a rapid loss of momentum occurs within a few meters from the release location, then the turbulent jet stage terminates. In the buoyant plume stage, the momentum is no longer significant relative to the buoyancy, which then becomes the driving force behind the movement for the remainder of the oil plume. The buoyancy eff ect causes the plume to keep rising toward the sea surface while ambient water is continuously entrained into the plume. As a result, the plume becomes enlarged and its density steadily approaches the ambient density. When the plume becomes fully developed, a considerable amount of water containing the oil droplets is pumped to the shallower region. Because ambient seawater is usually much denser in the deeper water, the oil plume may reach a level of neutral buoyancy at some depth. At this point, buoyancy no longer dominates the movement of the plume and the buoyant plume stage terminates. Subsequently, the plume dynamics become negligible and the oil continues to move as individual oil droplets passively following the ambient current and turbulence and rising due to droplet buoyancy. This is the so-called advection-diff usion stage. It should be noted that there is no obvious boundary between two successive stages. Actually, the three factors mentioned above exist in every stage; it is simply that the dominant factor varies in the diff erent stages. In the present study, oil transport during the three stages was simulated using two submodels. The first two stages were simulated with a plume dynamics model because the oil in both stages can be treated as an entirety(oil jet/plume), whereas the third stage was modeled with an advection-diff usion model. Many existing underwater oil spill models use both submodels(e.g., DeepBlow, ADMS/CDOG). There are some models(Lonin, 1999; Wang et al., 2013), however, that are based solely on the advectiondiff usion model and the possibility exists that they might be used in critical decision making and environmental impact assessment.

Fig. 1 Sketch of underwater oil spill process

In this study, a numerical model was developed to simulate the underwater transport of oil from a deepwater spill. Based on such simulations and given the spillage information and hydrodynamic conditions, it becomes possible to predict how the oil will be distributed within the water column, and when and where the oil will first appear at the surface. Such simulation results would be useful in contingency planning for oil spill emergency response.

The remainder of this article is organized as follows. We start by showing the formulation of the underwater oil spill model in Section 2. In Section 3, the verification of the model is demonstrated by a series of numerical tests in which the model predictions are compared with experimental data recorded in previous literature. Section 4 discusses the application of the model to a scenario in which a hypothetical oil spill occurs at the seabed of an oil/gas field in the northern SCS. Finally, we provide a summary and draw some conclusions in Section 5.

2 MODEL FORMULATION

The oil spill model presented in this paper comprises two submodels: the plume dynamics model(PDM) and advection-diff usion model(ADM). The PDM is used to simulate both the turbulent jet and buoyant plume stages, during which the mixture of water and a certain amount of spilled oil is treated as an entirety and the interaction between the oil and ambient water is considered. The remaining advection-diff usion stage is simulated by the ADM, where the spilled oil is divided into a large number of discrete particles. In the ADM, each particle represents a set of oil droplets of equal size, characterized by its spatial coordinates, velocity, volume, oil concentration, and droplet size(diameter). These particles are introduced into the water environment where the plume terminates(hereafter, the jet and plume are referred to simply as the plume), then they move in response to the shear current, turbulence, and buoyancy.

2.1 Plume dynamics model(PDM)

In the PDM, the Lagrangian integral technique is used to simulate the turbulent jet and buoyant plume stages of the oil spill. As shown in Fig. 2, the oil spill duration is divided into a series of intervals of equal length Δ t and each time interval corresponds to a small amount of spilled oil. Therefore, the oil/water plume can be represented as a series of non-interfering control elements. At a certain moment, the line connecting the center points of all the control elements is regarded as the oil trajectory. In this model, it is assumed that each control element is a cylindrical section of a bent cone with its bottom plane perpendicular to the plume trajectory. The local velocity(m/s), radius(m), and density(kg/m3)of the element are V, b, and ρ, respectively. Thus, the element thickness(m)is h = |Vt and its mass(kg)is m = ρ π b 2 h. The element is further characterized by a set of variables such as temperature, salinity, and oil concentration. These variables represent the average values for an element and they will change as the element moves within the 3D space. The following governing equations are applied to the control element.

Fig. 2 Sketch of control volume of plume dynamics model(Lee and Cheung, 1990)

(1)Conservation of mass

The control element may change in mass due to entrainment, oil dissolution, and turbulent diff usion, while rising because of buoyancy and becoming sheared over the ambient flow. The change of the element’s mass can be described by the following:

where ρ a is the density of the ambient fluid(kg/m3) and Q e is the volume flux(m3 /s)entrained through both shear-induced and forced entrainment. The second term of the right-h and side of Eq.1 represents the mass loss rate due to oil dissolution from the buoyant plume into the ambient environment and it can be computed based on the method of Cohen et al.(1980)as
where m i is the mass loss(kg)due to dissolution of the oil component i, K r is the dissolution mass transfer coefficient(0.01 m/hr), f S is an empirical constant(0.7, Rye, 1994), A =2π bh is the area of exposed plume surface, and Si is the freshwater solubility(kg/m3)of component i. Huang and Monastero(1982)suggested that for typical oil, the solubility could be calculated by
where S 0 is the solubility of fresh oil(0.03 kg/m3), α is a decay constant(0.1), and t is the time after spillage(hr).

The third term of the right-h and side of Eq.1 represents the mass loss rate due to turbulent diff usion, which can be related to the concentration gradient and written as

where m dif is the oil mass loss(kg)due to turbulent diff usion, K C is the oil concentration diff usivity(m2 /s), | ə C / ə r | is the norm of oil concentration gradient that can be approximated by(C - Ca)/ b, and Cand Ca are the oil concentrations(mass fraction of oil content)in the control element and ambient environment, respectively.

(2)Conservation of momentum

Once the oil is spilled into the ambient environment, the oil momentum will change because of interaction between the oil and ambient fluid. According to the law of conservation of momentum, the rate of change of momentum of a control element should satisfy the following equation:

where a V is the average velocity vector(m/s)of the ambient fluid flowing over the front surface of the plume, m i is the mass loss(kg)due to oil dissolution and turbulent diff usion, Δ ρ = ρ a - ρ is the density diff erence(kg/m3), g is the gravitational acceleration(9.81 m/s 2), C D is a drag coefficient, a a is the projection of a in the direction of, and k is the unit vector in the vertical direction. The four terms of the right-h and side of Eq.2 represent the momentum change due to entrainment, the buoyant force, oil dissolution and turbulent diff usion, and the drag force, respectively. Compared with the former two terms, the latter two have much smaller order of magnitude and therefore, they can be neglected during the simulation. Thus, Eq.2 is reduced to a simpler form:

(3)Conservation of heat, salinity, and oil mass

The buoyant plume will also exchange material and heat with the ambient environment. This leads to changes in the heat, salinity, and oil concentration of a control element, which can be described by the following conservation equation:

where the symbol I is a scalar parameter and denotes the heat C p T, salinity S, or oil concentration C ; C p is the specific heat(for water C p =4 216.3 J/(kg·K) and for crude oil C p =2 200 J/(kg·K); Jamaluddin et al.(1991)); T is temperature(K); and K represents the diff usivities. For heat, salinity, and oil mass, the term F I is written as F T =d [ m dis+dif(C p T)oil ]/ dt, F S = d(m dis+difs oil)/ dt, respectively. The diff usivities of heat, salinity, and oil concentration are represented by K T, K S, and K C, respectively. Equation 3 states that the changes of heat, salinity, and oil concentration of a control element are mainly due to the input contributed by the entrained liquid and the loss caused by turbulent diff usion and oil dissolution. In the present study, the eff ects of both turbulent diff usion and oil dissolution on the changes of heat and salinity are assumed insignificant and neglected in the simulation, i.e., K T = K S =0 and F T = F S =0.

(4)State equation

In the present model, the fluid within a control element is treated as a non-miscible mixture of seawater and oil. The size of the control element keeps growing via entrainment of seawater from the ambient environment, leading to changes to its temperature, salinity, and oil concentration and consequently, leading to the change to its density. The change of control element density in response to variation of its temperature, salinity, and oil concentration can be described by a state equation, which ultimately completes the governing equations. The functional form of seawater density can be found in several references(Pierson and Neumann, 1966). Bobra and Chung(1986)gave the functional form of density variation for many oils. According to the variation of temperature and salinity, Bemporad(1994)used a functional form of the state equation to calculate the jet density in a mathematical model for the numerical simulation of a submerged buoyant round jet. The functional form of the state equation presented by Bemporad(1994)is employed in the present PDM and it is written as

where ρ0 is a reference density(kg/m3), β T and β S are the thermal expansion coefficient(5×10 -4 /°C) and solutal expansion coefficient(8×10 -4), respectively, and Δ T(°C) and Δ S are variations from reference property values.

(5)Entrainment formulation

During migration of the buoyant plume through the water column, ambient fluid is entrained into the buoyant plume through the plume’s exterior surface and consequently, the plume mass is increased. A proper estimation of the rate of entrainment of the ambient fluid(Qe in Eq.1)plays an important role in the simulation of the buoyant plume. According to Frick(1984) and Lee and Cheung(1990), the entrainment can be separated into two distinct processes: shear-induced entrainment and forced entrainment. The former is due to shear between the buoyant plume and ambient fluid and it is present even when no ambient flow exists. The latter is due to the advection of ambient flow into the buoyant plume. In regions close to the discharge point or in very weak flow, shear-induced entrainment dominates. In general, however, forced entrainment dominates except very close to the source.

Accordingly, in the present model, the entrainment rate Qe is calculated as follows:

Qe = Qs + Qf,

where Qs and Qf represent the volume fluxes(m3 /s)due to the shear-induced and forced entrainment, respectively. The volume flux due to shear-induced entrainment is calculated with the functional form proposed by Fischer et al.(1979):

where α is an tunable entrainment coefficient. It is important to determine the value of the entrainment coefficient. From the kinetic energy and momentum equations, Schatzmann(1979, 1981)derived a hypothetical form of shear-induced entrainment as a function of the local densimetric Froude number and jet orientation. Subsequently, by comparing this hypothetical form with experimental data, Lee and Cheung(1990)presented a formulation to calculate the shear-induced entrainment coefficient, written as where ϕ is an angle between the jet trajectory and horizontal plane, and E is a constant of proportionality. By comparing the results of a numerical model with the asymptotic solutions of basic jet flow, Zheng and Yapa(1998)found that a value of E =2 was appropriate. Equations 7 and 8 provide the correct shear-induced entrainment coefficient for a pure source of a buoyant plume(Lee and Cheung, 1990).

If ambient flow exists, the ambient fluid will advect into the buoyant plume, causing forced entrainment. The forced entrainment of ambient fluid into the control element occurs on the upstream side of the plume as the control element moves along the centerline of the buoyant plume. The volume flux Qf caused by forced entrainment is formulated as

where Qfx, Qfy, and Qfz are the forced entrainment components in the x, y, and z directions, respectively. According to Frick(1984), forced entrainment consists of three terms. The first term represents forced entrainment due to the projected plume area normal to the ambient flow(cylinder term); the second is a correction because of the growth of the plume radius(growth term); and the third is a correction because of the curvature of the trajectory(curvature term). Based on the formulations of Frick(1984) and Lee and Cheung(1990), Yapa and Zheng(1997)derived an expression for forced entrainment that covers the wider range of flow conditions presented in the ocean; i.e., a 3D trajectory in a 3D ambient flow field. The derivation was presented in Yapa and Zheng(1997) and the results are given as follows:

where u a, v a, and w a are components of V a in the x, y, and z directions, respectively; θ is the angle between the x -axis and the projection of the plume trajectory in the horizontal plane(Fig. 2); is the displacement of a control element during one time step; and Δ x, Δ y, and Δ z are components of the displacement vector in the x, y, and z directions, respectively. In Eqs.10-12, the three terms in the brackets of the right-h and side represent the cylinder term, growth term, and curvature term, respectively. It is worth noting that Eqs.10-12 assume that all the ambient flow on the upstream side of the plume is entrained into the control element.

2.2 Advection-diff usion model(ADM)

When the oil plume travels upwards to a certain depth level, its density can be very close to that of the ambient water because a large amount of water will have been entrained into the plume. In a densitystratified water environment, the ambient seawater is denser in the deeper regions than nearer the surface. Denser seawater is entrained into the plume in the deeper regions, and as the oil plume rises to the depth where the seawater is less dense, the fluid of the plume may reach a neutral or even a negative buoyancy level below the sea surface. As a result, in the framework of the PDM, the control element may lose its upward driving force and attain a final rise height. Obviously, this is inconsistent with actual cases of underwater oil spill accidents, where the spilled oil can usually be observed at the sea surface. In fact, as the oil plume will have developed to a certain extent, the fluid property and turbulence intensity inside the plume may already have approached those outside.

Therefore, instead of following the mechanism of plume dynamics, the advection-diff usion mechanism of the ambient flow will dominate the behavior of the spilled oil and the buoyant PDM will no longer be applicable. Consequently, the ADM will be employed to simulate the transport of the spilled oil. The governing equation of the advection-diff usion mechanism can be written as

where C is the oil concentration(mass fraction of oil content)in the water column, V is the advection velocity vector(m/s),  is the gradient operator, K =(Kx, Ky, Kz)is the turbulent diff usion tensor in water(m2 /s), and Kx, Ky, and Kz are the diff usion coefficients in the x, y, and z directions, respectively. Generally, water turbulence can be treated as homogenously isotropic in the horizontal directions and thus, both Kx and Ky can be identified as a horizontal diff usion coefficient Kh. The sink/source for the ADM is represented by the term Si in Eq.13. In the present study, the final position of the oil plume obtained by the PDM serves as the source position for the ADM.

The advection-diff usion equation(Eq.13)can be solved using a variety of numerical methods based on the Eulerian concept, including finite diff erences and finite elements(Gomez et al., 2011; Wang et al., 2013). One difficulty with these methods is the stability problem; a numerical diff usion, which is unrelated to the physical diff usion process, may be generated during numerical discretization and a large numerical diff usion might introduce distortion to the model result. Additionally, these Eulerian methods could cause erroneous drifts in particle trajectories for non-uniform flow fields(Bennett and Clites, 1987) and they have difficulty in describing some processes, such as the non-Fickian transport in heterogeneous geological media(Berkowitz et al., 2006). In contrast, the Lagrangian particle-tracking(r and om walk)method is more commonly used in the simulation of oil transport within a water environment. The Lagrangian particle-tracking method was adopted very early on for use in ocean and atmosphere research(Csanady, 1973). Johansen(1984) and Elliott(1986)were the first to apply this method to oil spill simulations. This method not only avoids the numerical diff usion problem that besets the discretized advection-diff usion equation, but it also accurately describes the physical diff usion process including the breakup and separation of an oil slick. Numerous studies have shown that this technique is significantly more efficient in cases with higher dimensions, where there are relatively fewer particles and the contaminant resides in a patch that occupies only a portion of the total model area(Hunter, 1987; Al-Rabeh et al., 1989; Yapa et al., 1999; Wang et al., 2008; Dasanayaka and Yapa, 2009).

In the Lagrangian particle-tracking algorithm, the spilled oil is divided into a large number of discrete particles. Each particle contains a large number of oil droplets of equal size(diameter) and it is assigned with a set of spatial coordinates. These particles are introduced into the water environment at the source site at a rate corresponding to the spillage rate. It is assumed these particles advect with the surrounding water body, in response to dynamical sea processes such as currents, waves, and buoyancy, and diff use because of r and om processes due to turbulence. The advection velocity can be obtained by the deterministic method, while the diff usion velocity is represented by a 3D r and om walk technique.

In the present model, the Lagrangian particletracking algorithm is used to simulate oil transport in the advection-diff usion stage. The 3D movement of an oil particle is determined by the following relationship:

where S =(x, y, z)is the displacement vector of an oil particle; x, y, and z are the Cartesian coordinates; V is the advection velocity due to the combined eff ects of the current and waves within the water column; V is the diff usion velocity due to the turbulent diff usion process; wb is the buoyancy velocity of the oil particles; and k is the unit vector in the vertical direction. The advection velocity V is a significant factor in the transport of oil and it can be obtained from a 3D hydrodynamic model(e.g., POM, ROMS). The diff usion velocity V is a r and om variable, which can be calculated using the following formula: where Rx, Ry, and Rz are independent of each other and assumed uniformly distributed r and om numbers ranging from -1 to 1.

The buoyancy velocity wb contributes part of the vertical velocity of an oil particle and it is determined by the particle size, seawater viscosity, and density diff erence between the seawater and oil particle. According to Yapa et al.(1999), buoyant velocity can be calculated using either Stokes’ law or Reynolds’ law, depending on the oil droplet size. First, a critical droplet diameter d c is defined by

where μ is the seawater viscosity(about 10 -3 Ns/m2), and ρ and ρ a are the oil and seawater densities, respectively. For small particles where d < d c, wb is computed based on Stokes’ law and it is given by whereas for large particles where d ≥ d c, wb is computed based on Reynolds’ law and it is given by

From Eqs.17 and 18, it is evident that the oil droplet size is an important factor in the ADM. In deepwater oil spill simulations, the size distribution of oil droplets might have an important eff ect on the oil transport and concentration and , to a large extent, could determine the surface appearance of the released oil(slick formation)(Johansen, 2003). The oil droplet size distribution is greatly dependent on the conditions at the release point as well as on the oil properties(density, interfacial tension, viscosity). Unfortunately, few quantitative studies have been reported on droplet splitting in large volume oil jets emerging into water, and no reliable method for estimating the droplet size distribution in such systems exists(Masutani and Adams, 2001). However, in some underwater oil spill experiments, the oil droplet size distribution has been recorded and can be found in some literature(Johansen et al., 2001, 2013; Yapa et al., 2012; Brandvik et al., 2013). Meanwhile, some partially tested empirical or theoretical equations have also been developed for calculating droplet size(Delvigne and Sweeney, 1989; Delvigne, 1994; Chen and Yapa, 2007; Bandara and Yapa, 2011; Brandvik et al., 2013; Johansen et al., 2013; Zhao et al., 2014a, b). To some extent, both the data and the equations are applicable to practical oil spill simulations.

In the present model, the oil particles are introduced into the water environment at the end of the buoyant plume stage, when the transition from the plume dynamics stage to the advection-diff usion stage takes place. At the transition point, all particles in a control element are presumed uniformly distributed and their current positions are taken as their initial positions in the advection-diff usion stage. It is worth noting that the definition of this transition point is very important in this model. In the case where the water environment is stratified, Yapa et al.(1999)used a neutral buoyancy level as the transition point, i.e., the oil plume is considered to reach a terminal level once the plume density is equal to the density of the ambient seawater. Dasanayaka and Yapa(2009)examined three diff erent criteria that can be used as the transition point(in their paper the transition point was referred to as the Terminal Level of Plume Dynamics): the neutral buoyancy level criterion, droplet buoyant velocity criterion, and zero velocity criterion. Their parametric study suggested that the transition point could be based on the buoyant velocity corresponding to the median oil droplet size. Nevertheless, in this study, considering the stratification of the seawater, we followed Yapa et al.(1999)in using the neutral buoyancy level criterion to define the transition point between the PDM and ADM.

3 MODEL VERIFICATION

To examine the validity of the oil spill model, the model was applied to several cases of underwater buoyant jets and the model predictions compared with experimental data. The experimental data include seven sets of laboratory experimental data and one set of field experimental data that have been reported in earlier literature.

3.1 Comparison with laboratory experimental data

This section reports how the laboratory experimental conditions were parameterized and then substituted into the numerical model to perform the simulations. The comparison between the simulated and observed jet trajectories is presented here. There were seven tests corresponding to the seven buoyant jet cases, considering diff erent jet conditions and diff erent ambient conditions. To distinguish the characteristics of the various buoyant jets, three parameters are introduced:

Densimetric Froude Number

Jet/Crossflow Ratio

Stratification Number

Here, the subscript “0” denotes the initial value, g ′ = g(Δ ρ / ρ), Δ ρ = ρ a - ρ, and D is the diameter of the jet orifice(m). The values of the parameters used in seven experiments as well as in the seven numerical tests are listed in Table 1. According to the experimental setting described in the previous literature, the ambient flow is unidirectional and uniform. For simplicity and without loss of generality, the direction of the ambient flow in the numerical simulation is assumed to be in the x -axis, i.e., a V =(u a, 0, 0).

Table 1 Laboratory experimental parameters

Tests 1, 2, and 3 simulated three cases of vertically discharged buoyant jets into an unstratified stagnant flow. In Fig. 3, the model predictions of the jet trajectories are compared with the experimental data provided by Hirst(1972). In the three cases, the numerical results are in good agreement with the experimental data. As can be seen from Fig. 3, initially, the jets travel vertically. After some distance, the jets become significantly deflected downstream by the crossflow, as shown by the bent trajectory. A larger value of R 0 will lead to stronger jet deflection. In fact, the behavior of the jet near the orifice is dominated by its upward initial momentum. Although the buoyant force can accelerate the jet, a large amount of horizontally flowing water is also entrained into the jet. Consequently, a rapid loss of vertical momentum occurs and the upward movement of the jet slows. Meanwhile, the jet also obtains horizontal momentum from the entrained fluid, increasing its horizontal velocity. Thus, the combination of the decreasing upward velocity and increasing horizontal velocity generates the bent jet trajectory. As the jet grows, the momentum is no longer significant compared with the buoyancy, which becomes the principal driving force for the remainder of the jet. When the jet becomes fully developed, there is no longer significant variation in the horizontal and vertical velocities, then the jet trajectory approaches an oblique line. This interpretation appears in good qualitative agreement with the prediction shown in Fig. 3.

Fig. 3 Trajectories for oil spilled in an unstratified flowing environment(data from Hirst(1972))

Results of the other four tests, including 2D simulations in a stratified flow(Tests 4 and 5) and 3D simulations in an unstratified flow(Tests 6 and 7), are presented in the Appendix. Overall, good agreement between the model results and the data were obtained in all these cases, further supporting the validity of the PDM.

3.2 Comparison with field experimental data

It is worth noting that the above seven tests with which the validity of the PDM submodel PDM has been verified were based on small-scale laboratory experiments. To complete the model verification, in this section, another set of data from a larger-scale field experiment is employed to calibrate the other submodel ADM. The data used here are from a field experiment performed off the coast of Norway in 1995 by IKU, in conjunction with Norwegian Clean Seas(NOFO)in Norway, Esso Norge, and Norsk Hydro, with the aim of simulating oil release from a ruptured pipeline. Rye et al.(1996) and Rye and Brand vik(1997) reported this experiment in detail. Based on their reports, the experimental parametersare summarized in Table 2, recorded values of water salinity, temperature, and corresponding water density are plotted in Fig. 4, and recorded values of the sea currents are listed in Table 3. Based on sonar recordings, the width and position of the underwater oil plume, as well as the cloud of oil droplets projected to a vertical plane in the SE-NW direction, are shown by solid lines in Fig. 5. As can be seen, the spilled oil behaved as an oil plume in the deeper parts. After having risen closer to the sea surface, the oil began to be transported passively as a cloud of individual oil droplets under the combined eff ects of buoyancy and the sea currents. According to the video recording analysis in Rye et al.(1996), the transition zone was at a depth of 50-60 m, and the first oil droplets were observed at the sea surface 10 min after the oil release commenced.

Table 2 Field experimental parameters(Rye et al., 1996)

Fig. 4 Vertical profiles of salinity, temperature, and potential density in field experiment(Rye et al., 1996)

Fig. 5 Comparison between model results and field observations from Rye et al.(1996) a. projection of underwater oil to a vertical plane in the SE-NW direction. Dashed curve denotes the envelope of all oil particles and solid lines represent the field observations recorded at diff erent times; b. calculated and observed width and position.

The data reported in Tables 2 and 3, as well as those shown in Fig. 4, were substituted into the present model as input conditions to simulate the field experiment. The current velocities were calculated through linear interpolation of the data in Table 3. As mentioned in Section 2.2, droplet size distribution plays an important role in the ADM model. In the present simulation, the oil droplet diameter was assumed normally distributed within a range of 6-5 000 μm, with a mean of 250 μm and st and ard deviation of 75 μm, according to Delvigne and Sweeney(1989) and Delvigne(1994). In addition, the oil plume may reach a neutral buoyancy level below the water surface because the ambient seawater was stratified. Thus, the neutral buoyancy level criterion was employed to define the transition point between the two submodels, i.e., once the density of the plume equaled that of the ambient seawater, the transport of the spilled oil was simulated using the ADM submodel instead of the PDM. In the ADM, the horizontal and vertical diff usivities were 0.05 and 0.001 m2 /s at the sea surface, respectively. It was assumed that the diff usivities were constant in the top 10 m of the water column, but decreased with depth at a rate proportional to the distance between the oil particles and the water surface(Yapa et al., 1999).

Table 3 Current velocity profile of field experiment(Rye et al., 1996)

The model predictions of the width and position of the oil plume, and of the envelope of the droplet cloud projected to a vertical plane in the SE-NW direction, are compared with the field data in Fig. 5a. The lower part of the results(depth range 50-107 m)corresponds to the envelope predicted by the PDM and the upper part(depth range 0-50 m)corresponds to the cloud of oil particles predicted by the ADM. According to the simulation results, the transition from the lower part to the upper part occurs at the depth of 54.6 m, which is in good agreement with the field observations, i.e., the oil plume is trapped at this depth. Above this depth, the spilled oil appears to rise as individual droplets. The larger droplets rise more quickly than the smaller ones. Therefore, the larger droplets are expected to remain underwater for a shorter length of time and thus, the crossflow transports them over shorter distances than the smaller droplets. According to the model predictions, the largest droplets(with diameters of 5 mm)reach the sea surface first, and this happens 11.3 min after the oil release commenced, which is comparable with the field observations. The calculated and observed widths, as well as the center positions, of the spilled oil are compared in Fig. 5b. As can be seen, the overall agreement is satisfactory.

However, some discrepancies between the model predictions and field observations can be observed. On the one h and , in the bottom part, the model simulates the plume position well but it underpredicts the width, especially at the depth of 100 m. This discrepancy must be interpreted considering the vertical movement of the oil-release equipment caused by wave action on the boat during the experiment(Rye et al., 1996). This could have caused some extra initial mixing and thus, increased the entrainment and the consequent spatial extent of the buoyant jet. By ignoring this factor in the present model, a relatively small plume width is predicted in the bottom part. On the other h and , in the upper part, the calculated and observed widths of the oil droplet cloud are comparable, although a shift occurs between the calculated and observed center positions. This discrepancy can be attributed to insufficient detail regarding the unsteady current velocity profile. The advection velocity field, which plays a significant role in the ADM, is derived only from data measured at four depths and at two time points, whereas the oil transport occurs over a much larger space. Thus, the inaccuracy of the derived velocity field means that the simulated underwater distribution of the spilled oil deviates from the observations. Furthermore, the comparison is on a 2D plane, whereas the actual physical process occurs in 3D space. The observation plane might not be strictly in a SE-NW alignment and any slight shift of the plane might lead to somewhat diff erent results. In summary, the overall agreement is reasonably satisfactory. If the lack of information mentioned above were addressed, then this would enhance the rationality and validity of the model.

4 MODEL APPLICATION

To demonstrate the model capability further, the model was applied to a scenario of practical interest, i.e., a hypothetical oil spill at the seabed of an oil/gas field located in the north of the SCS(Fig. 6). The background hydrodynamic data, including the current velocity, temperature, salinity, and potential density, were provided by an operational 3D current forecasting system based on the σ-coordinate ROMS/ TOMS(owned by the First Institute of Oceanography, State Oceanic Administration). As the interest of the present study is focused on the underwater oil trajectory within the first few hours after release, wind data were not included. The time range of the hydrodynamic data was 48 hr(from12:00 on June 14 to 12:00 June 16, 2013), and the interval was 1 hr. The model domain was sufficiently large such that the spilled oil would not escape it. The horizontal resolution was 1/8°×1/8°. In the vertical direction, the water column was divided equally into 20 σ layers. The oil/gas field of interest is located at 19°55′N, 115°25′E, 300 km to the southeast of Hong Kong, where the local depth is 1 442 m. The vertical profiles of the local salinity, temperature, and potential density are illustrated in Fig. 7. As can be seen, the time variations in the three profiles are very small, especially at depths >500 m. Figure 8 presents the time-dependent vertical profile of current velocity. Note that the vectors in Fig. 8a have been normalized by their initial( t =0 hr)velocity magnitudes for easier visualization(Fig. 8b). Actually, considering the initial velocity profile(Fig. 8a), the currents in the shallower layers are much stronger than in the deeper layers. The oil is assumed spilled upwards from the seabed into the sea. The discharge conditions are summarized in Table 4. The horizontal and vertical diff usivities are the same as those described in Section 3.2. The oil spill model employs 100 000 Lagrangian particles to simulate the transport of the spilled oil. The oil droplet size distribution is the same as that used in Section 3.2. The series of insets in Fig. 9 show the simulated underwater particle locations projected to the W-E plane and S-N plane at four diff erent times. This helps the visualization of the spatiotemporal changes of the oil distribution.

Fig. 6 Bathymetry of the northern South China Sea and oil spill location denoted by the star

Fig. 7 Vertical profiles of salinity, temperature, and potential density In each panel, the three diff erent colored lines represent data at three diff erent times.

Fig. 8 Temporally varying vertical profile of current velocity at oil spill location a. vertical profile at t =0 hr. Blue, red, and black lines represent the eastern U, northern V, and velocity magnitude, respectively; b. time history of current velocity vectors, normalized by the magnitude at t =0 hr, in 20 layers. Time t =0 hr corresponds to 12:00 on June 14, 2013.

Table 4 Parameters used in the simulations

As can be seen from Fig. 9, the spilled oil initially behaves as a plume in the bottom part and then as a cloud of individual free droplets in the shallower part. The bottom part is based on the PDM. Note that no particle is actually employed in the PDM. The particles are plotted in the plume part only for easier visualization of the shape and position of the plume. In the plume part, the particles are not free; they are fixed uniformly and stochastically to a certain control element until that element reaches the terminal level. As the currents near the bottom are very weak(<2.5 cm/s, Fig. 8), the horizontal shift of the plume is insignificant compared with the large horizontal scale of the plots and therefore, the plume appears to rise almost vertically in Fig. 9. The oil plume rises about 76 m(1 366-m depth)until the buoyant plume reaches the neutral buoyancy level and terminates. Subsequently, the oil particles break away from the plume element and they are transported passively under the combined eff ects of buoyancy and the currents, while spreading in all three directions because of turbulent diff usion.

Fig. 9 Distribution of simulated underwater particles projected to the W-E plane(a, c, e, g) and S-N(b, d, f, h)at four diff erent times The droplet size ranges are 6-175, 175-325, and 325-5 000 μm for the small, medium, and large particles, respectively.

In Fig. 9, the part above the plume stage corresponds to the simulation results from the ADM. To facilitate the analysis, all particles were divided into three classes according to the oil droplet size. These classes are plotted in Fig. 9 as red, blue, and black dots corresponding to large particles(droplet diameters >250+75 μm), medium particles(droplet diameters between 250−75 and 250+75 μm), and small particles(droplet diameters <250-75 μm), respectively. As the droplet diameters are assumed normally distributed with a mean of 250 μm and a st and ard deviation of 75 μm, the medium particles account for 68.3% of the total oil volume, while the remaining 31.7% is divided equally between the larger and smaller particles. As shown in Fig. 9, the rising trajectories of the particles are deflected laterally because of the crosscurrents. As the smaller droplets move toward the surface more slowly, their horizontal deflection is greatest.

Furthermore, the droplets of diff erent sizes might reach the surface at diff erent times and locations. The larger droplets reach the surface at an earlier time and at a location nearer to the release point. According to the simulation from the present model, the first particles reach the surface 141.3 min after the release commenced. The smaller droplets might remain underwater for a longer time and they could be carried a larger distance downstream. From Fig. 9, it can be seen that within the first 24 hours, only a fraction of the large particles reach the surface. After 48 hours, a small part of the large particles and part of the medium particles remain underwater, while all the small particles remain trapped in a nearly suspended state below the depth of 750 m. The horizontal distance traveled from the release location by the surfacing droplets ranges from several kilometers for large particles to over 10 km for medium particles. Thus, it can be deduced that the smallest droplets will reach the surface much further away.

Influenced by the time-varying crosscurrents, the underwater oil distribution exhibits diff erent patterns and diff erent scales at diff erent times. In the first 12 hours, as can be seen in Fig. 9a, the underwater oil above the depth of 500 m is mainly distributed to the west of the release location, while the oil in the deeper water is mainly distributed to the east. Over the remaining 36 hours, the underwater oil becomes mainly distributed to the east of the release location and the pollution area extends further eastward. However, since northward currents are dominant in this region(Fig. 8), the underwater oil is mainly distributed to the north of the release location and the pollution area extends steadily northward over most of the 48 hours, as shown by Fig. 9b, d, f, h. This also leads to a pattern in which the extension scale in the S-N direction is increasingly larger than in the W-E direction.

Figs. 10 and 11 show the oil budget. Figure 10 shows the history of the volume fractions of oil reaching the sea surface and oil remaining underwater. As can be seen, the volume fraction of underwater oil increases steadily during the release in the first 24 hours, because oil was introduced continually into the sea. As the oil requires at least 141.3 min to reach the sea surface, the growth of the volume fraction of the surface oil lags behind that of the underwater oil for the same time length. As a result, the underwater oil volume fraction grows slightly faster in the first 141.3 min. At t =24 hr, when no further oil is released into the sea, the underwater oil volume fraction reaches its maximum at 52.5%, then it begins to decrease. After t =24 hr, the larger droplets continue to reach the surface until about t =27 hr, when only the remaining smaller droplets move towards the surface at a much slower speed. Consequently, both the surface and underwater oil volume fractions exhibit slower variation from 27 to 48 hr. At the end of 2 days, the surface oil volume fraction has increased to 60.6%, while the underwater oil volume fraction has decreased to 39.4%.

To examine the vertical distribution of the underwater oil, a statistical analysis was performed on the oil volume in three layers of the water column: 0-500, 500-1 000, and 1 000-1 442-m depth. The history of the volume fraction of underwater oil in the three layers is shown in Fig. 11. As can be seen, the variation of the oil volume fraction of the layers is similar to that of the total underwater oil volume fraction shown in Fig. 10, although the range of variation is smaller. As the spilled oil intrudes into the three layers at diff erent times, the growth of the oil volume fraction in the three layers also begins at diff erent times. For each curve in Fig. 11, the rapid growth at the beginning can be attributed to a quick intrusion of the largest droplets into the layer. Subsequently, the growth of the oil volume fraction of the layer slows, because the input rate of the largest droplets in each layer becomes almost equal to the output rate of the largest droplets. The continued growth of the oil volume fraction of the layer is attributable to the slow intrusion of the smaller droplets(because of their smaller buoyancy velocity). At the end of the first day, the oil volume fractions in the three layers are 4.4%, 6.4%, and 41.7%, respectively, which account for 8.4%, 12.2%, and 79.4% of the total underwater oil volume, respectively. After 24 hr, there is no new oil injected into the system and within a short time, the oil volume fraction of the layer drops suddenly, which can be interpreted as the final part of the largest droplets quickly floating upward out of the layer. Subsequently, because of the continuous movement of the residual underwater oil droplets toward the surface, the oil volume fraction in the third layer gradually decreases, whereas that of the other two layers remains almost constant. At the end of the second day, the oil volume fractions in the three layers are 2.2%, 4.5%, and 32.6%, respectively, which account for 5.6%, 11.4%, and 82.7% of the total underwater oil volume, respectively. From Fig. 11, it can be clearly seen that throughout the 2 days of the simulation, the third layer contains a considerably higher amount of the underwater oil than the other two layers. Although the spilled oil is continuously floating toward the sea surface, after 2 days of the simulation, more than 80% of the underwater oil remains in the deepwater area(>1 000- m depth), which accounts for less than one third of the total water depth. Thus, it will require a much longer time before this part of the oil reaches the sea surface.

Fig. 10 History of volume fraction of oil reaching sea surface and remaining underwater Volume fraction 1.0 corresponds to 1 000 m3 of oil.

Fig. 11 History of volume fraction of oil contained in three layers of the water column From sea surface to seabed, the depth ranges of the three layers are 0-500, 500-1 000, and 1 000-1 442 m, respectively. Volume fraction 1.0 corresponds to 1 000 m3 of oil.
5 SUMMARY AND CONCLUSION

In this study, a numerical model was developed to simulate the underwater transport of oil from a deepwater spillage. In this model, the underwater transport process of spilled oil is divided into three stages: the turbulent jet stage, buoyant plume stage, and advection-diff usion stage and these three stages are simulated by two submodels: the plume dynamics model and advection-diff usion model. The former is based on the Lagrangian integral technique and is used to simulate the first two stages, while the latter is based on the Lagrangian particle-tracking technique and is used to simulate the third stage. This model not only considers the 3D structure of the ambient currents and turbulence, but also the vertical profiles of temperature and salinity and the density stratification of the ambient environment. The oil transport can be simulated from the origin of the spill until the oil reaches the sea surface. Given the spillage information and hydrodynamic conditions, the simulations can predict when and where the oil will first appear at the surface and how the oil will be distributed within water column. Such simulation results would be useful in contingency planning for oil spillage emergency response.

Combined with an operational 3D current forecasting system, this model is capable of the quick and real-time prediction of the transport of oil accidentally released underwater and therefore, it is a useful tool for contingency planning during an accidental spill. The presented model also lays the foundation for further investigations into areas such as the formation mechanism of underwater oil plumes, oil transport after dispersant application, and deepwater blowouts of oil/gas mixtures. Already, some related works are underway or planned as future developments.

6 ACKNOWLEDGEMENT

We would like to deeply thank the editors and anonymous reviewers for their constructive suggestions whereby this work has been improved greatly. We also thank the First Institute of Oceanography, State Oceanic Administration for providing the hydrodynamic data used in the study.

References
Al-Rabeh A H, Cekirge H M, Gunay N. 1989. A stochastic simulation model of oil spill fate and transport. Appl .Math. Model ., 13 (6): 322-329.
Bandara U C, Yapa P D. 2011. Bubble sizes, breakup and coalescence in deepwater gas/oil plumes. J. Hy d raul .Eng ., 137 (7): 729-738.
Bemporad G A. 1994. Simulation of round buoyant jet in stratifi ed fl owing environment. J. Hy d raul. Eng ., 120 (5): 529-543.
Bennett J R, Clites A H. 1987. Accuracy of trajectory calculation in a fi nite-diff erence circulation model. J .Comput. Phys ., 68 (2): 272-282.
Berkowitz B, Cortis A, Dentz M, Scher H. 2006. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev. Geophys ., 44 (2):RG2003, http://dx.doi.org/10.1029/2005RG000178.
Bobra M A, Chung P T. 1986. A Catalogue of Oil Properties.Consultchem, Ottawa, Canada.
Boufadel M C, Abdollahi-Nasab A, Geng X L, Galt J, Torlapati J. 2014. Simulation of the Landfall of the deepwater horizon oil on the shorelines of the gulf of Mexico.Environ. Sci. Technol ., 48 (16): 9 496-9 505.
Brandvik P J, Johansen Ø, Leirvik F, Farooq U, Daling P S. 2013. Droplet breakup in subsurface oil releases—Part 1: experimental study of droplet breakup and eff ectiveness of dispersant injection. Mar. Pollut. Bull ., 73 (1): 319-326.
Chen F H, Yapa P D. 2003. A model for simulating deep water oil and gas blowouts—Part II: comparison of numerical simulations with “Deepspill” fi eld experiments. J .Hydraul. Res ., 41 (4): 353-365.
Chen F H, Yapa P D. 2004. Modeling gas separation from a bent deepwater oil and gas jet/plume. J. Marine Syst ., 45 (3-4): 189-203.
Chen F H, Yapa P D. 2007. Estimating the oil droplet size distributions in deepwater oil spills. J. Hygraul. Eng ., 133 (2): 197-207.
Cohen Y, Mackay D, Shiu W Y. 1980. Mass transfer rates between oil slicks and water. Can. J. Chem. Eng ., 58 (5): 569-575.
Csanady G T. 1973. Turbulent Diff usion in the Environment.Springer Science, Dordrecht, Boston.
Dasanayaka L K, Yapa P D. 2009. Role of plume dynamics phase in a deepwater oil and gas release model. J. Hydro-Environ. Res ., 2 (4): 243-253.
Delvigne G A L, Sweeney C E. 1989. Natural dispersion of oil.Oil Chem. Pollut ., 4 (4): 281-310.
Delvigne G A L. 1994. Natural and chemical dispersion of oil.J. Adv. Mar. Tech. Conf ., 11 : 23-40.
Doneker R L, Jirka G H. 1990. Expert system for hydrodynamic mixing zone analysis of conventional and toxic submerged single port discharges (CORMIX1). Technical Report EPA/600/3-90/012, U.S. EPA Environmental Research Laboratory, Athens, Ga.
Elliott A J. 1986. Shear diff usion and the spread of oil in the surface layers of the North Sea. D eu tsch. Hydrogr. Z ., 39 (3): 113-137.
Fanneløp T K, Sjøen K. 1980. Hydrodynamics of underwater blowouts. Norweg. Marit. Res ., 4 : 7-33.
Fischer H B, List J E, Koh R C, Imberger J, Brooks N H. 1979.Mixing in Inland and Coastal Waters. Academic Press,New York.
Frick W E. 1984. Non-empirical closure of the plume equations. Atmos. Environ ., 18 (4): 653-662.
Gomez S, Ivorra B, Ramos A M. 2011. Optimization of a pumping ship trajectory to clean oil contamination in the open sea. Math. Comput. Model ., 54 (1-2): 477-489.
Guo W J, Wang Y X. 2009. A numerical oil spill model based on a hybrid method. Mar. Pollut. Bull ., 58 (5): 726-734.
Hirst E. 1971. Buoyant jets discharged to quiescent stratifi ed ambient. J. Geophys. Res ., 76 (30): 7 375-7 384.
Hirst E. 1972. Buoyant jets with three-dimensional trajectories.J. Hydraul. Div ., 98 (11): 1 999-2 014.
Huang J C, Monastero F C. 1982. Review of the state-of-theart of oil spill simulation models. Final Report Submitted to the American Petroleum Institute, Raytheon Ocean Systems Company, East Providence, RI.
Hunter J R. 1987. The application of Lagrangian particletracking technique to modelling of dispersant in the sea.In : Noye J ed. Numerical Modelling: Application to Marine Systems. North-Holland.
Jamaluddin A K M, Kalogerakis N, Bishnoi P R. 1991. Hydrate plugging problems in undersea natural gas pipelines under shutdown conditions. J. Petrol. Sci. Eng ., 5 (4): 323-335.
Johansen Ø, Brandvik P J, Farooq U. 2013. Droplet breakup in subsea oil releases—Part 2: predictions of droplet size distributions with and without injection of chemical dispersants. Mar. Pollut. Bull ., 73 (1): 327-335.
Johansen Ø, Rye H, Cooper C. 2003. DeepSpill-fi eld study of a simulated oil and gas blowout in deep water. Spill Sci .Technol. B ull ., 8 (5-6): 433-443.
Johansen Ø, Rye H, Melbye A, Jensen H, Serigstad B, Knutsen T. 2001. Deep Spill JIP Experimental Discharges of Gas and Oil at Helland Hansen, Parts I, II, and III-Technical report. SINTEF Applied Chemistry, Norway.
Johansen Ø. 1984. The Halten Bank experiment-observations and model studies of drift and fate of oil in the marine environment. In : Proceedings of the 11th Arctic Marine Oil Spill Program (AMOP) Technical Seminar.Environment Canada, Edmonton, Alberta, Canada. p.18- 36.
Johansen Ø. 2000. DeepBlow -a Lagrangian plume model for deep water blowouts. Spill Sci. Technol. B ull ., 6 (2): 103- 111.
Johansen Ø. 2003. Development and verifi cation of deepwater blowout models. Mar. Pollut. Bull ., 47 (9-12): 360- 368.
Lee J H W, Cheung V. 1990. Generalized Lagrangian model for buoyant jets in current. J. Environ. Eng ., 116 (6): 1 085-1 106.
Li Y, Zhu J, Wang H. 2013. The impact of diff erent vertical diff usion schemes in a three-dimensional oil spill model in the Bohai Sea. Adv. Atmos. Sci ., 30 (6): 1 569-1 586.
Lonin S A. 1999. Lagrangian model for oil spill diff usion at sea. Spill Sci. Technol. B ull ., 5 (5-6): 331-336.
Masutani M S, Adams E E. 2001. Experimental study of multiphase plumes with application to deep ocean oil spills.Technical report, Hawaii Natural Energy Institute,University of Hawaii, Honolulu.
Mcdougall T J. 1978. Bubble plumes in stratifi ed environments.J. Fluid Mech ., 85 (4): 655-672.
Pierson G, Neumann W J Jr. 1966. Principles of Physical Oceanography. Prentice-Hall, Inc., Englewood Cliff s N J.
Rye H, Brandvik P J, Reed M. 1996. Subsurface oil release fi eld experiment-observations and modelling of subsurface plume behavior. In : Proceedings of the 19th Arctic and Marine Oil Spill Program (AMOP) Technical Seminar, Vol. 2. Environment Canada, Ottawa. p.1 417- 1 435.
Rye H, Brandvik P J. 1997. Verifi cation of subsurface oil spill models. In : International Oil Spill Conference Proceedings: April 1997. Environment Canada, Ottawa, 1997 (1): 551-577.
Rye H. 1994. Model for calculation of underwater blow-out plume. In : Proceedings of the Seventeenth Arctic and Marine Oil Spill Program (AMOP), Technical Seminar.Environment Canada, Ottawa. p.849-865.
Schatzmann M. 1979. Calculation of submerged thermal plumes discharged into air and water fl ows. In : Proc. 18th IAHR Congress. Int. Assoc. for Hydr. Res., Delft, The Netherlands, 4 : 379-385.
Schatzmann M. 1981. Mathematical modelling of submerged discharges into coastal waters. In : Proc. 19th IAHR Congress. Int. Assoc. for Hydr. Res., Delft, The Netherlands, 3 : 239-246.
Shen H T, Yapa P D, Petroski M E. 1986. Simulation of oil slick transport in Great Lakes connecting channels.Department of Civil and Environmental Engineering,Clarkson University, Potsdam, NY.
Socolofsky S A, Adams E E, Sherwood C R. 2011. Formation dynamics of subsurface hydrocarbon intrusions following the deepwater Horizon blowout. Geophys. Res. Lett ., 38 (9): L09602, http://dx.doi.org/10.1029/2011GL047174.
Topham D R. 1975. Hydrodynamics of an oilwell blowout. In :Beaufort Sea Technical Report no.33. Department of the Environment, Victoria, B. C.
Wang C H, Li X Y, Lv X Q. 2013. Numerical study on initial fi eld of pollution in the Bohai Sea with an adjoint method.Math. Probl. Eng ., 2013, Article ID 104591, 10p.
Wang S D, Shen Y M, Guo Y K, Tang J. 2008. Threedimensional numerical simulation for transport of oil spills in seas. Ocean Eng ., 35 (5-6): 503-510.
Yapa P D, Wimalaratne M R, Dissanayake A L, DeGraff J A Jr. 2012. How does oil and gas behave when released in deepwater? J. Hydro-Environ. Res ., 6 (4): 275-285.
Yapa P D, Zheng L, Nakata K. 1999. Modeling underwater oil/ gas jets and plumes. J. Hygraul. Eng ., 125 (5): 481-491.
Yapa P D, Zheng L. 1997. Simulation of oil spills from under water accidents I: model development. J. Hydraul. Res ., 35 (5): 673-688.
Zhao L, Boufadel M C, Socolofsky S A, Adams E, King T, Lee K. 2014a. Evolution of droplets in subsea oil and gas blowouts: development and validation of the numerical model VDROP-J. Mar. Pollut. Bull ., 83 (1): 58-69.
Zhao L, Torlapati J, Boufadel M C, King T, Robinson B, Lee K. 2014b. VDROP: a comprehensive model for droplet formation of oils and gases in liquids -Incorporation of the interfacial tension and droplet viscosity. Chem. Eng. J ., 253 : 93-106.
Zheng L, Yapa P D, Chen F H. 2003. A model for simulating deepwater oil and gas blowouts-Part I: theory and model formulation. J. Hydraul. Res ., 41 (4): 339-351.
Zheng L, Yapa P D. 1998. Simulation of oil spills from underwater accidents II: model verifi cation. J. Hydraul .Res ., 36 (1): 117-134.