# Key Technology of Oil Spill Early-warning and Forecasting

## Research content

China’s integrated marine forecast system is based on the internationally advanced marine numerical model FVCOM and atmospheric numerical model WRF. It carries out high-resolution accurate simulation of the current field and wind field of the researched sea area, completes seamless integration of current field and wind field with the OILMAP software, upgrades the OILMAP marine environment database, and realizes rapid forecast, early warning and trace-back of oil contaminants in the Bohai Sea and North Yellow Sea (Figure 4.1) [67, 68, 69]. Up to now, the system has been operating on the server of Shandong Maritime Safety Administration for one year. It is currently the only operational system of the Shandong Maritime Safety Administration which is used in forecast and emergency decision-making of oil spill accidents. The system has the function of unattended automatic forecasting, which is easy to use. The simulated marine dynamic environment field is timely and accurate, and the predicted oil spill drift trajectory is reliable.

FIGURE 4.1: Technical roadmap of China’s integrated marine forecast system

## Forecast module of the marine dynamics environment

### Introduction to numerical model

1) Three-dimensional Finite-Volume Coastal Ocean Numerical Model (FVCOM)

The model used in the system is Finite-Volume Coastal Ocean Model (FVCOM) with triangular grid, finite volume and three-dimensional primitive equations. The greatest feature and advantage of FVCOM. which is adopted for coastal and estuarine tidal circulation, is that it combines the finite-element method, which fits the boundary easily and enables the local refinement, and the finite-difference method, which facilitates the discrete calculation of the marine primitive equations [70, 71, 72].

2) Unstructured grid

Structured grid means that all internal points in the grid area have the same neighboring cells. In contrast to the definition of the structured grid, unstructured grid means that the internal points in the grid area do not have the same neighboring cells. In other words, the number of the grids connected to different internal points in the grid area is different. It can be seen from the definition that the structured grid and unstructured grid have overlapping parts, namely the unstructured grid may contain parts of the structured grid [73, 74, 75].

3) Design of triangular grid

Similar to the finite-element method, the calculation area is divided into some non-coincident triangular cells. Each triangular grid consists of three nodes, one center and three sides. In order to more accurately calculate the water level, velocity, temperature and salinity on the sea surface, the numerical

FIGURE 4.2: Distribution of variable space configuration calculation variables on the grid points calculation is carried out on a specially designed triangular grid. The variables on a node are calculated by the net flux of the line that connects such point, the center of the triangle and the center of the side, and the variables on the center are calculated by the net flux of the side of the triangle. The position of the variables is as shown in Figure 4.2. *и* and *v* represent the velocity, and F represents water level, temperature, salinity and density [76].

4) Examples of Bohai sea grid

The Bohai Sea is a near-closed shallow sea of continental shelf in the north of China. It communicates with the Yellow Sea through the Bohai Strait in the east. Its north, west and south sides are surrounded by land with complex terrain. With a large population and developed economy in coastal areas, the sustainable utilization of the Bohai Sea is one of the strategic concerns of the Chinese people, and the distribution and structural characteristics of marine current, temperature field and salinity field are also one of the important contents of physical marine studies in this sea area. It can be concluded by using surface-water modeling system (SMS) software to design the grid of the Bohai Sea (Figure 4.3) that the triangular grid can conveniently and correctly fit the complicated coastline and refine some parts of the grid. The minimum resolution at the open boundary is about 12,000 m and the complicated coastline, bay mouth and key port areas can be refined to make the maximum resolution 5 nr, thus meeting the requirements on resolution in general marine phenomenon processes.

FIGURE 4.3: Example of grid of the calculation sea area

5) Basic control equations

The FVCOM control equations consist of the momentum equation, continuity equation, temperature equation, salinity equation and density equation [77, 78].

The basic equations describing the hydrodynamic and thermal systems included the Navier-stokes equation about motion, continuity equation describing mass conservation, state equations of seawater and the conservation equation of temperature and salinity [63]. In order to apply these equations to the shallow sea of the continental shelf, some assumptions and approximations are made for these equations according to some characteristics of the shallow sea:

- a) In general fluid motion, the compressibility of fluid depends on the square of the velocity-to-sound velocity ratio. The speed of seawater movement is far slower than the sound velocity, so the seawater can be considered incompressible.
- b) The ratio of depth scale
*(D)*to horizontal scale (L) in coastal waters is a small scalar, namely*о*=*D/L*< l. Therefore, the vertical equation of motion can be approximated by static equilibrium. In other words, gravity equilibrates to vertical pressure gradient force.

In consideration of the right-hand Cartesian coordinate system, the eastward direction is the ж-axis positive direction, the northward is the «/-axis positive direction, and the vertical upward is the positive direction of the 2-axis. The free sea surface is *z = i}(x, y. t):* the bottom topography is *z = —H{x,y).*

*V* is the horizontal velocity vector, *V =* (*U, V),* and V is the horizontal gradient operator, then the continuity equation can be written as follows:

The momentum equation is:

The vertical static equilibrium equation is:

The conservation equation of temperature and salinity is:

The state equation of seawater is:

where *po* represents reference density of seawater, *p* represents the field density, *g* represents the gravity acceleration, *P* represents the pressure, Km represents the vertical eddy viscosity coefficient, / represents Coriolis parameter applying ,5 plane approximation, *в* represents potential temperature (in the shallow sea, it is spot temperature), s represents salinity, and *К и* represents vertical eddy diffusivity. Km and *Кц* are obtained through the Mellor-Yamada level-2.5 turbulence closure model. The pressure at the 2 depth *P[x. y, z, t*) can be obtained through integral equation.

It is assumed here that the atmospheric pressure *P _{a}tm* on the sea surface is constant.

Fx,Fy, Fe and *F _{s}* in the equations (4.2), (4.3), (4.5) and (4.6) can be written as follows:

It should he noted that *F* and *Fy* are invariant of coordinate rotation. Where Am is the horizontal eddy viscosity coefficient and is calculated with the equation proposed by Smagorinsky:

Based on needs, the parameter *a* can be 0.01 to 0.5 and is usually 0.1. By using the Prandtl number, the value of horizontal eddy diffusivity Ah can be obtained through Am-

The vertical eddy viscosity coefficient Km and the vertical eddy diffusivity *Кн* in Eqns. (4.2), (4.3), (4.6) and (4.7) can be determined through the Mellor-Yamada level-2.5 turbulence closure model (1974, 1982). In this way, the influence of artificial selection of eddy viscosity coefficient on the simulation of physical field can be overcome to some extent [79].

Such a model describes two physical quantities: turbulence kinetic energy *q*^{2}/2 and macroscopic scale of turbulence *l.* The equations are:

where *F _{q}* and

*F)*represent the horizontal diffusion terms of the turbulence kinetic energy and the macroscopic scale of turbulence, respectivety, and their form is similar to the temperature and salinity diffusion terms in Eqn. (4.11),

and *p _{s} = Km* [(f)

^{2}d

^{-}(

*yyr*) J is the turbulence kinetic-energy term caused by shear force;

*pi, = j^Kh^z*is the turbulence kinetic-energy term caused by buoyancy force;

*e*=

*q*is the turbulent kinetic-energy dissipation rate;

^{3}/BiL*ш*= 1 +

*E*

*-2*(jj^) is wall proximity function, and

*к =*0.4 is Von Karman constant.

At a location near the surface, *L* « *I/К.* hence *ui =* 1 +£2; at a location far away from the surface, *l* *ш* = 1. *Км, Кн* and *K _{q}* are determined by the following equations, respectively.

where, S’.у/, S// and *S _{q}* are stable functions and can be written as follows according to the research results of Mellor and Yamada as well as Galperin et ah:

where *S _{q} =* 0.20,

*G ц =*— is Richardson number, and

*N =*^ —

is Brunt-vaisala frequency.

Ai,A’*2**, В*i, *В**-2*, *Ci, Ei* and *E**2* are empirical constants and (Лу, Л2, *B, B**2**, Ci*, *Ei* and *E**2**) =* (0.92, 0.74, 16.6, 10.1, 0.08, 1.8, 1.33), respectively [79].

In the stable stratified fluid, the macroscopic scale of turbulence *l* should also satisfy the following [80]:

6) Boundary conditions

The above control equations can only be closed by giving appropriate initial conditions and boundary conditions. The boundary conditions include the conditions of sea surface, the seafloor boundary, the seafloor solid side boundary and the sea-side open boundary.

a) The sea surface boundary conditions

On the free sea surface, г = *r/(x, y, t):*

where *(т _{ох},т_{оу})* is the wind-stress vector on the sea surface, and its value is represented by

*u*The calculation equation is as follows:

_{Ts}.where *p _{a}* is air density and Co is wind drag coefficient. When the wind speed is less than 25 m/s, the wind drag coefficient is determined by the Large and Pond equation:

|W| is the sum of wind-speed vector. *H* is the net heat flux on the sea surface which includes four parts: solar radiation, long-wave radiation, sensible heat and latent heat. *S **= **S(0)(E **— **p)/po* is virtual salt flux, *(E* — *p)* is freshwater flux of sea surface evaporation or precipitation and 5(0) is the sea surface salinity.

b) Seafloor boundary conditions

On the seafloor, 2 = *—H(x,y):*

where *II(x, **y)* is the bottom topography, *(ть _{х},ть_{у})* is bottom friction stress vector on seafloor, and its value is represented by

*и*It is calculated with the following equation:

_{т}ь-*{Пх,П _{у}) = poC_{z}V_{b}Vb = poCzy/U*

^{2}+ P

^{2}(C,

*V),*(4.23)

Г -2 "I

where *C'z = MAX *j_{n}^/_{Zo})‘i_{;} 0-0025 is the bottom friction coefficient, *к* = 0.4 is Von Karman constant, zq is roughness parameters of seafloor.

c) The solid-side boundary conditions

It is generally assumed that the normal velocity vertical to the solid coast at the solid side of the researched area is zero. Namely: *U ■ n =* 0.

In addition, it is assumed that there is no heat and salt exchange between seawater and the solid boundary. Namely: щ = *щ =* 0

d) The open boundary conditions

On the open boundary of the sea area, the boundary conditions of the sea level are calculated from the harmonic constant of the main constituents at the boundary with the following equation:

where a* is the vibration amplitude of the constituent *i; oji* is the frequency of the constituent *i*; d>i is the angle of delay of the constituent г; and *E _{m}ean *is the residual water level at this point relative to the mean sea level.

In addition, the open boundary of velocity is given by a radiation condition or non-gradient boundary condition, and the open boundary condition of temperature and salinity is given by a non-gradient boundary condition or upwind convection scheme in the radiation inflow region of outflow region.

7) Vertical coordinate transformation

To better fit the bottom topography, *о* coordinate transformation is adopted in the model in the vertical direction. Namely:

where *H(x,y)* represents the bottom topography, *y(x,y*) represents the sea surface relief. At the point *z = y, a =* 0; at the point *z* = *—H,cr =* —1.

Assuming that *D = H +* 77, the transformation of the new and old coordinate is:

So the equation (4.19e) can be transformed into

In addition, any variable *G* which is integrated from the seafloor to the sea surface can be expressed as follows:

Therefore, Eqns. (4.1), (4.2), (4.3), (4.4), (4.5), (4.6) and (4.7) can be written as (for ease of representation, the asterisk is omitted):

where, uj is the vertical velocity vertical to the sigma layer after the coordinate transformation.

where *в,* represents *в, S,q ^{2}* and

*q*

^{2}l.8) Inner- and outer-mode splitting algorithm

In order to save computer time and improve calculation efficiency, FVCOM model adopts the method of inner- and outer-mode alternate calculation. The outer mode is two-dimensional, and the time step is shorter when calculating the water level and vertical average velocity. The inner mode is three- dimensional, and the time step is longer when calculating the physical quantities like current velocity, temperature, salinity and turbulence kinetic energy.

The equations for outer-mode calculation can be obtained by integrating Eqns. (4.34), (4.35) and (4.36) vertically from the seafloor to the sea surface:

where *— WU*(0), — *WV*(0) and *— WU*(—1), — *WV*(—1) are the sea surface wind stress and bottom stress component, respectively.

It is firstly assumed that all variables at the time points *t ^{n}~^{1}* and

*t"*have already been obtained, then the right-hand terms of the Eqns. (4.50) and (4.51) can be obtained, and it is considered that they remain unchanged from the time point

*t*to

^{n}*t*Equations (4.49), (4.50) and (4.51) and the time step of outer-mode DTE are used to show the several steps integral to the time point

^{n+1}.*t*thus the water level and vertical average velocity can be obtained for calculation of the inner mode.

^{n+l}.9) Lagrange residual current calculation module

Feng analyzed and pointed out that although the Euler residual current is a solenoidal vector field, it does not meet the condition of conservation of material surface composed of fluid micro-aggregates in the current field [81]. The mass transport speed is not only a solenoidal vector field, but can also be proved to satisfy the conservation equation of material surface in the current field in the most general sense. It is of practical significance in the lowest level to describe the circulation issues in the shallow sea by using the Lagrange residual current, which represents the velocity field constituted by mass transport velocity [82].

Euler is usually used to represent the shallow sea circulation, i.e., Euler residual circulation. It is defined as low-pass filtering at any selected point in the current field space, or simply time averaging over a tidal period, which is called Euler average velocity or Euler residual current.

Feng pointed out that mass transport velocity was derived as Lagrange residual current of the lowest level. The calculation process is: calculate the Euler average velocity first and then add the Stokes drift velocity to get the Lagrange average velocity [81, 82].

The calculation equations are:

Lagrange average velocity:

Euler average velocity:

Stokes drift velocity:

According to the above definition, the part for solving the Euler and Stokes average velocity is added to the model, and the Lagrange average velocity can be obtained by simply adding the simulation results of Euler and Stokes average velocity.

10) Dry and wet grid-processing technology

One of the most difficult problems in the numerical model of estuaries and coasts is to provide accurate simulation of the volume transport into or out of the tidal flat areas. This kind of volume transport may be several times larger than the volume transport of the main rivers. In the past 30 years, two methods have been used to solve this problem: one is the variable boundary method, and the other is the dry- and wet-point judgment method. In the first method, the calculation area is limited to a boundary between land and water. On such a boundary, the total water depth and vertical transport are equal to zero. Since this boundary moves in dry and wet areas, the model grid needs to be re-generated at each step. This method is suitable for ideal estuaries or coasts with simple shapes, but not for real estuaries or harbors with complex shapes including small ports, islands, obstacles and water bays. In the second method, the calculation area covers the maximum wet area. The numerical grid consists of dry and wet points and has a boundary defined by the boundary between land and water. The dry and wet points are distinguished with the local total water depth *D = H(x,y)* + *Q{x, y, t)* (of which, *H* is the reference water depth and £ is the surface relief). When *I)* > 0, it is wet point; when *I) =* 0, it is dry point. At the dry point, the velocity is automatically defined as zero, but the salinity value of the previous step is maintained at such point. Because this method is relatively simple, it is widely used in model of estuaries and coasts to simulate the volume transport of a waterbody in intertidal zone.

Whichever method is used to simulate the drying and wetting processes in estuaries or coastal intertidal zones, it must be ensured that mass conservation is observed. The dry-wet point in all these methods is determined by some empirical criteria, and the estimation of water transport in the dry-wet conversion area depends on: (1) criteria used to define the dry-wet point; (2) time step of numerical integration; (3) the horizontal and vertical resolution of the model grid; (4) amplitude of surface water level; (5) bottom t opography. In the *a* coordinate transformation model, it is also related to the thickness (*D _{m}i*

_{n}) of bottom viscous boundary layer.

A new dry and wet treatment method is used in FVCOM. This method has been proved to be effective in a series of tidal simulations in an ideal semi- closed harbor with an intertidal zone. The relationship between time step and grid resolution, externally forced amplitude, slope of dry and wet zone and thickness of bottom viscous boundary layer are discussed, and the criterion of selecting the time step is obtained. The criterion used in validation is mass conservation, which is a prerequisite for the objective assignment of dry-wet points in a bay.

The criteria for dry or wet judgement of nodes are as follows (Figure 4.4):

FIGURE 4.4: Definition of reference water depth (H), sea level (£) and topography height *(кв)*

For a triangular cell:

where *D _{m}i„* is thickness of the designated bottom viscous boundary layer, Iib is the height of the terrain connected to the boundary of the main waterways of the river, and

*i,j,k*are integers used to determine the three nodes of a triangular cell.

When a triangular cell is dried, the velocity at the mass center of the triangular grid cell is defined as zero, and there is no flux on the three boundaries of the triangular grid cell. This triangular grid cell will be removed from the flux calculation. For example, the integral form of the continuity equation in FVCOM is written as:

where *й* and *v* are the *x* and *у* component of vertical average velocity. In the dry-wet point system, only wet triangles are considered in the calculation, for the boundary flux of dry triangles is zero. In this way, the conservation of volume is ensured.

### Parameter setting

1) Parameter selection

The FVCOM model operation requires an initial temperature and salinity field, open boundary forcing and water depth data. The initial temperature and salinity field used for the calculation is the annual mean temperature and salinity field obtained from the historical observation data; the open boundary is forced by the water level forecasted by the tidal harmonic constant; the harmonic constants of six major semi-diurnal constituents and diurnal constituents *М**2**. S'**2*. A'2, *Кi. Оi* and Pi are used for forecast. The tidal harmonic constants used are obtained from the tidal model simulation of the Bohai Sea, Yellow Sea and East China Sea. The splicing of charts of different resolution provided by the Shandong Maritime Safety Administration is completed by using ArcGIS and other professional software. In addition, the related terrain and shoreline information extracted from the basic geographic information database of National Marine Information Center is supplemented, and the method of nearest-point interpolation of set thresholds is used to complete the production of high-precision terrain data products for the forecast area, thus it is ready for the configuration and development of the marine module of the forecast system.

2) Model operation

The sea areas calculated by the model are the Bohai Sea and Yellow Sea (37°-41°N, 117.42°-126.25°E). Grid refinement was carried out in key areas like Kiaochow Bay, coastal waters of Qingdao, Chengshan Cape routing waters and the Bohai Strait according to the contract requirements. The refined areas are divided into three levels. The primary refined area has a grid precision of 100 m and includes the following four sea areas: Kiaochow Bay, coastal waters of Qingdao, Chengshan Cape routing water and the Bohai Strait; the secondary refined area has a grid precision of 300 m and includes the following sea areas: the whole Bohai Sea and the sea areas west to 124°E; the rest of the sea areas belong to the tertiary refined area, which has a grid resolution of 1 bit. The resolution near the islands in the Bohai Strait reaches the order of 100 m, and basically all the islands are preserved and distinguished, avoiding the situation in the previous forecast that the marine current passes through the islands. The grid resolution in the non-key areas is reduced, which not only ensures the high-resolution simulation of the key areas, but controls the calculation cost and further optimizes the forecast efficiency of marine module.

The model adopts a coordinates, which are classified into 11a layers vertically and the upper-surface layer has higher resolution. The time step of the outer mode is 2 seconds and that of the inner mode is 12 seconds. The annual mean value obtained from the historical data is adopted as the initial value of temperature and salinity field. The actions of six constituents *Ml-* S'2. A2- *Оi*. *Кi* and *Pi* are considered in the model and introduced from the open boundary into the model through the water-level harmonic constant.

### Forecast results

The model automatically forecasts the conditions in the next 48 hours from 0:00 every morning and outputs the calculation result every hour. The output data is in the format of NetCDF, which is internationally accepted. The data format has a self-explanatory function, takes up small storage place, has good

FIGURE 4.5: Schematic diagram of researched sea area

FIGURE 4.6: Grid distribution of Bohai Sea and Yellow Sea (the gray represents the water depth)

FIGURE 4.7: Grid distribution of Bohai Sea

FIGURE 4.8: High-resolution grid distribution of Shandongtou versatility and is convenient for subsequent connection with OILMAP and visualized data processing.