Oil Spill Sensitivity Analysis of the Coastal Waters of Taiwan Using an Integrated Modelling Approach

Abstract

Pollution caused by marine oil spills can lead to persistent ecological disasters and severe social and economic damages. Numerical simulations are useful and essential tools for accurate decision making during emergencies and planning response actions. In this study, we applied the Princeton Ocean Model (POM) to determine current data, including seawater velocity, salinity, and temperature, and we obtained the fate and trajectory of spilled oil using OpenOil. Several probable oil slicks around Taiwan were simulated over time (12 months) and space (four spill locations in the marine area of each coastal city or county) using the model. The percentage risk under the effect of an oil spill is estimated. The risk zone of the coastal waters of Taiwan was identified based on the frequency of simulated oil slicks hitting the coast and sensitive resources. This information not only helps authorities guide the preparation of effective plans to minimise the impacts of oil spill incidents but could also be used to improve regulations related to shipping and vessel navigation in regional seas.

1. Introduction

Being an island nation, Taiwan began its industrial development in the 1960s. However, the lack of petrochemical energy sources in the country has resulted in the import of crude oil from oil-producing regions for the oil refining and petrochemical industries. Therefore, shipping transportation activities flourished. Owing to its unique geographical location, Taiwan is not only a crucial way-point for the West Coast of the United States and Asia Pacific routes for the international maritime transport system but also a transit centre for distant and near-sea shipping transport in East Asia, as shown in Figure 1. The increasing maritime traffic around Taiwan owing to recent economic growth in China and southeast Asia and the development of the greater East Asian trade system promotes industrial development, accompanied by frequent reports of accidents involving stranded ships or even sinking and collisions that may cause serious marine pollutions. The more booming the growth of technology in Taiwan, the more attention and awareness are being given to marine resources. To reduce the effects of pollution incidents on Taiwan’s aquatic ecosystem, implementing aggressive methods to prevent and control marine pollution and enhancing the prevention of marine pollution incidents and response capabilities is vital for the authorities. Assessments of oil slick risk not only help to better manage emergency plans, minimise economic loss [1,2], and mitigate the impacts of oil spills [3,4] but also support the formulation of emergency plans and the identification of mitigation policies [5]. However, for oil spills that occur in the water adjacent to Taiwan, most research has focused on using different numerical models to hinder particular events [6,7,8,9,10,11,12], analysing marine pollution management [13,14], and satellite observations [15]. Other efforts for better understanding the assessments of oil spills along the Eastern Mediterranean coasts can be found [16,17].

In this study, several probable oil slicks in Taiwan were simulated over time (12 months) and space (four spill locations in the marine area of each coastal city or county) using an integrated model. The integrated model included the Princeton Ocean Model (POM) [6,8,18] to determine three-dimensional current data and OpenOil (v1.1) [19] to obtain the fate and trajectory of the spilled oil. The percentage risks under the effect of an oil spill were then estimated. The risk zone of the coastal waters of Taiwan was identified based on the frequency of simulated oil slicks hitting the coast and sensitive resources.

Figure 1. The shipping transportation routes near Taiwan [20] (Google Map).

Figure 1. The shipping transportation routes near Taiwan [20] (Google Map).

2. Material and Method

2.1. Study Area

According to the Fisheries Agency of Taiwan, Taiwan’s Marine Protected Areas is “An area extending seaward from mean high tide line to a certain range, with special natural features, important cultural heritages and sustainable use of ecological resources, protected by law or other effective means”. The study area is placed along the coastline of the main island of Taiwan, which is about 1560 km (973.3 mi) long [21], covering 14 counties and cities, between 21°54′ and 25°18′ N latitude and 120°2′ and 122° E longitude, as illustrated in Figure 2. The authors used geographic information system (GIS) to create a boundary for the marine area within a 40 m radius extending seaward from the coast. Four sequential 10 km radial layers were created for the marine area of each county/city.

2.2. Methods

In this study, the risk magnitude of each zone under the effect of oil spills was determined based on the frequencies of modelled oil slicks hitting the coast and sensitive resources (i.e., coral reefs, marine/national parks, mangrove areas, scenic areas, artificial reefs, aquaculture, seagrass, wildlife refuges, sea turtles, sharks and dolphins). The distribution of these resources, as shown in Figure 3, is available online on the government’s public website, including the National Parks of Taiwan, the Fisheries Agency of Taiwan, ReefBase—The database of the Global Coral Reef Monitoring Network (GCRMN), and the International Coral Reef Action Network (ICRAN), or obtained from references [22]. The migratory routes of sea turtles and sharks were provided by Cheng [23] and Wang and Hsu et al. [24], respectively. The distribution of Taiwanese dolphins was based on that described by Wang et al. [25].

This work is based on two assumptions:

-

Coastal activities and sensitive resources are higher in number and diversity in shallow (closer to the shore) than in deep waters; therefore, the weights are different between layers.

-

Eleven sensitive resources have equal significance. Therefore, the weights are similar for different resources.

Percentage of risk caused by oil slicks hitting the coast and sensitive resources can be calculated as follows [26]:

$$PR=\left[\left({\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{w}_{ij}}}\times Nf{r}_{ij}\right)+\left({\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{w}_{ij}}}\times Nf{c}_{ij}\right)\right]\times \frac{100}{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n{w}_{ij}}}\times NFR{C}_{ij}}$$

(1)

where $w_{ij}$ are the weights given to an area in layer i =1, 2, 3, and 4 from a spill location in layer j = 1, 2, 3, and 4 in a city/county each month: 1, 0.5, 0.25, and 0.125, respectively; $Nf{r}_{ij}$ and $Nf{c}_{ij}$ are the normalised values (between 0 and 1) equal to the accumulative frequencies of simulated oil slicks hitting the coast (fc_ij) and sensitive resources (fr_ij) divided by the highest accumulative frequencies of simulated oil slicks hitting the coast (FC_ij) and sensitive resources (FR_ij) at one of the four layers in the city/county of interest monthly. $NFR{C}_{ij}$ is the normalised value (between 0 and 1) of FR_ij plus FC_ij.

2.3. Ocean Current Modelling

The Princeton Ocean Model (POM) [6,8,18,27] is a baroclinic model that uses a curvilinear orthogonal grid for the horizontal coordinate $\left(x,y\right)$ and a vertical grid for the sigma coordinate, $\sigma$. The governing equations of the models are

$$\frac{\mathsf{\partial }DU}{\mathsf{\partial }x}+\frac{\mathsf{\partial }DV}{\mathsf{\partial }y}+\frac{\mathsf{\partial }\omega }{\mathsf{\partial }\sigma }+\frac{\mathsf{\partial }\eta }{\mathsf{\partial }t}=0$$

(2)

$$\begin{array}{l}\frac{\mathsf{\partial }UD}{\mathsf{\partial }t}+\frac{\mathsf{\partial }{U}^{2}D}{\mathsf{\partial }x}+\frac{\mathsf{\partial }UVD}{\mathsf{\partial }y}+\frac{\mathsf{\partial }U\omega }{\mathsf{\partial }\sigma }-fVD+gD\frac{\mathsf{\partial }\eta }{\mathsf{\partial }x}\\ =\frac{\mathsf{\partial }}{\mathsf{\partial }\sigma }\left[\frac{K_M}{D}\frac{\mathsf{\partial }U}{\mathsf{\partial }\sigma }\right]+\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left[2A_{M}D\frac{\mathsf{\partial }U}{\mathsf{\partial }x}\right]+\frac{\mathsf{\partial }}{\mathsf{\partial }y}\left[A_{M}D\left(\frac{\mathsf{\partial }U}{\mathsf{\partial }y}+\frac{\mathsf{\partial }V}{\mathsf{\partial }x}\right)\right]\end{array}$$

(3)

$$\begin{array}{l}\frac{\mathsf{\partial }VD}{\mathsf{\partial }t}+\frac{\mathsf{\partial }UVD}{\mathsf{\partial }x}+\frac{\mathsf{\partial }{V}^{2}D}{\mathsf{\partial }y}+\frac{\mathsf{\partial }V\omega }{\mathsf{\partial }\sigma }+fUD+gD\frac{\mathsf{\partial }\eta }{\mathsf{\partial }y}\\ =\frac{\mathsf{\partial }}{\mathsf{\partial }\sigma }\left[\frac{K_M}{D}\frac{\mathsf{\partial }V}{\mathsf{\partial }\sigma }\right]+\frac{\mathsf{\partial }}{\mathsf{\partial }y}\left[2A_{M}D\frac{\mathsf{\partial }V}{\mathsf{\partial }y}\right]+\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left[A_{M}D\left(\frac{\mathsf{\partial }U}{\mathsf{\partial }y}+\frac{\mathsf{\partial }V}{\mathsf{\partial }x}\right)\right]\end{array}$$

(4)

where $U$ and $V$ are the horizontal velocity components; $\omega$ is the velocity component normal to sigma surfaces; $D$ is the total water depth; $\eta$ is the total water depth; $f$ is the Coriolis frequency; $g$ is the gravitational acceleration; $A_M$ is the horizontal diffusivity; and $K_M$ is turbulent closure parameter. Details can be found in Mellor [18].

In this study, the POM was used to simulate the ocean current around Taiwan with a horizontal grid resolution of 0.02° and 41 sigma layers over the entire year of 2018. The computational domain with bathymetry is illustrated in Figure 4.

The initial and boundary conditions’ driving forces of the surface elevation, temperature, and salinity are obtained from the Global Hybrid Coordinate Ocean Model [28], whose spatial and time resolution are about 1/12 degree and 1 day, respectively. And the atmosphere forcing is from the Weather Research and Forecasting model [29], which is operationally used to forecast the regional wind and weather by the Central Weather Bureau of Taiwan. In this study, the wind dataset covers a domain at 10° E to 35° E and 110° N to 137° N, and it has the spatial and temporal resolutions of 0.2 degree and 3 h, respectively. More detailed information can be found in the literature [6,7,8].

To validate the ocean currents, sea surface elevation was collected from six stations around Taiwan, and the data were compared with the sea surface elevation predicted by POM. These stations include Tamsui, Taichung Port, Donggang, Houbihu, Chenggong, and Hualien.

Table 1. Stations of validation.

No	Station	Longitude	Latitude
1	Tamsui	25°10′33″ N	121°25′29″ E
2	Taichung Port	24°17′16″ N	120°31′59″ E
3	Donggang	22°27′54″ N	120°26′18″ E
4	Houbihu	21°56′45″ N	120°44′43″ E
5	Chenggong	23°05′50″ N	121°22′49″ E
6	Hualien	23°58′50″ N	121°37′25″ E

and Figure 5 show the station positions.

2.4. Oil Spill Modelling

The oil slick is simulated using OpenOil [19]. In the model, the spilled oil drifts as a consequence of the combined action of wind, current and waves. Based on Zhang et al. [30], the drift velocity can be written as:

$$u={\alpha }_w\mathsf{{\rm Y}}{u}_w+{\alpha }_c{u}_c+u_{wave},v={\alpha }_w\mathsf{{\rm Y}}{v}_w+{\alpha }_c{v}_c+u_{wave},w=w_c$$

(5)

where $u_w$ and $v_w$ are the wind velocities at 10 m above the water surface; $u_c$, $v_c$ and $w_c$ are the surface water current velocities, which can be obtained from the POM; $u_{wave}$ is the wave-induced velocity, ${\alpha }_w$ is the wind drift factor, usually taken as 0.03 [31]; and ${\alpha }_c$ is the factor that accounts for the contribution of the drift of the oil slick on the water surface due to the current, and $\mathsf{{\rm Y}}$ is the transformation matrix of the wind angle [30].

For the rate of oil entrainment from the slick to the mixing layer, we refer to [9,32], which includes information on the relative importance of the inertial forces and oil–water interfacial tension and the ratio of viscous forces to the inertial and surface tension forces. In addition, the oil droplet size distribution can be expressed either as a number-based or a volume-based particle size distribution [33].

Other necessary data for the oil spill model are listed as follows:

Spill locations: For each county/city of interest, four locations with distances of 10, 20, 30, and 40 km to the coastline extending seaward at the boundaries of the four layers mentioned in Section 2.1. In this study, we used 56 spill locations, as illustrated in Figure 6.

Simulated periods: Five days/month/location during spring tide from January to December 2018.

Oil type and amount: Ten tons (10,000 particles) of heavy crude oil were used for every spill location.

Three-dimensional sea temperature and salinity: obtained from the POM model.

Wind: The magnitude and direction of wind during the simulated length of time 10 m above the land surface were obtained from the Global Forecast System (GFS) data of the National Centers for Environmental Prediction (NCEP).

The details of the three-dimensional hydrodynamic model (POM) combined with the OpenOil module were provided in the authors’ previous studies through the contents of the ocean circulation modelling and the hindcast of the oil spills [6,7,8,9].

3. Results

3.1. Tidal Current Validation

To validate ocean currents, we selected four months representing different seasons in Taiwan for each station. The time-history records of the water surface elevation simulated using POM were compared with data from the tide tables at four stations of the Water Resources Department of the Ministry of Economic Affairs from Figure 7, Figure 8, Figure 9 and Figure 10.

The calendar diagrams of the water elevation are shown in Figure 7, Figure 8, Figure 9 and Figure 10, where the red dotted line represents the numerical model simulation results, and the black solid line represents the actual measured values. The root mean square error (RMSE) used to estimate the simulation numerical error is defined as follows:

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _1^N{(X_{Obs}-X_{Sim})}^2}}$$

(6)

where $X_{obs}$ is the observed value, and $X_{sim}$ is the model simulation value. According to Equation (6), the RMSE values ranged from 0.048 to 0.063 m. In general, the elevations produced by the numerical model were in correlation with the observed elevations. The validations of the simulated current data from the POM model were intensively reported in previous studies [6,7,8,9] and therefore are neglected here.

3.2. Percentage Risk (PR) Caused by Oil Slicks Hitting the Coast and Sensitive Resources

The length of the coast and the number of sensitive resources dispersed by the simulated oil slicks from different release locations are shown in Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. The percentage of the risks of being influenced by oil spills is shown in Figure 20, the details of which will be discussed in the following sections.

In the northern part, including the marine area of New Taipei (Figure 11a), during winter (from October to March), the affected coast length (ACL) was mainly in the 10 and 20 m radii. Within a radius of 10 km, the value of the ACL varied from 35 km to approximately 120 km, and this fluctuation was 15–98 km within a radius of 20 km. The area had a radius of 30 km in October, the value of the ACL was approximately 40 km, and the value of the ACL in other months was insignificantly low. In contrast to winter, the ACL value was lower in summer (April to September), and it was no more than 22 km in a radius of 10 km (from a radius of 20 km onwards, almost unaffected). Compared to the monthly average ACL value in summer (9 km), the value in winter (approximately 72 km) was eight times larger. Regarding the number of sensitive resources, in winter, these values were 37 and 13, respectively, for two areas with radii of 10 and 20 km, respectively; in summer, this value was 15 (for a 10 km radius) and 1 (for a 20 km radius) (Figure 11b). It can be seen that the northern coastline was affected by the oil spill much more during winter (from October to March) than summer (from April to September). The surface currents in this area are approximately parallel to the coastline [34]. In winter, the northeasterly monsoon dominates the wind. Northerly wind forcing is one of the main reasons oil slicks hit the coastline.

The eastern part of the coast, including the Yilan, Hualien and Taitung areas, was the least affected by the oil spill. For Yilan, with a radius of 10 km (Figure 12a), the ACL value was highly affected in January, February, and June, with corresponding values of 54, 22, and 26 km, respectively. Its monthly average value is 9 km, which is approximately 8.2 times lower than that of the N. Taipei area (during winter). This is also evidenced by the number of sensitive resources: 26 and 4 for the 10 and 20 km radius areas, respectively, throughout the year. There was hardly any impact from the oil spill in the 30 and 40 km radius areas (Figure 12b).

In the Hualien area (Figure 13), throughout the year, the number of sensitive resources was seven and three for two areas with radii of 20 and 30 km, respectively. The calculation results showed that the impact of the oil spill only occurred in a radius of 10 m (with 22 sensitive resources), except for high ACL values in winter, falling in January, February, and October; the monthly ACL values in the Hualien area were low, less than 17 km. During winter, the average monthly ACL was approximately 6.3 km (10 times lower than that in the N. Taipei area).

In the Taitung area (Figure 14), compared to the areas at 10, 20, and 30 km, the number of sensitive resources in the 40 km area was the lowest, with a value of 5 for the whole year (occurring in January and February), whereas this number was 38, 26, and 23 for the whole year for the three areas closer to the shore (10, 20, and 30 km). However, compared to other regions (20 and 30 km radius areas), in the 10 km radius area, the monthly average ACL value was 18.5 km (nearly four times lower than that of the N. Taipei area). Notably, these impacts occurred during the first months of the year (January and February). Compared to other parts of Taiwan, the east coast, including the marine counties in Taitung, Hualien, and Yilan, had a low risk throughout the year. This may be due to the impact of a branch of the Kuroshio intrusion into the South China Sea, which dominates the area.

Meanwhile, in the southern region, including the coasts of Kaohsiung and Pingtung, almost all four areas (radii of 10, 20, 30, and 40 km) were affected by the oil spill. In the Kaohsiung area (Figure 15), this phenomenon occurred mostly in summer, from April to September (in contrast to the northern area, the N. Taipei Sea area, where the phenomenon mainly occurred in winter). During this period, the monthly average ACL was 43.5 km (for an area with a radius of 20 km), which was followed by 38.5 km (for a radius of 10 km) and approximately 30 km for the remaining two regions (radii of 30 and 40 km areas). In addition, the monthly average number of sensitive resources in all four areas (radius 10–40 km) ranged from 39 to 60, which showed that the oil spill phenomenon not only focused on one or two months but also affected all months throughout the summer.

In the Pingtung area (Figure 16), the impact of oil spills on coastal areas tended to gradually increase from February to March (winter), was strongest in June, July, and August (summer), and then gradually decreased. Simulations show that in July alone, although the number of sensitive resources for all four ranges (radii 10–40 km) was between 12 and 17, the ACL value was very large compared to other areas, ranging from 65 to 95 km. Therefore, the southern coast, including the coastlines of Kaohsiung and Pingtung, was affected by most of the simulated oil slicks dispersed from spill locations at different distances during the summer, whereas for the other counties/cities, this occurred in one or two months. This is due to the seasonal reversal of the monsoon wind, which shows opposite wind directions in summer and winter.

In the western region, from Taoyuan through Taichung to Tainan, we can see that in the central–western region, the impact of oil spills was greatest in the Taichung area (Figure 17). The number of sensitive resources was 46 in winter and 37 in summer, in which the ACL had a larger value in the 10 km radius region than in the 20, 30 and 40 km radius regions. January and February were the two months with the highest ACL values of 129 and 77 km, respectively (at a radius of 10 km), which were followed by a gradual decrease in the following months. The monthly average value in winter was higher than that in summer (46 km and 29 km, respectively).

Compared with the Taichung area, the number of sensitive resources was lower in Taoyuan, with values of 17 and 25 in winter and summer, respectively (Figure 18). Regarding the trend, similar to the Taoyuan area, the coastal area with a radius of 10 km was the most affected compared to areas with other radii. The ACL was highest at 71 km (in a 10 km radius area) in May; however, its monthly average value was not high: less than 10 km in winter and less than 30 km in summer.

The Tainan area (Figure 19) was the least affected area on the west coast. Coastal areas were only affected by oil spills in May. The values of the number of sensitive resources were 7, 14, and 14, corresponding to three areas with radii of 10, 20, and 30 km. Therefore, this phenomenon affected areas with radii of 10, 20, and 30 km, whereas areas with radii of 40 km were rarely affected. In May, the ACL values in these regions were 16, 86, and 51 km, respectively. Thus, for Tainan, during the summer, especially in May, two areas with radii of 20 and 30 km were affected by oil spills more than other areas and times of the year.

4. Discussion

Based on the analysis of the ACL and the number of affected sensitive resources (NASR), the cumulative frequencies of oil slicks hitting the coast were calculated. The calculation results showed that high risks of oil spills mainly appeared in the 10 m radius area for all observation areas around Taiwan Island during the year (Figure 20), except for the Hualien area, and the Tainan area. The Hualien area was not affected in June and December (Figure 20a), and the Tainan area was not affected in October and December (Figure 20c). Meanwhile, the Yilan area was only affected in (focused on) October, January, February, and March in winter and July in summer (Figure 20b).

Compared to the 10 km radius area, the 20 km radius area has a risk of intermittent oil spills throughout the year. The level of oil spill risk gradually decreased for areas with radii of 30 and 40 km. Within a radius of 40 km, the Hualien, Yilan, Tai-nan, N. Taipei, and Taoyuan areas had almost no risk of oil spills (Figure 20a–e).

Significant difference in oil spill risk can be observed between the southern and northern regions. While in the north (Figure 20d), the risk of oil spills only focuses on the 10 km area in both seasons, in Kaoshiung (south) (Figure 20f), the risk of oil spills affects all four areas (radii 10 to 40 km), and high risk occurs during the summer. In particular, the area with the highest risk was the area with a radius of 20 km, which was followed by areas with radii of 10, 30, and 40 km.

Another area in the south, Pingtung (Figure 20g), also had a high level of oil spill risk, similar to Kaoshiung, in the summer. The most notable was the high risk of oil spills occurring in May and July (within a radius of 20 km), which also occurred in areas with radii of 30 and 40 km (only in July). Another difference compared to all other areas is that in Pingtung, the risk level of oil spills not only occurs in summer but also in winter for all areas from a radius of 10 m to up to 40 km.

It should be noted that the difference in the simulated results among several marine continual/critical areas may be because of a variety of factors shaping the oil dispersion and fate patterns in those areas. Korotenko (2000) [35] highlighted that while the spill location, quantity and type of oil spilled, water currents, and local wind conditions influence oil dispersion and trajectory patterns in the sea, the meteorology of the area, properties of the spilled oil, and some physicochemical and biological processes, such as evaporation, emulsification, dissolution, oxidation, and biodegradation, affect the fate of oil slicks.

This information could support government officials, decision-makers, and local communities in limiting or avoiding oil-related activities in sensitive zones. More associated data related to the country’s coastline are needed, such as (1) the number of ports, frequency of oil spill incidents, number of ships, and coastal settlements; (2) local and international shipping routes; (3) oil seep locations; and (4) coastal settlements. This map allows for seasonally adjusted response measures and recommends that authorities alter traffic for shipping lanes that pass through sensitive areas.

The risk zone classification can be enhanced with more improvements to aid drawing up the sensitivity maps as depicted in Figure 21, Figure 22 and Figure 23 for the winter, summer, and whole year, respectively. In the figure, the significant accumulative frequencies are obtained by the average values of the larger one-third parts in Figure 20. In these figures, some conclusions can be drawn, as shown in the next section.

5. Conclusions

In this study, an integrated model based on POM for hydrodynamic simulation and OpenOil for oil spill prediction was applied to simulate oil slicks around the coastal waters of Taiwan over both space and time. Overall, high local risks occurred within a radius of 10 km in most counties/cities (except for Pingtung and Kaohsiung in some months). Because of the seasonal reversal of monsoon winds, while the northern coastline of Taiwan faces a much higher risk of oil spills during winter than summer, the southern coast was affected by most of the simulated oil slicks during summer. Compared to other parts of Taiwan, the east coast has a lower risk throughout the year. This information supports authorities in managing emergency plans, mitigating the impacts of oil spills, and improving shipping navigation regulations in regional seas. The sensitivity maps under the oil spill are drawn up using the significant frequency of simulated oil slicks hitting the coast and sensitive resources. More associated data related to the country’s coastline, such as the number of ports, frequency of oil spill accidents in the past, number of ships, coastal settlements, and local and international shipping routes, will be useful for future researchers.