Underwater spreading and surface drifting of oil spilled from a submarine pipeline under the combined action of wave and current

Applied Ocean Research 64 (2017) 217–235

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Underwater spreading and surface drifting of oil spilled from a submarine pipeline under the combined action of wave and current Hongjun Zhu ∗ , Jiahui You ∗ , Honglei Zhao State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China

a r t i c l e

i n f o

Article history: Received 13 June 2016 Received in revised form 8 March 2017 Accepted 8 March 2017 Keywords: Oil spill Underwater spread Surface drift Wave and current VOF

a b s t r a c t Tremendous economic loss and environmental damages are caused by oil-spilling accidents in sea. Accurate prediction of the underwater spreading and surface drifting of oil spills is important for the emergency response. In the present study, numerical investigation on the underwater spread and surface drift of oil spilled from a submarine pipeline under the combined action of wave and current was carried out to examine the effects of physical ocean environment, leaking flux and spilled oil density and viscosity. Reynolds-Averaged-Navier-Stokes (RANS) equations, realizable k-ε turbulence model and volume of fluid (VOF) model are employed to describe the multiphase flow, and velocity-boundary wavemaking technique combined with the sponge layer damping absorber technique realizes the numerical wave flume. Oil spill experiments were conducted to validate the numerical model. The calculation results indicate that compared with the environmental conditions of still water, only current and only wave, a larger scope of underwater spreading and relatively slower rising rate and relatively faster drifting rate of oil droplets are observed under the combined action of wave and current. The leaking flux affects the floating time and dispersion concentration, while the ocean environment affects the horizontal migration and surface drifting. Under the specific conditions of present work, oil density has obvious effect on the underwater spread but limited effect on the surface drifting, while oil viscosity has little effect on both the two processes. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Worldwide exploitation of marine oil resources has been an effective way to solve the energy crisis caused by the depletion of onshore oil resources. However, offshore oil drilling and production activities and oil transportation by cargo ships and pipelines across the open and coastal seas increase the potential of oil spills. Both the marine ecosystems and environment would be damaged by oil pollution. Especially for the coastal areas, house large populations and important economic activities are vulnerable and sensitive regions [1]. Even though the accident may take place away from coast, the oil pollution could be transported onshore by wind drift, wave drift and currents. The recent oil spills tragedies in the Gulf of Mexico and Bohai Bay have brought us profound and painful lessons. Therefore, forecasting of oil spills and their trajectories is important for risk assessment and emergency response.

∗ Corresponding authors. E-mail addresses: [email protected], [email protected] (H. Zhu), [email protected] (J. You). http://dx.doi.org/10.1016/j.apor.2017.03.007 0141-1187/© 2017 Elsevier Ltd. All rights reserved.

A great many studies have illustrated the location and size of oil spills for a specific oil incident based on the satellite or radar observations [2–8]. Among the remote sensing technology, SAR (synthetic aperture radar) is the most commonly used instrument, which is not affected by clouds. Although the radar images can bring us a valuable information of oil spills detection, it just capture the real-time images of the surface oil spills. In order to track oil movement and predict its fate, some oil spill models have been proposed. As reviewed by Galt [9], Reed et al. [10] and Hackett et al. [11], commonly used operational oil spill models include OILMAP (oil spill model and response system), GNOME (general NOAA operational modelling environment), OSCAR (oil spill contingency and response), MOTHY (French operational oil spill drift forecast system), OSIS (oil spill identification system) and OD3D (oil drift 3 dimensional model). Their complexity has a wide range. Some are more dependent on environmental input parameters from atmospheric, ocean and wave forecast models or observations. The influence of oceanic and atmospheric physical variables and processes such as evaporation, dispersion, dissolution, emulsification, photo-oxidation and biodegradation, may change significantly the oil spill track [12].

218

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 1. Schematic of the computational domain and boundary conditions.

Table 1 Simulation cases. Case

Environmental condition

Oil density, ␳o (kg/m3 )

Oil viscosity, ␮o (Pa·s)

Oil leaking rate, uo (m/s)

Diameter of leak, d(m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Still water Current Wave Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current Wave and current

950 950 950 950 950 950 950 950 950 950 950 950 750 850 910 930 970 990 950 950 950 950 950 950

1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1 2 2.5 3 4.5 6

3 3 3 3 3 3 3 3 1 2 4 5 5 5 3 3 3 3 3 3 3 3 3 3

0.06 0.06 0.06 0.06 0.02 0.04 0.08 0.1 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06

A number of researchers have conducted effective attempts to simulate the process of oil spill by computational fluid dynamics (CFD) [13–18]. Chao et al. [13] developed a 2D oil spill model for surface oil slick to predict the major movement direction and polluted area. Wang et al. [14] proposed a two layers model consisting of surface oil slick and suspended oil droplets entrained over the depth of the flow based on Eulerian-Lagrangian approach. However, the oil is spilled from a surface source. Tkalich [15] used a consistent Eulerian approach to predict rate of oil entrainment by solving the layer-averaged Navier-Stokes equations. Sayol et al. [18] have proposed a lagrangian model for tracking surface spills. Li et al. [16] employed CFD software to simulate the oil spill drift-diffusion on the water surface with floating ices. Zhu et al. [17] have investigated the process of oil spilled from a submarine pipeline to free surface under shear current using CFD software. However, the actual marine environment is more complex. Wind and wave may also contribute to the spread of spilled oil as stated in Wang et al. [19]. For oil leaking from a damaged submarine pipeline, both the predictions of underwater spread and surface drift are important to determine a quick response. When and where the spilled oil can reach to the free surface and temporal variation in the horizontal migration distance of the spilled oil are key issues. Factors affecting

the spread and drift of oil spills include the physical characteristics of the spilled oil (density and viscosity), the physical ocean environment (wave, current and wind) and the leakage parameters (leaking point size and leaking pressure or velocity) etc. Few studies have comprehensively considered these factors. Since CFD method has been well applied in environment modelling [20–23] including oil spill issue, in this paper, CFD approach coupling with the volume of fluid (VOF) model is used to investigate the floating and drifting process of oil spilled from a submarine pipeline under the combined action of wave and current. The appearance of wind generated waves is common in the real ocean. It is assumed that the wave has been generated by the wind but the wind is stopped. The aim of this work is to examine the effects of the coupling of wave and current comparing with only wave and only current. The shallower the water depth, the shorter is the time required for oil floating, and the shorter is the time left to the emergency response. Therefore, the oil spill in shallow water is more challenging. In this work, the water depth is set to 18 m as Bohai Sea, which is a typical semi-enclosed shallow sea with average water depth of 18 m, significant wave height of 1 m and average wave length of 54 m [24,25]. Oil spill is usually affected by a many factors such as the climatology, the wave and current directions, the

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

219

Fig. 2. Schematics of the experimental setup: (a) experimental arrangement in the water flume; (b) high speed camera on one side of the flume (top view); (c) ADV (Acoustic Doppler Velocimeter); (d) experimental apparatus.

vertical density distribution of the ocean and the water depth etc. In this work, only Bohai average climatology was considered. The effects of marine environment, leaking point size and leaking rate, oil density and viscosity on the process of oil floating and drifting are discussed in detail. 2. Mathematical model and numerical approach 2.1. Governing equations In order to capture the surface wave and oil drift, three fluids including oil, water and air are considered in simulation, and they are treated as incompressible flow. At the interface of fluids,

no phase change is assumed. Because the vertical floating and the horizontal migration are the two concerns in this study and huge CPU cost is required for a 3D calculation, two-dimensional (2D) simulations are performed with x-axis representing the horizontal direction and z-axis representing the vertical direction. Spilled oil spreading directly downstream can reach the maximum horizontal migration distance. Therefore, 2D simulation is accurate enough to capture the maximum horizontal migration distance, and can give a relatively fast prediction that is welcome to the practice. Under the combined action of wave and current, the spilled oil is transported in a complex turbulent and time-dependent hydrodynamic field. Therefore, Reynolds-Averaged-Navier-Stokes (RANS) equations and realizable k-ε turbulence model are employed to

220

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 3. Comparison of the rising height of oil droplets and underwater spread at different leaking pressures between experimental results (ER, performed in a water flume with water depth of 18 cm) and simulation results (SR, contours of simulation results are presented in oil volume fraction): (a) Po = 99 hPa; (b) Po = 88 hPa; (c) Po = 75 hPa.

describe the multiphase flow. The RANS equations are expressed as [26,27]

∂ui =0 ∂xi

(1)



∂ (k) ∂ (kui ) ∂ = + ∂t ∂xi ∂xj



∂ (ε) ∂ (εui ) ∂ = + ∂t ∂xi ∂xj



∂ui uj ∂ui ∂ui uj 1 ∂p =− + ∇ 2 ui + gi +  ∂t ∂xj ∂xi ∂xj

Realizable k-ε turbulence model [28,29] including two equations for turbulent kinetic energy and turbulent kinetic energy dissipation rate is written as

(2)

+C1 Sε − C2 in which,

where ui represents the instantaneous velocity component in i direction, for example u and w are the velocities in x and z direction, respectively, while ui is the fluctuation velocity component in i direction, xi is the space coordinate in i direction, gi is the gravitational acceleration in i direction (0 and g for x and z directions, respectively), t is time, p is pressure and  and  are the density and kinematic viscosity of fluid, respectively.



C1 = max 0.43, =S



k ε

S = 2Sij · Sij

t + k



 ∂k 

t + ε

∂xj

+ Gk + Gb − ε

 ∂ε  ∂xj

ε2 ε + C1ε (1 − C3ε ) Gb √ k k + ε +5

(3)

(4)

 (5) (6)

1/2

(7)

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

221

Fig. 5. Computational mesh: (a) overall grid distribution of the computational domain; (b) close-up of the grids near the water surface; (c) close-up of the grids near the leak.

Fig. 4. The process of wave generation (contours are presented in water volume fraction) and the variation of water surface amplitude at four representative times with time interval of 24 s (comparison between the theoretical and numerical results): (a) t = 24 s; (b) t = 48 s; (c) t = 72 s; (d) t = 96 s.

Sij =

1 2



∂ui ∂uj + ∂xj ∂xi

∂uj Gk = −ui uj ∂xi 

Gb = −gi



t ∂ Pr t ∂xi

t = C

k2 ε

(8)

(9)

(10)

(11)

where k and ε represent turbulent kinetic energy and turbulent kinetic energy dissipation rate per unit mass, respectively, is the relative strain parameter, S is the strain rate, Gk and Gb represent production term of turbulent kinetic energy due to the average velocity gradient and production term of turbulent kinetic energy due to the lift, respectively,  =  is the dynamic viscosity of fluid, t is the turbulent viscosity, Prt is Prandtl number taken as 0.85, C , C1ε , C2 and C3ε are empirical model constants taken as 0.09, 1.44, 1.9 and 0.9, respectively, and  k and  ε are turbulent Prandtl numbers taken as 1.0 and 1.2, respectively. VOF model is an Eulerian multiphase flow model, which consists of one momentum equation for the mixture of the phases and one equation for the volume fraction of fluid. In the present work, ␣

is used to represent the fractional volume of a cell occupied by the liquid phase [30,31], i.e., ˛w and ˛o are the region occupied by water and oil, respectively. A unit value of ˛w corresponds to a cell full of water, while a zero value indicates that the cell contains no water. The fraction functions ˛w and ˛o are defined as ˛w =

Vw Vc

(12)

˛o =

Vo Vc

(13)

where Vc , Vo and Vw represent volume of a cell, volume of oil inside the cell and volume of water inside the cell, respectively. The two-dimensional transport equations for the fraction functions are given by

∂˛w ∂u˛w ∂w˛w + + =0 ∂t ∂x ∂z

(14)

∂˛o ∂u˛o ∂w˛o + + =0 ∂t ∂x ∂z

(15)

Then, the density and viscosity of the fluid in a cell can be expressed as  = (1 − ˛w − ˛o ) a + ˛w w + ˛o o

(16)

 = (1 − ˛w − ˛o ) a + ˛w w + ˛o o

(17)

where subscript a, o and w represent air, oil and water, respectively. Based on an analysis of the statistical wave data of the Bohai Sea, the second-order Stokes wave may be a good representative reflecting the statistical average feature. Therefore, the secondorder Stokes wave is adopted in simulation to examine the effect of

222

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

wave on surface drift. The potential function is expressed as [32,33]

cosh [2k(z + H)]

2 cosh 3Hw [2k(z + H)] + sin [2(kx − t)] 16Tw sinh4 (kH)

−Hw k + cos(kx − t) 2 sinh 2(kh)

Hw k cosh(kh)(2 cosh2 (kh) + 1) + cos 2(kx − t) 4 sinh3 (kh)

sinh4 (kH)

(18)

where Hw , Lw and Tw are the wave height, wavelength and wave period, respectively, which are taken as 1 m, 54 m and 6 s, respectively [24,25], H is the still water depth, k is the wave number ),  is the circular frequency of the wave ( = 2 (k = 2 ). Wave Lw Tw surface function = (x, t) is defined as



2 32 Hw ∂ Hw cosh [k(z + H)] cos(−t) + = Tw 4Tw Lw sinh(kH) ∂x

uw =

Hw Lw cosh [k(z + H)]

(x, z, t) = sin(kx − t) 2Tw sinh(kH)

(x, t) = Hw

in which

ww =

(19)

2.2. Numerical approach The model is implemented in the well-known CFD code FLUENT (V.15.0) that utilized finite volume method (FVM) to discretize the relevant partial differential equations. The coupling between the pressure and velocity fields is solved by the SIMPLE algorithm [34] to satisfy the conservation law of mass. The first-order (1st) accurate backward scheme in the time domain and the secondorder (2nd) accurate upwind scheme in the spatial domain are applied in simulations due to their stability and veracity [35]. The non-dimensional time-step is defined as uin t/D = 0.001 to meet the requirement of Courant-Friedrichs-Lewy (CFL) number being smaller than one. The convergent criteria for all calculations are defined as the residual in the control volume for each equation is smaller than 10−6 . 2.3. Computational domain As shown in Fig. 1, the computational domain is a rectangle region, and the overall is 300 m in the streamwise direction (x direction) and 25 m along the vertical direction (z direction). The coordinate origin is located at the lower left corner. Water occupies the lower region with height of 18 m, while air occupies the upper region. The wave absorption region is a relatively small rectangle region adjacent to the right exit, and is symmetric about the still water surface. The length is equal to the wave length (LD = Lw = 54 m), and the height is twice the wave height (hD = 2Hw = 2 m). The length of the computational domain is large enough to capture more than four complete waves. The leaking point of the damaged submarine pipe is located on the sea bed 100 m downstream of the left inlet. The size of the leakage hole (d) is a variable ranging from 0.02 m to 0.10 m with increment of 0.02 m, in order to examine the effect of leak size.

(22)

2 ∂ Hw sinh [k(z + H)] 32 Hw sin(−t) + = Tw 4Tw Lw sinh(kH) ∂z

sinh [2k(z + H)] sinh4 (kH) |x=0 = Hw



cos 2(−t)



sin 2(−t)

(23)

−Hw k + cos(−t)+ 2 sinh 2(kh)

Hw k cosh(kh)(2 cosh2 (kh) + 1) cos 2(−t) 4 sinh3 (kh)

 (24)

where uc , uw , ww are the velocity of current, the horizontal velocity of wave and the vertical velocity of wave, respectively, and H is the water depth. In this study, uniform current flowing along the horizontal direction with uc = 0.1 m/s is considered. The two velocity components of wave in the inlet boundary are obtained by the derivative of the potential function (Eq. (18)) and substituting x = 0 into the formulas. The direction of the wave is the same as the current. Under the water surface, the inlet velocity includes the velocity of wave and the velocity of current, while above the water surface, the velocity is equal to zero due to the still air. Waves are generated at the boundary with the defined inlet velocity (Eqs. (22) and (23)). Eqs. (20)–(24) are translated by a program into a series of statements written in C++ language, which is called User-Defined Function (UDF) embedded into the code. At the right outlet, a linear static pressure profile is employed to meet the hydrostatic pressure distribution as pout = g(H-z) for z ≤ H + , while the pressure is equal to 0 hPa (relative pressure, absolute pressure is 1013.25 hPa) for z ≤ H + . This outlet boundary is also defined by UDF. For the top boundary of the computational domain, the pressure is defined as 0 hPa so that the region above the water has the same pressure as the atmospheric pressure. Velocity inlet boundary condition is also employed in the oil leak with the flowing direction normal to the boundary (uo ). In order to reflect the release coming from a point, hydraulic diameter (equals to the diameter of the leak point) is defined in the velocity inlet boundary of the oil leak. The leaking rate is affected by the leak size, the pressure of water above the leak, the pressure within the damaged pipe and the pressure drop for oil flowing through the leak. Therefore, leaks in different pipes or in different locations at one pipe have different leaking rates. In this work, five rates ranging from 1 m/s to 5 m/s with increment of 1 m/s are selected to examine the effect of leaking rate. In the wave absorption region, the sponge layer damping absorber technique is applied by adding a flow resistance source item in the momentum equations [36,37]. 2

2

2.4. Boundary and initial conditions

∂u ∂u ∂u ∂ u ∂ u 1 ∂p + ) − a u +u +w =− + (  ∂x ∂t ∂x ∂z ∂x2 ∂z 2

The boundary conditions are also plotted in Fig. 1. The velocity inlet boundary condition is defined for the left inlet of the computational domain, which is expressed as

∂w ∂w ∂w ∂ w ∂ w 1 ∂p + ) − a w +u +w =− + g + (  ∂z ∂x ∂z ∂x2 ∂z 2 ∂t

 

u = uw + uc , u = 0,

z ≤H+ z>H+

w = ww ,

z ≤H+

w = 0,

z>H+

2

(25)

2

(26)

in which ˇ(x − x1 ) LD

(20)

a =

(21)

where a is the damping coefficient of wave absorption, x1 is the starting point coordinate of wave absorption region and ˇ is an empirical coefficient taken as 8 [38].

(27)

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

223

Fig. 6. Five selected mesh systems adopted to conduct the mesh independence check (the grid size decreases and the amount of cells increases gradually from mesh 1 to mesh 5): (a) mesh 1 (N = 60868); (b) mesh 2 (N = 89600); (c) mesh 3 (N = 193140); (d) mesh 4 (N = 382087); (e) mesh 5 (N = 573897).

Besides the boundary conditions, initial conditions are also defined before simulation. At initial time, the lower region with height of 18 m is filled with water and the upper region is filled with air. Both the water and air region are static so that the free surface (the interface of the two regions) is level. Since the water depth of Bohai is small, no ocean stratification is considered in the simulations. The densities of air and water are kept constant at 1.225 kg/m3 and 1025 kg/m3 , respectively, while the density of oil is a variable ranging from 750 kg/m3 to 990 kg/m3 in order to analyze the effect of oil density. The viscosities of air and water are also kept constant at 1.8 × 10−5 Pa·s and 1.003 × 10−3 Pa·s, respectively. Seven oil viscosities ranging from 1.0 Pa·s to 6.0 Pa·s are adopted to examine the effect of oil viscosity. In this study, the density and viscosity of the spilled oil and the leak size and leaking rate are variables. Twenty-four cases are simulated with the information listed in Table 1. A wide range of leaking flux has been considered with the minimum volume flux occurring in case 9 (uo = 1 m/s and d = 0.06 m) and the maximum volume flux appearing in case 8 (uo = 3 m/s and d = 0.10 m). 2.5. Model validation Oil spill experiments was carried out in a water flume to validate the numerical model. The water flume has a test section of 150 cm (length) × 10 cm (width) × 30 cm (height). In the tests, the water depth was set as 18 cm. An Acoustic Doppler Velocimetry (ADV) was employed to measure the flow profile, which was arranged 30 cm upstream, as shown in Fig. 2. An organic glass tube with diameter of 15 mm was arranged at the center of the channel bottom along the flow direction. A round hole with diameter of 4 mm was opened upwardly at the middle of the tube to model the leak. Two pressure sensors installed on both sides of the leak were used to monitor upstream and downstream pressures. The model oil was dyed by oil soluble dye. The density, viscosity and interface tension of model oil

are 905.7 kg/m3 , 0.0616 Pa s and 21.5 N/m, respectively. As shown in Fig. 2, a high speed camera on one side of the flume was used to record the oil spill process with sampling frequency of 100 frames per second. The size of each frame is 2048 × 1088 pixels, and the resolution is 4 pixels per mm. By changing the oil pump speed, oil spills under three different leaking pressure were tested. Fig. 3 shows the comparison of underwater spread of oil at different leaking pressures between experimental and numerical results. It is seen that both the dispersion pattern and migration distance of spilled oil coincide well with the experimental results. The rising height of oil droplets at different leaking pressures also agree with the test data, as shown in Fig. 3. The average error is less than 4%. Therefore, the numerical model can provide reasonable results for oil spill. At initial time, the water in the computational domain is static. The second-order Stokes wave is generated by the incident velocity in the left inlet, and then it propagates downstream until the wave absorption region, where the wave is disappeared due to the large dampness. Thus, it requires enough time to form wave in the domain. Fig. 4 shows the process of wave generation and propagation, which is presented in the contour of the volume fraction of water. The process of wave generation is displayed with time interval of 24 s, which equals to four wave periods. It is clearly seen that the inflow promotes the movement of water forming waves, and spread to the downstream gradually. The results show that the adopted damping absorber technique can avoid the reflection generated at the outlet boundary. The variation of water surface amplitudes at four representative times are compared between the numerical and theoretical results. As shown in Fig. 4, the generated wave almost reaches to the wave absorbing zone at t = 72 s. After the arrival, a stable periodic wave is formed. At t = 96 s, the wave is still stable, which further indicates the outstanding performance of the wave damper in the absorption

224

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 7. The rising height of oil droplets during the floating process at different mesh systems (mesh 1 (N = 60868), mesh 2 (N = 89600), mesh 3 (N = 193140), mesh 4 (N = 382087) and mesh 5 (N = 573897), contours are presented in water volume fraction).

of the coming wave. The coincidence of water surface amplitudes indicates the wave generation technique is effective.

2.6. Computational mesh As shown in Fig. 5, the computational domain is divided into five blocks. The water region above the leak from x = 97 m to x = 103 m is discretized with triangular cells, and progressive mesh is used to capture the near-leak flow properties as shown in Fig. 5(c). The closer the leak, the smaller is the grid size. The other four blocks of the computational domain are discretized with quadrilateral cells. As shown in Fig. 5(b), in order to capture the properties of surface wave, the grids of the narrow area below and above the still water level with height of 2 m (z = 17–19 m) are refined with half size. The required time of oil floating and the horizontal migration distance of spilled oil are two key parameters concerned in this work. Therefore, the microscopic details of dispersed droplets are not important. The mesh independence check has been conducted by identifying the changes in the two key parameters. Five selected mesh systems are compared as shown in Fig. 6. From mesh 1 to mesh 5, the grid size decreases and the amount of cells increases gradually. The grid size of mesh 5 is the same as the scale used in the

model validation. Therefore, the dispersed oil droplets are captured clearly in mesh 5 (Fig. 7). However, the overall trend of underwater spread and the rising height are mostly the same. Since the CPU time increases significantly with the increasing of computational cells, mesh 1 with a small number of cells (N = 60868) and the same calculation results for the concerned parameters is chosen in the following simulations.

3. Results and discussions 3.1. Effect of marine environment The oil spills are under the action of complex ocean environment including wind, wave and current, which makes the underwater spread and surface drift of oil spills intricate. The combined action of wave and current is considered in this paper, and the environments including still water, only current, only wave are also considered to make a comparison. Fig. 8 shows the floating process of oil spilled from the leak to free surface at different marine environments. In still water, there is a continuous vertical oil stream just releasing from the leak with the height of about 3 m. When the rising height further increases,

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

225

Fig. 8. The rising height and underwater spread of oil droplets during the floating process at different marine environments (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave, contours are presented in water volume fraction).

the oil stream becomes diffuse in the form of oil droplets or droplet groups. However, the diffusion is symmetrical about the vertical line across the leak center. It requires 40 s for oil droplets to reach the free surface. When it comes to the only current environment, the continuous oil stream releasing from the leak is also presented, but it is inclined

to the downstream. Under the action of gravity, inertia force, buoyancy and shear stress, the oil stream is tore apart into droplets at a certain height (about 2.5 m) by the current. The oil droplets become more dispersed with the increase in rising height, and the droplets are migrating to the downstream while floating, which is also found in our previous study [17]. Due to the longer length of rising path, it

226

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 9. The horizontal migration distance and surface drift of oil droplets during the drifting process at different marine environments (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave): (a) still water; (b) only current; (c) only wave; (d) wave and current (contours are presented in oil volume fraction).

spends more than forty seconds to reach the free surface. As shown in Fig. 8, under the action of wave, the underwater spread presents more dispersed due to the oscillation of water particles. When there is only wave, the rising path of oil droplets sometimes is tilted to the upstream and sometimes is tilted to the downstream in the floating process. However, taken as a whole, it is still symmetrical about the vertical line across the leak center. Fifty seconds is required for the spilled oil to present on the water surface. Under the combined action of wave and current, the dispersed oil droplets are migrating to the downstream while floating and oscillating. Although the time required for the whole floating process is increased to 60 s, the horizontal migration distance is larger than other cases. Therefore, there is a larger scope of underwater diffusion as the superposition of wave and current, which becomes more difficult to control. Fig. 8 also plots the comparison of the rising height of oil droplets during the floating process under different marine environments. We can see that the spilled oil under the uniform current has the same rising height as that in still water in the initial 5 s. The reason is

the oil just released from the leak has relatively larger inertia, which is beyond the shear action of current. However, with the increasing of rising height, the inertia of spilled oil is gradually weakened by the viscous shearing stress. As the rising height increased to 4 m, the shear action of the current begins to appear. Then the oil droplets are migrating to the downstream by the current while floating as mentioned in [14], and the slope of the rising height curve becomes smaller, resulting in more time required to reach the surface. For second-order Stokes wave in shallow water, there is water particle motion in the bottom of the sea. Under the wave shearing action, the rising rate of oil droplets becomes slower. Since the horizontal maximum velocity of wave particle (0.7m/s) is far greater than the current rate, the shear stress imposed by the wave is greater than that imposed by the current. Therefore, the slope of the rising height curve for spilled oil in wave environment is smaller than that in current, as shown in Fig. 8. The rising rate of oil droplets is further slowed under the combined shearing action of wave and current. So the corresponding curve has the smallest gradient about 0.3 m/s.

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

227

Fig. 10. The rising height and underwater spread of oil droplets during the floating process at different leak sizes under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave, contours are presented in water volume fraction).

After the oil droplets reaching the free surface, they will continue to drift on the water surface. Fig. 9 shows the drifting process of spilled oil at different marine environments. It is clearly seen that the size of oil films drifting on water surface without wave is larger than that subjected to the wave. The fluctuation of water surface

and the oscillation of water particles contribute to the dispersion of oil films. As shown in Fig. 9, under the combined action of wave and current, the spilled oil has the fastest surface drifting velocity, which takes 161 s to arrive at the wave absorbing zone.

228

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 11. The horizontal migration distance and surface drift of oil droplets during the drifting process at different leak sizes under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave): (a) d = 0.02 m; (b) d = 0.04 m; (c) d = 0.06 m; (d) d = 0.08 m; (e) d = 0.10 m (contours are presented in oil volume fraction).

The horizontal migration distances of spilled oil during the overall floating and drifting process at different marine environments are depicted in Fig. 9. The values in the brackets of the legend represent the coordinates of the turning points between the two processes, which provide the information about when and where the spilled oil reaches the water surface. Making overall leaking time (floating time and drifting time) as abscissa and horizontal displacement as vertical ordinate, the curve slope represents the horizontal migrating velocity. In still water, the spilled oil stream spreads like a fountain. The horizontal migration distance linearly increases with time during the drifting process, and the surface diffusion rate is about 0.59 m/s. When there is a uniform current, the surface drifting rate becomes to be 0.70 m/s, which can be approximately considered as the superposition of the diffusion velocity in still water (0.59 m/s)

and the velocity of the current (0.10 m/s). Under the action of wave, the surface drifting rate increases to 1.16 m/s. The velocity of water particles in the wave surface contributes to this increase in drifting rate. The oil droplets have the largest horizontal migration distance of 18 m during the floating process under the condition of current combining with wave, and the spilled oil migrates with the speed of 1.27 m/s during the drifting process. Therefore, the spilled oil spreads very quickly under the superposition of wave and current and fast responses are required in this condition. 3.2. Effect of leaking flux The volume flux of spilled oil is mainly determined by the size of the leak and the leaking rate. Fig. 10 displays the floating process of oil spilled from different sizes of leaks at the same leaking rate

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

229

Fig. 12. The rising height and underwater spread of oil droplets during the floating process at different leaking rates under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave, contours are presented in water volume fraction).

uo = 3 m/s by the wave-current combined effect. The larger the leak size, the faster is the rising rate of spilled oil. This phenomenon is attributed to the great leaking momentum of oil spilled from a large

leak. For oil spilled from the leak d = 0.10 m, 49 s is spent to complete the floating process, which is about 60% of the time required for oil spilled from a small leak d = 0.02 m. Meanwhile, due to the large

230

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 13. The horizontal migration distance and surface drift of oil droplets during the drifting process at different leaking rates under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave): (a) uo = 1m/s; (b) uo = 2m/s; (c) uo = 3m/s; (d) uo = 4m/s; (e) uo = 5 m/s (contours are presented in oil volume fraction).

leakage, the oil spilled from the large leak spreads faster and wider during the floating process, bringing about a more serious threat. The rising height of oil droplets during the floating process at different leak sizes are also given in Fig. 10. Within 6 m height above the seabed, the rising height of spilled oil increases linearly with time, and the curve in Case 8 (d = 0.1 m) has the largest slope. The inertia mainly controls the floating of spilled oil near the seabed. As the oil droplets further rise, the underwater spread begins to be affected by the wave-current, and the rising-height curves exhibit a certain fluctuations. When the spilled oil is close to the surface, there is an obvious decrease in the rising speed. The severe oscillation of wave particles in the surface contributes to the reduction in

the slope of the rising height curves. In general, the oil spilled from the largest leak has the highest rising height at the same time. Fig. 11 shows the surface drifting process of oil spilled from different sizes of leaks. The size and shape of dispersed oil droplets drifting in the surface are similar in the five cases. However, the oil spilled from large leak has relatively dense distribution on the water surface. Because the times consumed in floating process are different, the total times required for spilled oil to reach the wave absorption zone are also different. However, there is no obvious difference for the consumed time in the drifting process, which ranges from 101 s to 115 s. The horizontal migration distances during the overall process at different leak sizes are shown in Fig. 11. Although the horizon-

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

231

Fig. 14. The rising height and underwater spread of oil droplets during the floating process at different oil densities under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave, contours are presented in water volume fraction).

232

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 15. The horizontal migration distance and surface drift of oil droplets during the drifting process at different oil densities under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave): (a) o = 750 kg/m3 ; (b) o = 850 kg/m3 ; (c) o = 910 kg/m3 ; (d) o = 930 kg/m3 ; (e) o = 950 kg/m3 ; (f) o = 970 kg/m3 ; (g) o = 990 kg/m3 (contours are presented in oil volume fraction).

tal migration velocities are different in the underwater floating process, they become similar in the surface drifting process. Besides the leak size, the leaking rate under wave-current combined effect is also considered in this paper. The floating process of oil spilled at different leaking rates is plotted in Fig. 12. We can see that the faster the leaking rate, the shorter is the required time

to reach the water level. Eighty-two seconds is spent for the floating process of oil spilled with uo = 1 m/s, while it is reduced to 52 s for oil spilled with uo = 5 m/s. Moreover, large leaking rate means a large amount of leakage, leading to high dispersion concentration under the water surface. The rising height of oil droplet increases linearly near the seabed. Then a certain fluctuation appears until

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

233

Fig. 16. The rising height and underwater spread of oil droplets during the floating process at different oil viscosities under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave, contours are presented in water volume fraction).

H = 16 m. The action of wave is highlighted near the water surface, which reduces the rising rate of spilled oil. As shown in Fig. 13, the drifting process and surface distribution are similar to Fig. 11. Large volume flux (uo = 5 m/s) has relatively dense distribution of oil droplets in the surface. However, the time spent in the drifting process is similar to the other four cases with relatively small large volume flux. The same horizontal migration velocities during the drifting process are confirmed in Fig. 13, which exhibits the same growth rate of horizontal migration distances. 3.3. Effects of oil density and viscosity The physical characteristics of the spilled oil including density and viscosity are also taken into consideration in this work. Fig. 14 shows the rising height of oil droplets versus the time at different oil densities. Under the same wave-current condition, the greater the oil density, the longer is the required time for rising. The relatively light oil spreads faster and wider at the same time. For the relatively heavy oil with density of 990 kg/m3 , a part of spilled oil spreads close to the seabed under the action of gravity. Therefore, for heavy oil leakage, not only the surface recovery but also the seabed pollution should be considered in emergency response. The rising height increases linearly within 6 m height above the seabed.

As the height further increases, the curve slope presents a reduction accompanied by a certain fluctuation, and it has an obvious reduction near the water surface. The greater the oil density, the smaller is the curve slope. The horizontal migration distances during the overall process at different oil densities are shown in Fig. 15. We can see that low density oil has relatively dense distribution on the water surface during the drifting process, while high density oil has wide distribution under the water. The consumed time for spilled oil to reach the wave absorption zone increases with the increasing of oil density. For oil with density of 750 kg/m3 , 82 s is spent in the drifting process, while it takes 112 s for oil with density of 990 kg/m3 . The horizontal displacement increases near linearly, and the growth is relatively slow for high density oil. The difference is significant with gradients of 1.63 m/s and 1.14 m/s for oil densities of 750 kg/m3 and 990 kg/m3 , respectively. As show in Fig. 16, seven curves of different viscosities are almost superposition. It is seen that there are no obvious differences in rising rate and underwater diffusion pattern. For the seven cases, one minute is consumed in the floating process. Under the shearing action of wave and current, spilled oil presents in dispersed oil droplets, which follow with the water flow. In this condition, the

234

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

Fig. 17. The horizontal migration distance and surface drift of oil droplets during the drifting process at different oil viscosities under the combined action of wave and current (uc = 0.1 m/s, Hw = 1 m, Lw = 54 m and Tw = 6 s, the current is in the same direction as the wave): (a) o = 1.0 Pa·s; (b) o = 1.5 Pa·s; (c) o = 2.0 Pa·s; (d) o = 2.5 Pa·s; (e) o = 3.0 Pa·s; (f) o = 4.5 Pa·s; (g) o = 6.0 Pa·s (contours are presented in oil volume fraction).

viscous force is not the main force, which has very little effect on the spreading process. Fig. 17 shows the rising height of oil droplets versus the time at different oil viscosities. Similarity to the floating process, the surface drifting barely has difference. The consumed time for the drifting process ranges from 101 s to 106 s. The wave breaking and

turbulent mixing are not considered in the model. This limitation may be one reason for the little effect of viscosity. 4. Conclusions The present study aims to investigate the effects of physical ocean environment, leakage parameters and oil physical character-

H. Zhu et al. / Applied Ocean Research 64 (2017) 217–235

istics on the underwater spread and surface drift of oil spilled from submarine pipeline. The 2D RANS approach and realizable k-ε turbulence model coupling with VOF model are utilized to solve the multiphase flow, which are validated by experiments. Although the results are obtained in specific conditions, they can still provide a certain references, and the numerical method can be employed for other conditions. Based on the results of the present investigation, the following conclusions can be made: (1) In still water, the spilled oil stream spreads like a fountain, and the diffusion is symmetrical about the vertical line across the leak center. It spends the least time for oil to reach the surface. Under the action of current, the spilled oil stream is tore apart into droplets at the height of 2.5 m, and the oil droplets migrate to the downstream while floating, resulting in the increase in rising path. Under the action of wave, the underwater spread of oil appears more dispersed due to the oscillation of wave particles. However, the wave particle velocity in the surface contributes to the increase in drifting rate. As the superposition of wave and current, there is a larger scope of underwater diffusion, and the oil droplets have the slowest rising rate but the fastest surface drifting rate. (2) The inertia of spilled oil controls the floating near the seabed. With the rising of oil droplets, the underwater spread becomes to be affected by the wave-current. The action of wave is highlighted near the water surface, and the oscillation of wave particles contributes to the reduction in the rising speed of oil near the surface. (3) The leakage parameters including the leak size and the leaking rate affect the floating time and dispersion concentration, while the ocean environment affects the horizontal migration and surface drifting. In the floating process, the oil spilled with large leakage spreads faster and wider, leading to a more serious threat. However, there are no obvious differences in the consumed time and horizontal displacement in the drifting process for different leaking fluxes. (4) Oil density has obvious effect on the underwater spread but limited effect on the surface drifting. Low density oil rises quickly, leaving short time for response. On the other hand, not only the surface recovery but also the seabed pollution should be considered in emergency response for heavy oil leakage. Both the underwater spread and surface drifting are almost unaffected by the oil viscosity. One reason is the viscous force is not the main force for the spilled oil presented in dispersed droplets. Another reason may be the limitation of the model. Acknowledgements The research work was supported, in part, by National Nature Science Foundation of China (No. 11502220) and Scientific Research Starting Project of SWPU (No. 2014PYZ001). References [1] M. Reed, P. Daling, A. Lewis, M. Kristin-Ditlevsen, B. Brørs, J. Clark, D. Aurand, Modelling of dispersant application to oil spills in shallow coastal waters, Environ. Modell. Softw. 19 (2004) 681–690. [2] M. Gade, W. Alpers, Using ERS-2 SAR for routine observation of marine pollution in European coastal waters, Sci. Total Environ. 237 (4) (1999) 441–448. [3] I. Keramitsoglou, C. Cartalis, C.T. Kiranoudis, Automatic identification of oil spills on satellite images, Environ. Modell. Softw. 21 (5) (2006) 640–652. [4] G. Ferraro, S. Meyer-Roux, O. Muellenhoff, M. Pavliha, J. Svetak, D. Tarchi, K. Topouzelis, Long term monitoring of oil spills in European seas, Int. J. Remote Sens. 30 (3) (2009) 627–645.

235

[5] O. Garcia-Pineda, B. Zimmer, M. Howard, W. Pichel, X.F. Li, I.R. MacDonald, Using SAR image to delineate ocean oil slicks with a texture classifying neural network algorithm (TCNNA), Can. J. Remote Sens. 35 (5) (2009) 411–421. [6] P. Liu, X. Li, J. Qu, W. Wang, C. Zhao, W. Pichel, Marine oil spill detection with fully polarimetric UAVSAR data, Mar. Pollut. Bull. 62 (12) (2011) 2611–2618. [7] M. Migliaccio, F. Nunziata, C. Brown, B. Holt, X.F. Li, W. Pichel, M. Shimada, Polarimetric synthetic aperture radar utilized to track oil spills, Eos Trans. Am. Geophys. Union 93 (16) (2012) 161–162. [8] M. Konik, K. Bradtke, Object-oriented approach to oil spill detection using ENVISAT ASAR images, ISPRS J. Photogramm. 118 (2016) 37–52. [9] J.A. Galt, Trajectory analysis for oil spills, J. Adv. Mar. Technol. Conf. 11 (1994) 91–126. [10] M. Reed, Ø Johansen, P.J. Brandvik, P. Daling, A. Lewis, R. Fiocco, D. Mackay, R. Prentki, Oil spill modeling towards the close of the 20th century: overview of the state of the art, Spill Sci. Technol. Bull. 5 (1) (1999) 3–16. [11] B. Hackett, E. Comerma, P. Daniel, H. Ichikawa, Marine oil pollution prediction, Oceanogr 22 (3) (2009) 168–175. [12] M. Marta-Almeida, M. Ruiz-Villarreal, J. Pereira, P. Otero, M. Cirano, X.Q. Zhang, R.D. Hetland, Efficient tools for marine operational forecast and oil spill tracking, Mar. Pollut. Bull. 71 (1–2) (2013) 139–151. [13] X.B. Chao, N.J. Shankar, H.F. Cheong, Two- and three dimensional oil spill model for coastal waters, Ocean Eng. 28 (12) (2001) 1557–1573. [14] S.D. Wang, Y.M. Shen, Y.H. Zheng, Two-dimensional numerical simulation for transport and fate of oil spills in seas, Ocean. Eng. 32 (13) (2005) 1556–1571. [15] P. Tkalich, A CFD solution of oil spill problems, Environ. Modell. Softw. 21 (2) (2006) 271–282. [16] W. Li, X. Liang, J.G. Lin, Mathematical model and computer simulation for oil spill in ice waters around island based on FLUENT, J. Comput. 8 (4) (2013) 1027–1034. [17] H.J. Zhu, P.Z. Lin, Q. Pan, A CFD (computational fluid dynamic) simulation for oil leakage from damaged submarine pipeline, Energy 64 (1) (2014) 887–899. [18] J.M. Sayol, A. Orfila, G. Simarro, D. Conti, L. Renault, A. Molcard, A Lagrangian model for tracking surface spills and SaR operations in the ocean, Environ. Modell. Softw. 52 (2014) 74–82. [19] H.X. Wang, J.L. Xu, W.K. Zhao, J.Q. Zhang, Effects and risk evaluation of oil spillage in the sea areas of Changxing Island, Int. J. Environ. Res. Public Health 11 (8) (2014) 8491–8507. [20] J. Papadimitrakis, M. Psaltaki, M. Christolis, N.C. Markatos, Simulating the fate of an oil spill near coastal zones: the case of a spill (from a power plant) at the Greek Island of Lesvos, Environ. Modell. Softw. 21 (2006) 170–177. [21] P. Neofytou, A.G. Venetsanos, D. Vlachogiannis, J.G. Bartzis, A. Scaperdas, CFD simulations of the wind environment around an airport terminal building, Environ. Modell. Softw. 21 (2006) 520–524. [22] B.A. Wols, J.A.M.H. Hofman, W.S.J. Uijttewaal, L.C. Rietveld, J.C. Van-Dijk, Evaluation of different disinfection calculation methods using CFD, Environ. Modell. Softw. 25 (2010) 573–582. [23] G. Maggiotto, R. Buccolieri, M. Antonio-Santo, L. Sandra-Leo, S. Di-Sabatino, Validation of temperature-perturbation and CFD- based modelling for the prediction of the thermal urban environment: the Lecce (IT) case study, Environ. Modell. Softw. 60 (2014) 69–83. [24] X.C. Lv, D.K. Yuan, X.D. Ma, J.H. Tao, Wave characteristics analysis in Bohai Sea based on ECMWF wind field, Ocean Eng. 91 (2014) 159–171. [25] W.J. Guo, G.X. Wu, B.C. Liang, T.J. Xu, X.B. Chen, Z.W. Yang, M.X. Xie, M.R. Jiang, The influence of surface wave on water exchange in the Bohai Sea, Cont. Shelf Res. 118 (2016) 128–142. [26] J.O. Hinze, Turbulence, McGraw-Hill, USA, New York, 1975. [27] U. Frisch, Turbulence, Cambridge University Press UK, London, 1995. [28] V.M. Karim, M. Bart, Application of two buoyancy-modified k-␧ turbulence models to different types of buoyant plumes, Fire Saf. J. 41 (2) (2006) 122–138. [29] P. Rohdin, B. Moshfegh, Numerical predictions of indoor climate in large industrial premises: a comparison between different k-␧ models supported by field measurements, Build. Environ. 42 (11) (2007) 3872–3882. [30] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (81) (1981) 201–225. [31] P.D. Hieu, T. Katsutoshi, V.T. Ca, Numerical simulation of breaking waves using a two-phase flow model, Appl. Math Model. 28 (11) (2004) 983–1005. [32] P. Lin, Numerical modeling of water waves, Comput. Phys. 3 (4) (2008) 141–147. [33] P. Lin, L.F. Liu, Internal wave-maker for Navier-Stokes equations models, J. Waterw. Port Coast. 125 (4) (1999) 207–214. [34] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, USA, New York, 1980. [35] W.B. Jeffrey, D.H. William, Upwind schemes for the wave equation in second-order form, J. Comput. Phys. 231 (231) (2012) 5854–5889. [36] Z. Hu, W.Y. Tang, H.X. Xue, X.Y. Zhang, J.T. Guo, Numerical simulations using conserved wave absorption applied to Navier-Stokes equation model, Coast. Eng. 99 (2015) 15–25. [37] T.W. Hsu, S.H. Ou, B.D. Yang, I.F. Tseng, On the damping coefficients of sponge layer in Boussinesq equations, Wave Motion 41 (1) (2005) 45–57. [38] Z. Hafsia, M.B. Hadj, H. Lamloumi, K. Maalel, Internal inlet for wave generation and absorption treatment, Coast. Eng. 56 (9) (2009) 951–959.