Minimizing the Environmental Risk from Oil Tanker Grounding Accidents in the High North

Brice Assimizele^1,2,3, Robin T. Bye⁴

¹Department of Ocean Operations and Civil Engineering, NTNU—Norwegian University of Science and Technology, Ålesund, Norway.
²Department of Logistics, Molde University College, Molde, Norway.
³University Institute of Technology, Douala, Cameroon.
⁴Cyber-Physical Systems Laboratory, Department of ICT and Natural Sciences, NTNU—Norwegian University of Science and Technology, Ålesund, Norway.

DOI: 10.4236/ajor.2020.103005                

Abstract

To address the high environmental risk related to the increased oil tanker traffic in the High North, the Norwegian Coastal Administration (NCA) manages one of its vessel traffic service (VTS) centers in the town of Vardø, Norway. The fleet of tugboats, controlled by the VTS center operators, patrols the coastline to hook-up with any potential drifting oil tanker in the region of interest, before it runs ashore. Presently, the tugboats are controlled manually, which is not only challenging but less effective. In this paper, we develop two alternative binary integer programming models that give better tugboat policies in less computational time compared to previous work. Promising results with historical data illustrate great potential for optimal environmental risk reduction along the northern coast of Norway

Keywords

OR in Maritime Industry, Mixed Integer Programming, Search and Rescue, Oil Spill, Risk Management, Decision Support Systems, Dynamic Resource Allocation

1. Introduction

Maritime transportation plays an essential role in the international trade as it provides a cost-effective means to transport large cargo volumes. It is, however, characterized by a high level of uncertainty, which creates various risks in terms of fatalities, environmental pollution, and loss of property. In particular, oil spills from oil tankers grounding accidents have a devastating effect on the marine ecosystem [1], they involve prohibitive clean-up operations costs [2] and have a significant impact on the economic activities of the local communities [3]. Ship grounding accidents are generally caused by technical and mechanical failures, environmental factors and human errors. Presently, there is no consensus on the statistical distribution of the causes of shipping accidents [4], due to the different viewpoints of accident analysis. Thus, prevention remains the primary way of addressing the environmental issues related to maritime oil transport.

In pursuit of sustainable sea transportation in the High North, the Norwegian coastal administration (NCA) administers one of its vessel traffic service (VTS) centers in the town of Vardø, Norway. About 200 vessels are monitored daily by the VTS center of which five to six oil tankers receive special attention due to their size or risk of pollution. Through the automatic identification system (AIS), the VTS center obtains static information (cargo, identity, dimensions) and dynamic information (heading, position) from oil tankers moving in the region. Additionally, dynamic models of wind, ocean currents, wave heights and weather forecasts are used to predict potential drift trajectories and grounding locations of vessels. Moreover, the oil tankers are required by law to move along a predefined corridor approximately 50 km away from the coastline. Any oil tanker that losess its maneuverability through steering or propulsion failure is immediately assigned to the closest patrol tugboat for rescue operation before it runs ashore. The size of the zone of interest is about 1100 km of coastline, and the number of tankers entering the region makes it difficult to effectively move the tugboats at the right place in time.

Previous works [5] - [9] consider a one-dimension modeling approach and focus on the minimization of the distances between potentially drifting vessels and the nearest tugboat by means of genetic algorithms and a mixed integer programming (MIP). Their model and algorithms allocate tugboats to oil tankers, but do not give information on the probability of successful hook-up. In addition, the implementation of a one-dimension modeling approach is problematic, as it would give inaccurate geographical positions. Moreover, the dynamic risk model developed by [10] prioritizes oil tankers based on their potential oil spill volume in case of accidents, but does not suggest tugboat positions. [11] develops a two-dimension mathematical model, with hook-up probabilities, that minimizes the expected cost of grounding accidents. Despite these improvements, they do not account for uncertainty on weather conditions and use only one drift trajectory to predict the path followed by the potential drifting vessel. Moreover, the discretization of the region into cells of 5 by 5 km in their models is very large for optimal tugboat policies. In this paper, we address these issues by developing two alternative mathematical models that use more than one drift trajectory and smaller cells size for optimal decisions on tugboat positions in less computational time.

The remainder of this paper is organized as follows. Section 2 formulates the tugboat positioning problem and presents the two linear integer models that minimize the environmental risk from oil tanker grounding accidents. We discuss the integration of a receding horizon control (RHC) scheme into the mathematical models in Section 3. In addition, we present the numerical results with realistic test instances as well as case studies with historical events in Section 4. Finally, conclusions and further research are provided in Section 5.

2. Model Formulation

Following the formulation approach from [11], we discretize the time horizon into a finite set of time periods $\mathcal{T}=\left\{\mathrm{1,2,}\cdots \mathrm{,}T\right\}$ and subdivide the region of interest controlled by the VTS center into a finite number of cells $\mathcal{C}=\left\{\mathrm{1,}\cdots \mathrm{,}C\right\}$ . Each tugboat or oil tanker occupies one cell at each time period and can move to neighborhood cells depending on the speed, which is influenced by the weather conditions. The non-drifting oil tankers move on cells defined in the corridor and tugboats in the zone close to shore and approximately parallel to the corridor. Furthermore, the cells are constructed in a way that any tugboat will not need more than a time period to move from a given cell, except cells with very bad weather conditions at specified time periods.

For every oil tanker (vessels) v, in the set $\mathcal{V}$ , entering the region of interest, we consider independent potential drift trajectories at each time period in the defined time horizon. Let $\bar{\Omega }$ represent the set of possible scenarios in the planning horizon. Obviously, $\bar{\omega }\in \bar{\Omega }$ is a combination of drift trajectories (vessel scenarios), or normal routes in absence of an incident, followed by each vessel. That is, $\bar{\omega }=\left({\omega }_1\mathrm{,}\cdots \mathrm{,}{\omega }_v\right)$ , where ${\omega }_v\in {\Omega }_v$ denotes the vessel scenario for vessel v in a given time period and ${\Omega }_v$ is the set of all possible scenarios for vessel v. In case of drift, a vessel will follow a path denoted by $p=\left(c_1\mathrm{,}{c}_2\mathrm{,}\cdots \mathrm{,}{c}_T\right)$ , $c_t\in \mathcal{C}$ , which is a succession of cells followed by the drifting vessel. Although the model inputs are updated every time period as discussed in Section 3, uncertainty on drift trajectories is addressed by predicting more than one single potential path. That is, ${\omega }_v=\left(t\mathrm{,}{p}_i\right)\forall p_i\in {\mathcal{P}}_{{\omega }_{v}t}$ , where ${\mathcal{P}}_{{\omega }_{v}t}$ represents the set of all predicted paths for vessel scenario ${\omega }_v$ at time of distress call t and we denote by N the cardinality of ${\mathcal{P}}_{{\omega }_{v}t}$ . Thus, ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ represents the potential scenario for vessel v, where $t\in \mathcal{T}\cup \left\{T+1\right\}$ is the time the VTS center notices or is alerted to the distress of vessel v and $p_i$ are the predicted paths followed by the drifting vessel. In the absence of incident, t is set to $T+1$ .

Let $\mathcal{G}$ be the set of tugboats run by the VTS center in the town of Vardø. At the beginning of the planning, each tugboat $g\in \mathcal{G}$ is positioned at an initial cell $c_{0g}\in \mathcal{C}$ . The tugboats can only transit between neighborhood cells at each time period, which is determined by their maximal speeds in the planning horizon. Accordingly, let $\mathcal{F}\left(c\right)\subset \mathcal{C}$ be the set of cells that are adjacent to $c\in \mathcal{C}$ . Thus, $\mathcal{F}\left(c\right)$ represents cells that are reachable from cell c within one time period.

The main objective is to determine the position of tugboats at each time period such that the expected environmental consequence of oil tanker grounding accidents is minimized. Thus, let $K_{{\omega }_v}$ denote the environmental consequence associated with oil tanker v if vessel scenario ${\omega }_v$ occurs and no tugboat manage to rescue it before it runs ashore. In the next subsection, we present the risk

model $Ris{k}_{{\omega }_v}$ for any vessel scenario ${\omega }_v$ that helps to derive the risk for all

scenarios $\bar{\omega }$ , which represents the main function to be minimized in the two binary integer programming (BIP) models presented in Subsection 2.2.

2.1. Environmental Risk Modeling for Drift Grounding Vessels

A risk is a combination of the probability of an event and its consequence. In drift grounding accidents, the risk model for each potential vessel scenario ${\omega }_v$

is the product of the probability of failure $R_{{\omega }_v}$ , the probability of grounding given that it is adrift, $1-Q_{{\omega }_v}$ and the environmental consequence, $K_{{\omega }_v}$ :

$Ris{k}_{{\omega }_v}=R_{{\omega }_v}\left(1-Q_{{\omega }_v}\right)K_{{\omega }_v}$ (1)

An oil tanker might start drifting at any time period with a certain probability,

$R_{{\omega }_v}$ , that depends on the internal factors from the oil tanker itself as well as wind and current forces, and wave heights.

In Equation (1), $Q_{{\omega }_v}$ represents the probability of successful hook-up of the

drifting vessel with the nearest tugboat. The VTS center detects every drifting oil tanker and informs the nearest tugboat. Practically, the tugboat response time is determined by three main factors: 1) preparation time (reaction time and mobilization time), 2) sailing time and 3) connection or towing time. [10] illustrates the need for further analysis on the weather dependent towing time which is about 2 hours. Once the drifting vessel is reached by the tugboat, the time $t_l$ left before it runs ashore will determine the probability of successful hook-up. Thus, the probability of successful hook-up with the nearest tugboat given that

the vessel is adrift, denoted by $Q_{{\omega }_v}$ , mainly depends on $t_l$ . Accordingly, let $Q_{gc{\omega }_v}$ denote the probability of successful hook-up by tugboat g with drifting

vessel v, given tugboat g is in cell c at time of distress call t and vessel v follows scenario ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ . This probability depends on the position of the nearest tugboat at time of distress call, currents, wind, waves, distance of the vessel to shore and property of the drifting vessel such as type, draft, size and loading

condition. All these dependencies are captured in ${\lambda }_{gc{\omega }_v}$ , which is the predicted time left once tugboat g, in cell c at time of distress call t, reaches the drifting vessel in scenario ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ . As in [11], we determine ${\lambda }_{gc{\omega }_v}$ using the maximal

operational speed of the nearest tugboat and its location relative to the drifting vessel’s trajectory, and set

$Q_{gc{\omega }_v}=\frac{{\beta }_{{\omega }_v}\exp \left({\delta }_{{\omega }_v}\left({\lambda }_{gc{\omega }_v}-t_{\min }\right)\right)}{1+\exp \left({\delta }_{{\omega }_v}\left({\lambda }_{gc{\omega }_v}-t_{\min }\right)\right)}.$ (2)

The parameter $t_{\text{min}}$ represents the minimal remaining drift time required to attempt a hook-up. If ${\lambda }_{gc{\omega }_v}$ is less than $t_{\text{min}}$ , $Q_{gc{\omega }_v}$ is set to 0. In addition, ${\beta }_{{\omega }_v}\in \left[\mathrm{0,1}\right]$ and ${\delta }_{{\omega }_v}\ge 0$ represent the influence of weather conditions [11].

The environmental consequence of a drift grounding accident depends on the expected oil spill size (S) and the impact (I) of one tonne of oil on the environment

[10], such that $K_{{\omega }_v}=S_{{\omega }_v}{I}_{{\omega }_v}$ . It is important to note that spill size and spill impact

include both bunker and cargo spill. The spill size depends on the vessel type, size, loading condition and on whether the ship is single or double hulled. It is

found by combining the probability of an oil spill ${\tau }_{{\omega }_v}$ , given that the vessel run aground with the expected oil outflow in the event of oil spill, $O_{{\omega }_v}$ , in scenario ${\omega }_v$ : $S_{{\omega }_v}={\tau }_{{\omega }_v}{O}_{{\omega }_v}$ . Moreover, $O_{{\omega }_v}={\alpha }_{{\omega }_v}{\gamma }_{{\omega }_v}\text{Dwt}$ , where ${\alpha }_{{\omega }_v}$ is the expected outflow rate given as a percentage of the tank content volume and ${\gamma }_{{\omega }_v}$ is the

volume of cargo and bunker oil as a percentage of vessel dead-weight tonnage Dwt.

The oil spill impact per tonne depends on the type of oil spilled and the vulnerability of the affected area. This is modeled as environmental sensitivity index,

$E_{{\omega }_v}$ and oil type significance index, $L_{{\omega }_v}$ ( $I_{{\omega }_v}=E_{{\omega }_v}{L}_{{\omega }_v}$ ). The value $E_{{\omega }_v}$ depends on oil type and incorporates the vulnerability and ecological significance of the geographical area. In addition, $L_{{\omega }_v}$ describes the significance of the oil

type spilled. In case of drift for a given vessel scenario ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ , the impact of an oil spill will depend on the distance to shore, the weathering processes, the chemical composition of the oil, and the drift trajectory, which depends on the local wind and current condition.

2.2. Binary Integer Programming Models

We present two different binary integer models that minimize the expected environmental consequences from oil tanker grounding accidents. The first model, BIP-1, allocates the potential drifting vessels to the nearest tugboat, while the second model, BIP-2, focuses on the number of vessels that could not be rescued within a predefined threshold.

2.2.1. BIP-1 Model

We denote by $z_{gc{\omega }_v}$ a binary variable taking the value 1 if tugboat g is in cell c

and is the nearest tugboat at time of distress call t of vessel scenario ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ , and 0 otherwise. In addition, we assume that the probability of vessel scenarios ${\omega }_v$ is mutually independent. This assumption may not always be reasonable, however, we justify it by the fact that vessels in distress are usually spatially separated with few common environmental factors [11]. Thus, the probability for

a scenario $\bar{\omega }$ is given by $R_{\bar{\omega }}={\displaystyle {\prod }_{{\omega }_v\in \bar{\omega }}}{R}_{{\omega }_v}$ . In addition, we define $x_{gct}$ as a

binary variable taking the value 1 if tugboat g is in cell c at time t, and 0 otherwise. The environmental risk function to be minimized is then written as followed:

$f\left(\bar{z}\right)={\displaystyle \underset{\bar{\omega }=\left({\omega }_1\mathrm{,}\cdots \mathrm{,}{\omega }_v\right)\in \bar{\Omega }}{\sum }}{R}_{\bar{\omega }}{\displaystyle \underset{{\omega }_v\in \bar{\omega }}{\sum }}{\displaystyle \underset{g\in \mathcal{G}}{\sum }}{\displaystyle \underset{c\in \mathcal{C}}{\sum }}\left(1-Q_{gc{\omega }_v}\right)K_{{\omega }_v}{z}_{gc{\omega }_v}\mathrm{,}$ (3)

where $\bar{z}$ denotes the vector with components $z_{gc{\omega }_v}$ . The binary integer programming model below is developed to optimally minimize the objective function $f\left(\bar{z}\right)$ subject to some constraints.

Formulation

$minf\left(\bar{z}\right)$ (4)

s.t.

${\displaystyle \underset{c\in {\mathcal{F}}_{tg}\left(c^{\prime }\right)}{\sum }}{x}_{gct-1}\ge x_{g{c}^{\prime }t}\forall g\in \mathcal{G}\mathrm{,}\forall c^{\prime }\in \mathcal{C}\mathrm{,}\forall t\in \mathcal{T}\backslash \left\{0\right\}$ (5)

${\displaystyle \underset{c\in \mathcal{C}}{\sum }}{x}_{gct}=1\forall g\in \mathcal{G},\forall t\in \mathcal{T}$ (6)

$x_{g,c_{0g},0}=1\forall g\in \mathcal{G}$ (7)

$z_{gc{\omega }_v}\le x_{gct}\forall g\in \mathcal{G},\forall c\in \mathcal{C},{\omega }_v=\left(p,t\right)\in {\Omega }_v,v\in \mathcal{V}$ (8)

${\displaystyle \underset{g\in \mathcal{G},c\in \mathcal{C}}{\sum }}{z}_{gc{\omega }_v}=1\forall {\omega }_v\in \bar{\omega },v\in \mathcal{V}$ (9)

$x_{gct},z_{gc{\omega }_v}\in \left\{0,1\right\}\forall g\in \mathcal{G},\forall c\in \mathcal{C},\forall t\in \mathcal{T},\forall {\omega }_v\in {\Omega }_v,v\in \mathcal{V}$ (10)

Constraints (5) ensure tugboats move only between neighborhood cells. In addition, constraints (6) make sure tugboats are located in only one cell at each time period. The initial positions of tugboats are given in constraints (7), such that cell $c_{0g}$ is the position of tugboat g at the beginning of the time horizon. Constraints (8) and (9) allocate nearest tugboats to vessel scenarios and ensure that each vessel scenario ${\omega }_v$ is allocated to only one tugboat.

2.2.2. BIP-2 Model

The main objective of this model is to minimize the expected environmental consequence associated with the potential drifting vessel scenarios that could not be rescued within a predefined threshold $\rho$ . From the previous approach,

$Q_{gc{\omega }_v}=\frac{{\beta }_{{\omega }_v}exp\left({\delta }_{{\omega }_v}\left({\lambda }_{gc{\omega }_v}-t_{\text{min}}\right)\right)}{1+exp\left({\delta }_{{\omega }_v}\left({\lambda }_{gc{\omega }_v}-t_{\text{min}}\right)\right)}\mathrm{,}$ (11)

where ${\lambda }_{gc{\omega }_v}$ represents the estimated drift time left once the vessel is reached by the tugboat. Thus, we define $H_{gc{\omega }_v}$ as a binary parameter taking the value 1 if tugboat g is at cell c at time of distress call t and is not able to hook-up with vessel v under scenario ${\omega }_v=\left(t\mathrm{,}{p}_i\right)$ , within a predefined threshold time $\rho$ and 0 otherwise. Additionally, let $y_{{\omega }_v}$ be a variable that takes the value 1 if no tugboat

is able to hook-up with vessel v, doing scenario ${\omega }_v$ , within a predefined threshold $\rho$ and 0 otherwise. The expected environmental consequence to be minimized in then written as follow.

$g\left(\bar{x}\right)={\displaystyle \underset{\bar{\omega }=\left({\omega }_1\mathrm{,}\cdots \mathrm{,}{\omega }_V\right)\in \bar{\Omega }}{\sum }}{R}_{\bar{\omega }}{\displaystyle \underset{{\omega }_v\in {\Omega }_v\mathrm{,}v\in \mathcal{V}}{\sum }}{K}_{{\omega }_v}{y}_{{\omega }_v}$ (12)

where $\bar{x}$ denotes a vector with components $x_{gct}$ .

In addition to $g\left(\bar{x}\right)$ , we define $u\left(\bar{x}\right)$ as a function that gives incentive to tugboats to optimally position themselves once the threshold is reached:

$u\left(\bar{x}\right)={\displaystyle \underset{\bar{\omega }=\left({\omega }_1\mathrm{,}\cdots \mathrm{,}{\omega }_V\right)\in \bar{\Omega }}{\sum }}{\displaystyle \underset{{\omega }_v\in {\Omega }_v\mathrm{,}v\in \mathcal{V}}{\sum }}{\displaystyle \underset{c\in \mathcal{C}}{\sum }}{\displaystyle \underset{g\in \mathcal{G}}{\sum }}\frac{K_{{\omega }_v}{x}_{gct}}{{\lambda }_{gc{\omega }_v}}\mathrm{.}$ (13)

The BIP-2 model that minimizes the objective functions $g\left(\bar{x}\right)$ and $u\left(\bar{x}\right)$ is presented below.

Formulation

$min\mathrm{<mi>\; </mi>}g\left(\bar{x}\right)+u\left(\bar{x}\right)$ (14)

s.t.

constraints (5)-(7) in BIP-1 and

${\displaystyle \underset{c\in \mathcal{C}}{\sum }}{\displaystyle \underset{g\in \mathcal{G}}{\sum }}{H}_{gc{\omega }_v}{x}_{gct}\le y_{{\omega }_v}+card\left(\mathcal{G}\right)-1\forall {\omega }_v\in \bar{\Omega }$ (15)

$x_{gct}\mathrm{,}{y}_{{\omega }_v}\in \left\{\mathrm{0,1}\right\}\forall g\in \mathcal{G}\mathrm{,}\forall c\in \mathcal{C}\mathrm{,}\forall t\in \mathcal{T}\mathrm{,}\forall {\omega }_v\in {\Omega }_v\mathrm{,}v\in \mathcal{V}$ (16)

Constraints (15) capture the vessels that could not be reached within the predefined threshold and the other constraints are the same as in BIP-1.

3. Integrating the BIP Models with the RHC

In this section, we integrate the RHC algorithm with the BIP models to account for uncertainty in weather conditions and dynamic changes of the input parameters. A RHC is a class of algorithms that make use of explicit process models to predict future response of a system, with optimizations as intermediate steps. The main idea is to dynamically run the BIP model in real time, while implementing only the first time period over the whole planing horizon (see [11] and references therein). Indeed, updating the parameters with new accurate values improves the output quality of the BIP model. For instance, information about vessels entering and leaving the region of interest as well as available operational tugboats need to be updated at every time period. Moreover, it is certainly efficient to pro-actively include next time periods in the BIP model. We then implement the RHC algorithm as present in Table 1.

4. Test Cases

We present the numerical settings in this section and discuss the quality and performance of the BIP models, compared with previous work, run with realistic test cases. In addition, the promising results with a historical event highlight the important features derived from the BIP models as a decision support tool for the NCA managers.

4.1. Computational Settings

The region of interest covers about 1100 km of coastline and the corridor is on average 50 km away from the coast. We discretize the region by collecting the

Table 1. RHC-BIP algorithm.

center position of each cell and transform them into Cartesian coordinates for input to the model. Once the optimal solution of the BIP model is obtained, the drift trajectories, oil tanker and tugboat positions are transformed back to geographical coordinates. The drift trajectories are obtained using the AIS and Norwegian Meteorological Institute (NMI) information with the algorithm presented in Section 3 (see [11] for details). Presently, the VTS center operates a fleet of two tugboats with an average operating speed of 12 knots and vessels typically have an operating speed of 14 - 15 knots. The VTS center subdivided the region into two zones, where each zone is assigned to one tugboat. The first zone spans from the border to Russia to Torsvåg and the second zone from Torsvåg to Røst. Moreover, the hook-up probabilities are computed using the formula presented in Section 2.1. We set the threshold $\rho =5$ hours, which is larger than the 2 hours average towing time. This threshold value for BIP-2 can be changed according to the VTS center operators needs as discussed in Section 4.2.4.

Previous research, on the same region of interest, conducted by [10] presents a resource specific environmental sensitivity index ranging from 1 to 9. The coastal segments used for summarizing index values are presented in Figure 1. Because of data accessibility, we randomly generate the environmental consequence according to a uniform random variable on $\left[\mathrm{1,9}\right]$ times an absolute normal random variable with a mean value of 5 and a standard deviation of 2. This is a reasonable assumption as the main goal is to dynamically assess tugboat positions according to different consequence levels. Moreover, the only unavailable input data required to compute the environmental consequence is the Deadweight Tonnage (Dwt) of each vessel, which is accessible to the VTS center operators for practical implementation of the BIP model. The failure probability is randomly chosen between $\left[\mathrm{0.01,0.09}\right]$ as in [11]. Specific settings to each case are presented in the corresponding subsections.

Figure 1. Coastal segments used for summarizing index values, shown in color codes on 10 by 10 km squares [10].

All computations are carried out on a personal computer with an Intel® Pentium® IV 3.0 CPU and 4.0 GB of RAM. The optimization software Gurobi 6.0.5 is used as a solver, with Python 2.7.3 and Pyomo 4.2, on Microsoft Windows 7.

4.2. Test Cases with Realistic Data

This subsection discusses the numerical results for three different cases. For each case, the models are run for a total of 24 hours with 100 different instances and the environmental risk associated with each scenario is computed according to the tugboat positions from each model policy. A total number of 6 oil tankers, which correspond to the current average daily number, are used with random geographical positions, directions and speeds.

4.2.1. Case 1: Large Cell Size

In this case, we use a large cell size of 5 by 5 km as in [11] and compare the quality and performance of their MIP-U model with those of BIP-1 and BIP-2 models. In order to allow comparison with the MIP-U model from previous work, we only use one predicted drift trajectory for vessel scenario ${\omega }_v$ . The numerical results in Table 2 present the statistical values of the computational time and environmental risk related to each scenario from the test instances for each model.

As presented in Table 2, the average environmental risk is almost the same for the three models with a value of ≈11.8. In addition, the standard deviations from the BIP-1 and BIP-2 models are slightly higher than that of MIP-U. The computational time is, however, very high for the MIP-U model. In fact, the BIP-1 is about three times faster than MIP-U model. Moreover, the BIP-2 performance is by far better than those of MIP-U and BIP-1. Nevertheless, the performances of these three models are acceptable for this case with small scenarios number and

Table 2. Numerical results for 5 by 5 km cells size.

large cells size. Essentially, the models are run dynamically, where only the first step of one hour is implemented, as described in Table 1, and each of these models can be run every hour. The very small standard deviation of the computational time in BIP-2 model is a good indication of its great performance when considering lager scenarios number discussed in the next subsection.

4.2.2. Case 2: Smaller Cell Size

In this subsection, we consider smaller cells size, of 2 by 2 km, compared to that of Case 1. We also use a single path as in Case 1 to predict the drift trajectory of each potential drifting vessel. The computational results for MIP-U and BIP-2 models are presented in Table 3. The BIP-1 model is not included because it could not yield an optimal solution after two hours of run time with each test instance of this case.

Noticeably, the computational time for both MIP-U and BIP-2 models have considerably increased compared to the values in Case 1. Indeed, smaller cells size increases the overall number of cells, which consequently expand the problem size. The average run-time for the BIP-2 model is equal to 5.5 minutes with a standard deviation of only 3.7 minutes. These values are significantly smaller than those of MIP-U model. The average performance of almost 40 minutes, with a maximum of 91.5 minutes, in the MIP-U model makes it impossible to be run dynamically and account for uncertainty with the algorithm in Table 1. Thus, the BIP-2 model is well suited for large number of cells as well as possible extension of the current region of interest. Additionally, the environmental risk for this case with small cells size has considerably decreased compared to that of Case 1. In fact, a large number of cells increase the flexibility of tugboats and allow for better positions that minimize the environmental risk.

4.2.3. Case 3: Large Scenarios Number

This case study uses the same cells size of 2 by 2 km as in Case 2. The main difference, however, is the total number of scenarios. In this case, we use three drift trajectories to predict the path followed by each vessel scenario ${\omega }_v$ . In addition, we consider up to two possible vessel scenarios ${\omega }_v$ for each scenario $\bar{\omega }$ . This gives more than 50,000 total number of scenarios.

For this case, none of the MIP-U and BIP-1 models are able to provide solutions in less than three hours. Thus, the numerical results presented in Table 4

Table 3. Numerical results for 2 by 2 km cells size.

Table 4. Results for 2 by 2 km cells size with BIP-2 model.

are those of the BIP-2 model only. The initial risk column in Table 4 represents statistical values of the potential environmental risk if no action is taken from tugboats. Fortunately, the computational time is less than 15 minutes for all the instances considered in this case, which makes it possible to combine the BIP-2 model with the receding horizon control algorithm described in Table 1. Obviously, the increase in the computational time compared to Case 2 is due to the larger number of scenarios and cells, which also increase the complexity of the model. Additionally, the average and standard deviation of the risk have increased from 2.253 to 3.769 and 5.070 to 5.762, respectively (see Table 3 and Table 4). This could be misleading, but it is important to note that in Case 2, only one drift trajectory is used to predict the path followed by a drifting vessel, which does not account for uncertainty such as wave heights, ocean currents and wind forces, and underestimates the real potential risk.

4.2.4. Case 4: Different Values of the Threshold

In order to assess the effect of the threshold parameter on the solution values, we run the BIP-2 model with a real world instance. The test case consists of 9 vessels that moved along the coast over a time period of 15 hours. In addition, we discretize the region of interest into small cells with large scenarios number as in Case 3. The expected potential risk for each value of the threshold $\rho \in \left\{\mathrm{2,5,8,11}\right\}$ is presented in Table 5.

As presented in Table 5, the average potential environmental risk increases with higher values of the threshold. For $\rho =2$ the expected risk is equal to 4.39, whereas that of $\rho =11$ is equal to 4.92. We notice, however, that the standard deviation of the risk decreases with higher values of the threshold. For $\rho =2$ the standard deviation of the risk is equal to 4.62, whereas that of $\rho =11$ is equal to 4.12. Moreover, the maximum value of the risk decreases from 24.85 to

Table 5. Results for different values of the threshold with BIP-2 model.

18.71 while the minimum risk value increase from 1.04 to 3.11 for a threshold value of 2 and 11, respectively. The threshold parameter gives more options with regards to tugboats policy and level of risk. That is, the managers at the VTS center will have to make a trade-off between having smaller standard deviation of the risk and avoiding the high risk of worst case scenarios at the expense of higher expected potential risk.

The histograms in Figure 2 present the distribution of the potential environmental risk for different decisions that are computed using different values of $\rho$ . As shown in Figures 2(a)-(c), the tails of the distribution reduces with higher values of $\rho$ . This is mostly seen on the upper tails, which represent the very rare but high risk scenarios. The policy for $\rho =11$ in Figure 2(c) considerably shapes the distribution of the risk by reducing the maximum, but also decreases the probability of smaller risk values.

4.3. Test Case with Historical Data

This case is based on real-world data collected from the AIS and the NMI. In addition, we use the basemap library in python to plot and draw the maps with drift trajectories, oil tanker and tugboat positions.

On the 21st of Mars 2014 at 11:10 pm, a vessel ran aground at N71˚01.06'N - 028˚27.46'E after 15 hours of drifting time. At the time of distress call, the nearest tugboat was located at N70˚40'N - 023˚40'E and unsuccessfully tried to reach the drifting vessel. The tugboat was located about 142.8 km away from the vessel at the time of grounding. We ran the BIP-2 model for 15 hours prior to the time of distress call for more than 50,000 scenarios. A total number of 7 vessels, including the one that ran ashore, sailed along the region during the considered planning horizon. Their corresponding directions, latitudes, longitudes and speeds over ground (SOG) at the beginning of the planning horizon are presented in Table 6.

The results for the first and last time period are presented in Figure 3, where the initial positions of vessels are presented in red cycles. In addition, the drift trajectories for vessel scenarios that could be rescued within the threshold of $\rho =5$ hours are in blue solid lines while those of the vessels that could not hook-up with tugboats within the threshold are in red solid line. Moreover, the actual drift trajectory followed by the grounded vessel is represented by green solid lines. Furthermore, the two red directed lines in Figure 3(b) are the actual

Table 6. Case 3A. Initial vessel directions, positions and SOG.

(a)(b)(c)

Figure 2. Case 4. Distribution of the risk for different value of $\rho$ . The histograms show the reduction of the tail of the distribution as $\rho$ increases. (a) $\rho =2$ vs. $\rho =5$ ; (b) $\rho =2$ vs. $\rho =8$ ; (c) $\rho =2$ vs. $\rho =11$ .

Figure 3. Results of the first and last time period. The dashed green lines represent the suggested movements of tugboats by the BIP-2 model and the predicted drift trajectories that can be rescued within the threshold $\rho =5$ hours are in blue solid lines while those in red solid line represent the drift trajectories of the vessels that could not be hook-up within the threshold. Additionally, the actual drift trajectory of the vessel that ran aground is represented by green solid lines. The two directed red lines close to shore in (b) are the actual positions of tugboats from the time of distress call to the time of grounding. (a) t = 1; (b) t = 15.

paths followed by the tugboats from the time of distress call to the time of grounding while the green solid lines linked with small squares represent the suggested movement of the tugboats by the BIP-2 model. The risk values associated with each vessel scenario is not presented in the figure because of the small visibility.

A zoomed-in view of the grounded location in Figure 4 shows actual and predicted drift trajectories of the grounded vessel. Additionally, the small square in red represents the suggested position of the nearest tugboat by the BIP-2 model at the time of distress call. The probability of successful rescue of the grounded vessel by the nearest tugboat with the BIP-2 model is about 0.70 while that of the actual tugboat policy is equal to 0.2. That is, the grounded vessel had 70% chance to be rescued if the BIP-2 model was implemented. For this particular case, the hook-up probability is slightly smaller than that of the MIP-U model of 0.86 in [11]. However, their model does not account for the uncertainty on drift trajectory. Thus, in the long run, more vessels in distress will have very low probability of successful hook-up with tugboats as very few number of scenarios are considered.

5. Conclusions

In this paper, we address the environmental risk related to oil tankers traffic along the northern Norwegian coast. We propose two alternative models that could be used as a decision support tool at the vessel traffic service center in the town of Vardø, for a better rescue operation of vessels in distress. These models are combined with a receding horizon control algorithm to account for uncertainty in weather conditions and to dynamically update the constantly changing input parameters. For a large cells size of 5 by 5 km and smaller scenarios number,

Figure 4. Zoom on the grounded vessel. The predicted drift trajectories are represented in blue solid lines and the actual drift trajectory of the vessel that ran aground is represented in green solid lines. In addition, the red directed path represents the positions of the nearest tugboat from the time of distress call to the time of grounding. The dashed green lines are the suggested movements of the nearest tugboat by the BIP-2 model.

the BIP-1 and BIP-2 models outperform the MIP-U model from previous work. In addition, the BIP-A model is by far faster than the other models for large scenarios number and small cells size, which considerably adds complexity to the models. Moreover, the BIP-2 model gives flexibility to the operators at the VTS center by allowing different threshold levels. The results with a historical event indicate better decisions on tugboat patrol operations.

It is recommended that further research is done to determine the optimal fleet of tugboat required as well as extension of the BIP models to consider other search and rescue operations. Additionally, more research is needed to better assess the failure probabilities of vessels, oil spill rates, probability of oil spill given that an accident has occurred, and environmental consequence of the region of interest. Furthermore, the hook-up probability formula could be better estimated with empirical data including new features such as “ship arrestors” newly acquired by the NCA to reduce the speed of the drifting vessels.

Acknowledgements

This work was carried out as part of the project Dynamic Resource Allocation with Maritime Application (DRAMA), grant no. ES504913, with partial funding from the Norwegian research Council and NTNU—University of Science and Technology, Ålesund, Norway. We sincerely acknowledge Trond Ski and colleagues at the Norwegian Coastal Administration for their support and input during the work with this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Lecklin, T., Ryömä, R. and Kuikka, S. (2011) A Bayesian Network for Analyzing Biological Acute and Long-Term Impacts of an Oil Spill in the Gulf of Finland. Marine Pollution Bulletin, 62, 2822-2835. https://doi.org/10.1016/j.marpolbul.2011.08.045
[2]	Montewka, J., Weckström, M. and Kujala, P. (2013) A Probabilistic Model Estimating Oil Spill Clean-up Costs—A Case Study for the Gulf of Finland. Marine Pollution Bulletin, 76, 61-71. https://doi.org/10.1016/j.marpolbul.2013.09.031
[3]	Crotts, J.C. and Mazanec, J.A. (2013) Diagnosing the Impact of an Event on Hotel Demand: The Case of the BP Oil Spill. Tourism Management Perspectives, 8, 60-67. https://doi.org/10.1016/j.tmp.2013.07.002
[4]	Goerlandt, F. and Montewka, J. (2014) A Probabilistic Model for Accidental Cargo Oil Outflow from Product Tankers in a Ship-Ship Collision. Marine Pollution Bulletin, 79, 130-144. https://doi.org/10.1016/j.marpolbul.2013.12.026
[5]	Bye, R.T., van Albada, S.B. and Yndestad, H. (2010) A Receding Horizon Genetic Algorithm for Dynamic Multi-Target Assignment and Tracking: A Case Study on the Optimal Positioning of Tug Vessels along the Northern Norwegian Coast. Proceedings of the International Conference on Evolutionary Computation (ICEC), Vol. 1, 114-125.
[6]	Bye, R.T. (2012) A Receding Horizon Genetic Algorithm for Dynamic Resource Allocation: A Case Study on Optimal Positioning of Tugs. In: Madani, K., Dourado, C.A., Rosa, A., Filipe, J., Eds., Computational Intelligence, Springer, Berlin, Heidelberg, Vol. 399, 131-147. https://doi.org/10.1007/978-3-642-27534-0_9
[7]	Assimizele, B., Oppen, J. and Bye, R.T. (2013) A Sustainable Model for Optimal Dynamic Allocation of Patrol Tugs to Oil Tankers. Proceedings of the 27th European Conference on Modeling and Simulation, Aalesund, 27-30 May 2013, 801-807. https://doi.org/10.7148/2013-0801
[8]	Bye, R.T. and Schaathun, H.G. (2014) An Improved Receding Horizon Genetic Algorithm for the Tug Fleet Optimisation Problem. Proceedings of the 28th European Conference on Modelling and Simulation (ECMS 2014), Brescia, 27-30 May 2014, 682-690. https://doi.org/10.7148/2014-0682
[9]	Bye, R.T. and Schaathun, H.G. (2015) Evaluation Heuristics for Tug Fleet Optimisation Algorithms: A Computational Simulation Study of a Receding Horizon Genetic Algorithm. Proceedings of the 27th European Conference on Modeling and Simulation, Alesund, 27-30 May 2015, 270-282.
[10]	Eide, M.S., Endresen, O., Breivik, O., Brude, O.W., Ellingsen, I.H., Roang, K., Hauge, J. and Brett, P.O. (2007) Prevention of Oil Spill From Shipping by Modelling of Dynamic Risk. Marine Pollution Bulletin, 54, 1619-1633. https://doi.org/10.1016/j.marpolbul.2007.06.013
[11]	Assimizele, B., Royset, J.O., Bye, R.T. and Oppen, J. (2018) Preventing Environmental Disasters from Grounding Accidents: A Case Study of Tugboat Positioning along the Norwegian Coast. Journal of the Operational Research Society, 69, 1173-1792. https://doi.org/10.1080/01605682.2017.1409157