A Modeling Study on the Oil Spill of M/V Marathassa in Vancouver Harbour

Abstract

The M/V Marathassa oil spill occurred on 8 April 2015 in the English Bay. In the present study, the trajectory and the transport mechanism of the spilled oil have been studied by using the three-dimensional and particle-based Oil Spill Contingency and Response (OSCAR) model forced by the Finite-Volume Community Ocean Model (FVCOM). FVCOM provided the hydrodynamic variables used by the oil spill model of OSCAR. The results showed that the fraction of the oil on the water surface and on the shoreline, as well as the amount of oil recovered were affected by the time of the initial release, the overall duration of the discharge, wind and recovery actions. The hindcast study of the M/V Marathassa oil spill showed that the likely starting time for the discharge was between 14:00 and 15:00, on 8 April 2015. The release may have lasted for a relatively long time (assumed to be 22 h in this study). The results of modeling in this study were found reasonably acceptable allowing for further application in risk assessment studies in the English Bay and Vancouver Harbour. The trajectory of the spill was mainly controlled by the tidal currents, which were strongly sensitive to the local coastline and topography of First Narrows and that in the central harbour. The model results also suggested that a high-resolution model, which was able to resolve abrupt changes in the coastlines and topography, was necessary to simulate the oil spill in the harbour.

1. Introduction

Canada has the world’s largest reserves of oil sands, which are deposits of bitumen in sand or porous rock [1]. The bitumen extracted from oil sands can be upgraded into various petroleum fuels (such as gasoline, diesel and aviation fuel) via proper hydro-treating processes. Due to the increasing bitumen and heavy oil production in Canada, the Trans Mountain Expansion Project (TMEP) was proposed to increase the capacity of bitumen and heavy oil transportation via pipeline from the province of Alberta, which has the majority of oil sands in Canada, to the west coastal province, British Columbia (BC). TMEP intends to triple the pipeline transportation capacity, which will consequently increase the oil tanker traffic by seven times on the BC coast, as well [2].

The biggest import and export port on the BC coast, Port of Metro Vancouver (PMV), consists of 34 anchorages (20 in the English Bay, 8 in Vancouver Harbour, 4 in the Indian Arm, 1 in Robert Bank, and 1 in Sand Heads), as shown in Figure 1. In 2017, PMV handled about 142 million tonnes of cargo, which is 5% more than the previous year (2016) [3]. This busy and growing vessel traffic in PMV increases the potential risk of oil spill. The Canadian Coast Guard (CCG) receives about 600 pollution reports on the BC coast every year, nearly 40 of which occur in the PMV [4]. For instance, a small oil spill took place in the English Bay (one of the PMV anchorages) on 8 April 2015, which resulted in at least 2800 L of oil released from the cargo vessel, M/V Marathassa [4].

The oil spill was first reported by the public at 16:48 Pacific Time on 8 April 2015 [6]. It was suspected that IFO-380 (Intermediate Fuel Oil 380) was spilled from the MV Marathassa vessel, which anchored at the location of latitude: 49°17.5167′ N, longitude: 123°11.2333′ W (Anchorage #12) [7]. During this spill event, several aerial overflights, including the National Aerial Surveillance Program (NASP) flights provided by Transport Canada, were conducted to estimate the pollutant on the water surface and shoreline as shown in Table S1. At 12:20 on 9 April, it was estimated that about 2800 L of spilled oil remained on the water surface [4]. This estimate did not include any weathered and previously recovered fuel oil [4]. It was estimated that the Western Canada Marine Response Corporation (WCMRC) recovered 1400 L of spilled oil by using three vessels with skimmer equipment (Table S2) [4]. However, the type and efficiency of the skimmers were not clearly recorded. Later on, the Shoreline Cleanup Assessment Technique (SCAT) teams surveyed over 85 km of shoreline between 9 April and 23 April 2015 and determined that the most contaminated shoreline was the west side of Stanley Park, North Vancouver, and West Vancouver [7]. On 14 April 2015, the City of Vancouver provided the spilled oil distribution map shown in Figure 2, which clearly showed the observed spilled oil on the water surface and the contamination on the shoreline [8]. Unfortunately, the specific cause of this spill was not clear, and the exact volume of spilled oil was unknown.

To understand the fate/trajectory of spilled oil in the marine environment, oil spill models may be used. Examples of this type of model includes: the SPILLCALC by Tetra Tech [9], the GNOME (General NOAA Operational Modeling Environment) model by NOAA (National Oceanic and Atmospheric Administration) [10], the OSCAR model by SINTEF (Stiftelsen for INdustriell of TEknisk Forskning ved NTH—Foundation for Industrial and Technical Research) [11], the OILMAP and SIMAP models by RPS-ASA (Applied Science Associates, Inc.) [12,13], the MOHID water model by MARETEC (Marine and Environmental Technology Research Center) [14] and the MIKE Oil Spill model by DHI (Dansk Hydraulisk Institut) [15]. For application to the TMEP, several organizations and consulting companies have simulated the potential risk of the oil spill in the Burrard Inlet, which geographically includes the English Bay, Vancouver Harbour, as well as in the Salish Sea (the mouth of Burrard Inlet opens onto the Salish Sea) by using various oil spill models. For example, the SPILLCALC model was used to simulate the possible trajectory of spilled diluted bitumen (dilbit) in 2013 [16]; the GNOME model was used to simulate the potential dilbit spill trajectory in 2015 [17]. However, these previously used models were limited by the following aspects: the stochastic model in SPILLCALC was 2D, which only tracked the surface transport of oil and did not provide the probability of water column contamination, and the study using the GNOME model simulated the trajectory of oil based on rough wind conditions and currents’ information, but not the fate/weathering processes.

Oil spill modeling typically incorporates the modeling of hydrodynamic forcing. H3D is a 3D hydrodynamic model that has been used in several studies of the oil spill in the Salish Sea and Burrard Inlet [16,18,19]. However, the resolution of this H3D model was relatively low in the study area, with a 1 km × 1 km horizontal grid space. In order to get a more accurate hydrodynamic forcing for the Salish Sea, the NEMO (Nucleus for European Modeling of the Ocean) model has been applied. The horizontal grid space of the NEMO model was almost uniform from 440 m × 440 m to 500 m × 500 m in the Salish Sea [20]. Unfortunately, this model was unable to simulate currents in the English Bay and Vancouver Harbour. Therefore, a high-resolution hydrodynamic model was needed for the modeling of the oil spill in the English Bay and Vancouver Harbour.

This study aims to validate a three-dimensional (3D) and high-resolution hydrodynamic model (the Finite-Volume Community Ocean Model (FVCOM)) as the first step. Then, the validated FVCOM output was incorporated into a 3D oil spill model (the Oil Spill Contingency and Response model (OSCAR)) to model the oil spill in the English Bay. Forty numerical simulations were carried out to test this coupled oil spill model based on historical information from the MV Marathassa oil spill. Specifically, the mass balance and trajectory of MV Marathassa spilled oil were simulated by varying different factors, including the oil start of release time, discharge duration, wind forcing and recovery action.

2. Materials and Methods

2.1. Hydrodynamic Forcing: FVCOM

2.1.1. FVCOM Description

The hydrodynamic forcing used for this study was generated using the Finite-Volume Community Ocean Model (FVCOM). It is a 3D, finite-volume and unstructured grid ocean model, which was first developed by Chen et al. [21] and further upgraded by joint efforts from researchers at the University of Massachusetts, Dartmouth and Woods Hole Oceanography Institution [22,23,24,25]. FVCOM allows the use of different resolutions to fit complex coastline and topography by using the triangle mesh system. The model used in the present study was based on the model set up by Wu et al. [26]. The model was capable of achieving relatively high resolution in the region of interest (English Bay and Vancouver Harbour in this case), as shown in Figure 3. For instance, the horizontal grid spacing is about 10 m in Vancouver Harbour and about 2 m around the bridge bases in the Second Narrows. The vertical grid has twenty-one sigma levels that were stretched gradually, in order to gain higher resolution in the surface and bottom layers. More detailed information of the model can be found in Wu et al. [26].

2.1.2. FVCOM Validation

The overall validation of the model has been done in Wu et al. [26] using tidal gauge water elevations and the ship-mounted Acoustic Doppler Current Profiler (ADCP) current data. Here, we further evaluate the model using surface drifter data, which were obtained from two Surface Current Tracker drifters (SCT). SCT is comprised mainly of wood for the structural support and cellulose sponge for floatation [27]. Four aluminum fins are mounted below the sponge to increase the surface area, and a zinc weight is installed at the very bottom of the unit to act as ballast [27]. There is also a thin aluminum disk installed above the cellulose sponge to facilitate labeling of the SCT with drifter ID and contact information [27]. SCT is a low-cost, low-impact, easily deployable drifter, tracks the surface currents and reports its location and timestamp [27]. Two SCT, named SCT1 and SCT2, were released in Vancouver Harbour (SCT1: 49°17.8812′ N, 122°57.6414′ W; SCT2: 49°17.8788′ N, 122°57.6432′ W) at 15:11 on 8 November 2015. The drifter’s locations and velocities were recorded every 2–6 min. It is notable that a time step of five minutes was applied during simulations.

The modeled trajectory was compared with the observed drifters’ trajectory. In addition, the prediction ability of FVCOM was statistically assessed by computing the following measures as shown in Equations (1)–(4) [28,29]: the Root-Mean-Square-Error (RMSE):

$$RMSE={\left\{\frac{1}{N} {{\displaystyle \sum }}_{i=1}^N{(X_{mod}-X_{obs})}^2\right\}}^{\frac{1}{2}}$$

(1)

the relative average error (E):

$$E=100\%\frac{{{\displaystyle \sum }}_{i=1}^N{(X_{mod}-X_{obs})}^2}{{{\displaystyle \sum }}_{i=1}^N({|X_{mod}-{\overline{X}}_{obs}|}^2+{|X_{obs}-{\overline{X}}_{obs}|}^2)}$$

(2)

the correlation coefficient:

$$R=\frac{{{\displaystyle \sum }}_{i=1}^N(X_{mod}-{\overline{X}}_{mod}) (X_{obs}-{\overline{X}}_{obs})}{{[{{\displaystyle \sum }}_{i=1}^N{(X_{mod}-{\overline{X}}_{mod})}^2 {{\displaystyle \sum }}_{i=1}^N{(X_{obs}-{\overline{X}}_{obs})}^2]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$$

(3)

and the quantitative agreement between model and observations:

$$Skill=1- \frac{{{\displaystyle \sum }}_{i=1}^N{|X_{mod}-X_{obs}|}^2}{{{\displaystyle \sum }}_{i=1}^N{(\left|X_{mod}-{\overline{X}}_{obs}\left|+\right|X_{obs}-{\overline{X}}_{obs}\right|)}^2}$$

(4)

where $X$ is the variable being compared with a time mean $\overline{X}$. The subscripts “mod” and “obs” represent the model results and observations, respectively.

After validating FVCOM, it was run for 10 days (from 5 to 15 April 2015) with a time step of 1 s and saved every 1 h to generate the hydrodynamic forcing, which did not include the waves, because the study period was reported as “very calm”, the surveillance photo showed no signs of breaking waves (white caps) and the non-breaking wave would also be very low due to low wind. The wave height used in the OSCAR model was computed from winds.

2.2. Oil Spill Model: OSCAR

The OSCAR model was used to simulate the mass balance and trajectories of the oil spill, based on the MV Marathassa oil spill’s observation data in the English Bay. This is a 3D particle-based model, which is designed based on SINTEF’s experimental field and laboratory data to support oil spill contingency and response decision making. The general structure of the OSCAR model is similar to most oil spill models as shown in Figure 4. The OSCAR model is capable of calculating the oil contamination on the sea surface and shorelines, in the water column and sediment, along with several oil weathering processes. Various oil weathering processes can be simulated by using the OSCAR model, including spreading, drifting, natural dispersion, chemical dispersion, evaporation, stranding, dissolution, adsorption, settling, emulsification and biodegradation of spilled oil. Overall, the OSCAR model has broad applications in oil spill modeling and has been validated in many related studies [30,31,32,33].

2.3. Wind Forcing: HRDPS Model

The present study used the wind forcing from the High-Resolution Deterministic Prediction System (HRDPS), which has been employed for weather prediction on the West Coast of Canada [34]. HRDPS is a set of the nested and Limited-Area Models (LAM) with forecast grids from the non-hydrostatic version of the Global Environmental Multiscale (GEM) model. This GEM has a 2.5-km horizontal grid spacing. The example of wind speed and direction at 16:00 on 8 April 2015 is shown in Figure 5. The dominant wind directions are south, southwest and southeast with speeds below 7 m/s near the release point from 5–12 April 2015, as shown in Figure 6.

2.4. Hindcast Study of the M/V Marathassa Oil Spill

2.4.1. Identification of Spilled Oil

PMV collected the polluted samples and identified that the spilled oil had an API gravity degree of 13 and a density of 978–979 kg/m³ (979 kg/m³ at 15 °C). The oil chemical compositions were tested as well, and the results showed that the spilled oil contained about 96–99% of bunker fuel. Further testing on oil physiochemical properties illustrated that the spilled oil had comparable physical and chemical properties as IFO-380 [4].

IFO-380 is typically classified as a heavy fuel oil with an API gravity of 10–17.1 degrees (density of 950–1000 kg/m³) [7,35]. It has a relatively high viscosity (maximum viscosity of 380 cSt [36,37,38]) and behaves as a semi-solid product at ambient temperature, which leads to a low rate of dispersion and evaporation [7,35]. The detailed chemical composition of IFO-380 was adapted from OSCAR’s oil database as shown in Table S3.

2.4.2. Potentially Influential Factors

As reported, the oil probably began to spill between 11:00 and 16:48 on 8 April 2015, but the exact start time of the release is still unknown. Five possible starting times (12:00, 13:00, 14:00, 15:00 and 16:00) were explored in this study. Although the wind forcing can be obtained via the HRDPS model as illustrated in Section 2.3, the wind speed was reported as quite low (<2.6 m/s) during 8–11 April 2015 [7]. It is therefore interesting to study the spilled oil fate and trajectory without taking the influence of the wind into consideration. Because of the lack of information on the duration of oil release and the lack of documentation on the details of recovery actions, the duration of discharge and recovery actions was included in the model study as two additional factors. The discharge duration was assumed as 2 h (a case of a relative instantaneous release) and 22 h (a case of slow release over a long period of time). The case with or without recovery actions was studied to investigate the impact on the fate and trajectory of the spilled oil. It is notable that the assumptions of oil recovery actions were made based on the CCG’s report, as shown in Table S4 [4] and the Western Canada Marine Response Corporation’s (WCMARC) website. A summary of the above-mentioned factors that might influence the fate and trajectory of spilled oil is presented in

Table 1. The studied factors and their corresponding settings.

Factor	Setting
Starting-releasing time	12:00	13:00	14:00		15:00	16:00
Wind forcing	With			Without
Discharge duration	2 h			22 h
Recovery action	Yes			No

, and detailed setup information for each simulation is shown in Table S5.

2.4.3. Deterministic Approach

The oil spill modeling can use both deterministic and stochastic approaches. A deterministic approach is used to simulate the fate and behaviour of oil from a single model run. This approach is helpful when studying a known historical oil spill event. A stochastic approach, on the other hand, is used to analyze the probability of oil contamination in the area of concern by overlaying a great number (tens to thousands) of individual deterministic simulations.

In this study, a deterministic approach was employed to study the mass balance and trajectories of the oil spill occurring on 8 April 2015. For each simulation, the oil was assumed to be released at Anchorage #12 (latitude: 49°17.5167′ N, longitude: 123°11.2333′ W) in the English Bay and then tracked for 3 days. A track duration of 3 days was used, because only a trace amount of spilled oil (5.9 L) remained on the water surface after 3 days, as reported by Transport Canada [4]. A time step of 20 min was selected to run the model. Since the hydrodynamic forcing was hourly, the use of a 20-min time step based on interpolation of current data helps to simulate a relatively smooth particle trajectory with less computation requirement compared with smaller steps (such as 1 min). The mass balances and trajectories for each individual simulation were saved every 1 h and represented by using 5000 particles. The chosen number of particles would affect the simulation to some extent. The use of 5000 is based on the preliminary test using 1000, 5000 and 10,000 particles. While the use of 1000 can produce a trajectory similar to that of 10,000, the use of a large number retains more details of the concentration field. Using 5000 can provide better details with less computational demand. This was also discussed in Reed and Hetland [39].

2.5. Statistical Analysis on Mass Balance

A full factorial design that incorporates the studied factors and their corresponding settings (Table 1) was generated by using Minitab software (version 18.1), resulting in 5 × 2 × 2 × 2 = 40 combinations in total. The mass balance (%) of oil calculated for the water surface, shoreline, water column, atmosphere, biodegradation and recovery was selected as the studied response. Analysis Of Variance (ANOVA) was carried out to evaluate the influence statistically of the studied factors on the mass balance. A p-value < 0.05 indicates that a certain factor has a significant influence on the mass balance. The normal distribution and constant variance on the error terms were assured during analysis, as well.

3. Results

3.1. FVCOM Validation

In order to validate the FVCOM, the simulated trajectory and velocities (U-velocity and V-velocity) from the model were compared with the observational data from the SCT drifters (SCT1 and SCT2). The simulated and observed trajectory are plotted in Figure 7 and Figure 8 for SCT1 and SCT2, respectively. It can be noticed that both SCT1 and SCT2 moved from east to west (Central Harbour to Second Narrow to Vancouver Harbour), and the modeled trajectory was comparable to the observed trajectories for both SCT1 and SCT2.

To quantify the prediction ability of FVCOM for velocities, statistical analysis, including the RMSE, relative average error (E), correlation coefficient (R) and skill, was carried out as presented in

Table 2. Results of statistical analysis between model simulations and observations.

Statistical Measures	SCT1		SCT2
	U-Velocity	V-Velocity	U-Velocity	V-Velocity
Root-Mean-Squared-Error (RMSE) (m/s)	0.149	0.056	0.158	0.052
Relative average error (E) (%)	48.65	76.38	73.42	68.44
Correlation coefficient (R)	0.618	0.256	0.310	0.383
p-value for R	0.000	0.001	0.003	0.000
Skill	0.719	0.514	0.577	0.551

. Both RMSE and E for FVCOM were satisfactorily low (less than 0.16 m/s and 77%, respectively) indicating that only a slight difference existed between modeled and observed velocities. The correlation between the modeled and observed velocities was represented by R values, and their significance levels were indicated by p-values. As shown in Table 2, all p-values were lower than 0.05, which again demonstrated the satisfied correlation between modeled and observed velocities. The skill values were all greater than 0.51, which further verified the agreement between modeled and observed velocities.

The time series velocities from simulations and observations for SCT1 and SCT2 are plotted in Figure 9 and Figure 10, respectively. In general, the simulated velocities matched well with that of observed data, even though some data were not recorded for unknown reasons. Overall, FVCOM was validated by using data from observed drifters, including trajectory and velocities in this study.

3.2. Impacts of Studied Factors on Oil Mass Balance and Trajectory

After validating FVCOM for hydrodynamic forcing, it was incorporated into the OSCAR model to study the mass balance and trajectory of the oil spilled from the M/V Marathassa. Four potential factors mentioned in Section 2.4.2 might influence the spilled oil mass balance, including the release start time, oil discharge duration, wind forcing and recovery actions. The raw data on their influence on the mass balance of oil (e.g., water surface, shoreline, water column, sediments, atmosphere, biodegraded and recovered) are presented in Table S6. Since the mass balance for the water column, sediments, atmosphere and biodegraded were all less than 3% due to the very weak wind/waves, only the oil components at the water surface, on the shoreline and the oil recovered were statistically analyzed in this study. Analysis of Variance (ANOVA) was carried out, and the p-values for the influence of studied factors on the oil mass balance are presented in

Table 3. The p-values for the influence of studied factors on oil mass balances. Significant influence (p-value < 0.05) is shown in bold.

Source	Water Surface	Shoreline	Recovered
Start-releasing time	0.000	0.001	0.008
Wind	0.000	0.000	0.196
Discharge duration	0.000	0.000	0.760
Recovery action	0.003	0.179	0.000

. The detailed mass balance distributions (after three days of tracking) are provided in Figure 11. In addition, the examples of trajectory comparison are shown in Figures S1–S4.

3.2.1. Influence of Release Start Time

From Table 3, it can be clearly seen that the oil start of release time had a significant impact on the mass balance of water surface, shoreline and recovered, as their p-values were less than 0.05. About 32.7% of spilled oil remained on the water surface and heavy contamination on the shoreline (63%) when the oil started spilling at 12:00, as shown in Figure 11a. In comparison with the 12:00 start of release time, the shoreline contamination was reduced (52.8%) along with an increased amount of spilled oil on the water surface (37.9%). If recovery was conducted, 7.46% of the oil was removed when it was released at 13:00. Interestingly, much more spilled oil (36.4%) can be recovered when the start of release time was 14:00, along with 32.9% contamination on the shoreline. The similar contamination on the shoreline was also observed if the oil spill started at 15:00 and/or 16:00, but larger amounts of spilled oil remained on the water surface (54.7% for and 70.5% for 16:00), rather than being recovered. In terms of oil trajectory, overall, the earlier the oil spill occurred, such as 12:00, the greater the contamination of the water surface and shoreline. However, this difference was not significant when the oil spill started at 14:00 and 15:00 (Figure S1).

3.2.2. Influence of Wind

The effect of the wind was to increase the amount of oil on the shoreline and decrease the amount of oil on the water surface compared to the no-wind simulation (Table 3 and Figure 11b). Specifically, the fraction of oil remaining on the surface decreased from 67% to 23%, and the amount of oil on the shore increased from 21% to 56%. This can be seen in the trajectory results (Figure S2), which illustrate the heavy contamination on the shoreline of West Vancouver and the western side of Stanley Park. The amount of oil recovered stayed roughly constant at 15–20% and was not influenced by wind forcing in this study.

3.2.3. Influence of Discharge Duration

A short discharge duration (2 h) led to more serious shoreline contamination (54.1% vs. 22.2%) than that of long discharge duration (22 h) and resulted in less oil on the water surface (29.6% for 2 h vs. 54.1% for 22 h). Most of the contaminant was still concentrated on the water surface around the release location after 3 days tracking when a long discharge duration was taken into consideration (Figure S3). The discharge duration did not play a significant role in the amount of oil recovered.

3.2.4. Influence of Recovery Action

Whether the recovery action did not significantly influence the mass balance of the oiled shoreline (Table 3), as well as the oil trajectory (Figure S4), only about 34.4% of the spilled oil remained on the water surface, if the recovery action removed 30.2% of the oil, and 33.5% ended up on the shoreline, as shown in Figure 11d. When no recovery action was taken (0% of recovered oil), 55.6% of spilled oil remained on the water surface, and 42.8% contaminated the shoreline.

4. Discussion

4.1. FVCOM Validation

In general, the simulated trajectory and velocities from the FVCOM were comparable with that of SCT drifters in this study. However, it was relatively less capable of predicting SCT2 U-velocity, as shown in Table 2. This is likely due to the following three reasons: (1) The SCT drifter used in this study was a shallow water drifter that worked close to the water surface. This type of shallow drifter was therefore susceptible to the surrounding windage, which could potentially cause higher uncertainty on recorded data. This was supported by a similar statement that was proposed by Halverson et al. [40] to explain the inconsistency of radial and observed velocities. (2) Relatively more observed data of SCT2 velocity were missed, which resulted in a less thorough comparison of modeled and observed data. (3) The difference between the model and the drifters may also be due to the winds and waves, which are not included in FVCOM.

4.2. Hindcast of the MV Marathassa Oil Spill

4.2.1. Comparison of Oil Trajectory

The model simulations of oil trajectory were evaluated and compared with the observed oil distributions, as shown previously in Section 1 (Table S1 and Figure 2). The oil distribution map indicated the observed oil trajectory on the water surface and the contamination on the shoreline in the English Bay and Vancouver Harbour from 8 April 2015 to 10 April 2015 [8]. The contaminated water surface area was labeled as 1–10, and the contaminated shoreline area was labeled as A–P. The comparison of the results of modeled and observed for water surface and shoreline contamination are listed in Tables S7 and S8, respectively. Four scenarios achieved the highest matches with the observation data. The studied factors’ setting in these four scenarios was: (1) oil started to release at 14:00, discharged continuously (22 h), with wind and without recovery actions (labeled as Scenario #4 in Table S5); (2) oil started to release at 14:00, discharged continuously (22 h), with wind and recovery actions (labeled as Scenario #8 in Table S5); (3) Scenario #4, which started to discharge at 15:00; (4) Scenario #8, which started to discharge at 15:00.

As described in Section 3.2.1 and Section 3.2.4, there was almost no difference between the oil trajectories whether or not the recovery action was used, and the difference of trajectories was not significant when oil started to discharge at 14:00 and 15:00. Therefore, as shown in

Table 4. Examples of water surface contaminant comparison. The simulated results were compared with observation data.

Start-Releasing Time	Scenarios #	Labels of Surface Contaminant										Matches (%)
		1	2	3	4	5	6	7	8	9	10
14:00	4	√	×	√	√	×	√	√	√	√	×	70
	8	√	×	√	√	×	√	√	√	√	×	70
15:00	4	×	√	√	√	√	×	√	√	√	×	70
	8	×	√	√	√	√	×	√	√	√	×	70

Scenario #4 represents oil discharged continuously (22 h), which then moves with the wind and without recovery actions; Scenario #8 represents oil discharged continuously (22 h), which then moves with wind and recovery actions. Detail factors’ setting in each scenario is shown in Table S5. “×” means the simulated results do not match with the observed data; “√” indicates the simulated results match the observed data.

and

Table 5. Examples of shoreline contaminant comparison. The simulated results were compared with observation data.

Time to Start Spill	Scenarios #	Labels of Shoreline Contaminant																Matches (%)
		A	B	C	D	E	F	G	H	I	J	K	L	M	N	O	P
14:00	4	×	√	√	√	√	×	√	√	√	√	×	√	×	√	×	×	62.5
	8	×	√	√	√	√	×	√	√	√	√	×	√	×	√	×	×	62.5
15:00	4	×	√	√	√	√	×	√	√	√	√	×	√	×	√	×	×	62.5
	8	×	√	√	√	√	×	√	√	√	√	×	√	×	√	×	×	62.5

Scenario #4 represents oil discharged continuously (22 h), which then moves with the wind and without recovery actions; Scenario #8 represents oil discharged continuously (22 h), which then moves with wind and recovery actions. Detail factors’ setting in each scenario is shown in Table S5. “×” means the simulated results do not match with the observed data; “√” indicates the simulated results match the observed data.

, those four scenarios mentioned above have achieved 70% and 62.5% matches in the comparison of surface contamination and shoreline contamination, respectively.

The oil trajectory in Scenario #4 that started at 14:00 is plotted in Figure 12. It can be seen that oil was first transported east of the oil release point and then moved to the southwest in the next twelve hours under the forcing of hydrodynamics and wind. Spilled oil was forced and moved into the First Narrow and eventually entered into Vancouver Harbour, which resulted in heavy oil contamination on the water surface and the shoreline around Vancouver Harbour and the First Narrow. There was no oiled shoreline until 19:00 (9 April 2015) when the oil reached English Bay Beach, which conformed well to the observed information [4,7]. The majority of shoreline contamination was on the west side of Stanley Park, West Vancouver, and North Vancouver, which matched the observation data well [4,7].

4.2.2. Comparison of Mass Balance

The oil mass balance in the simulations of the above-mentioned four scenarios was compared with that from the 2D, Automated Data Inquiry for Oil Spills (ADIOS, version 2.0) model (in CCG’s report) [7]. In CCG’s report, the ADIOS2 model was employed to study the mass balance of spilled IFO-380. Three metric tons (about 3067 L) of IFO-380 were assumed to be spilled at 4 Coordinated Universal Time (UTC) on 9 April 2015 (18:00 Pacific time on 8 April 2015) and then tracked for five days. A constant wind speed of 10 knots (about 5.14 m/s) was selected in CCG’s modeling [7].

The modeling results from the ADIOS2 model (Figure 13) indicated that approximately 11% and 2% of the oil was expected to evaporate and disperse, respectively after three days post spill. The other 87% of spilled oil was expected to remain on the water surface. This proportion of evaporation and remaining oil was totally different from OSCAR’s modeling. In Scenario #4, as shown in Figure 14 and Figure 15, around 1.4% of spilled oil was predicted to evaporate after three days of tracking. About 24.9% and 43.1% of the spilled oil were expected to contaminate the shoreline in Scenario #4 that started at 14:00 and 15:00, respectively. There are a number of reasons that could contribute to the different mass balance in the ADIOS2 and OSCAR models. The first main reason is that oil was trapped on shorelines as observed by oil spill responders. This process was included in the OSCAR model because it is a 3D fate/transport model that uses geographic and bathymetry data. By contrast, ADIOS2 is a weathering only model, and it has a limitation in accurately representing the significant onshore component. Another main difference is the evaporation, of which the rate is affected by wind, wave, currents and temperature [7]. The wind and currents conditions are very different in this study and CCG’s modeling.

By comparison, a lesser proportion of the spilled oil remained on the water surface in Scenario #8 due to recovery actions. Nearly 61.8% (1730 L) and 65.5% (1834 L) of the spilled oil were recovered when oil was discharged at 14:00 and 15:00, respectively. The modeled recovered oil was more than the actual volume of spilled oil recovered, which was probably due to the lack of information reported by the response vessels.

Overall, among the total number of forty studied scenarios, the results from four scenarios agree well with observations. The results indicate that the M/V Marathassa oil spill was most likely started between 14:00 and 15:00 on 8 April 2015. This spill was most likely a continuous slow release for an unknown period (assumed to be 22 h in this study) instead of an instantaneous release.

5. Conclusions

The FVCOM implementation for English Bay and Vancouver Harbour was further validated in this study by comparing the simulated trajectory and velocities with that of observed data from SCT drifters (SCT1 and SCT2). This validated FVCOM was then used to generate the hydrodynamic forcing in English Bay and Vancouver Harbour, which was input in the state-of-art OSCAR model to simulate the M/V Marathassa oil spill.

The M/V Marathassa oil spill event was numerically simulated to assess the ability of the coupled oil spill model. Forty scenarios were performed using the OSCAR model to study the effects of various input parameters on the fate and transportation of spilled oil. The results were compared with the available data of the M/V Marathassa oil spill. The trajectories from four scenarios match well with the observed data. The assumed recovery actions were performed better in the scenario of oil discharged continuously (22 h) with winds at 14:00 than that in the other simulations. The combined results of trajectory and mass balance indicated that the M/V Marathassa oil spill probably started between 14:00 and 15:00 (8 April 2015) and kept discharging oil for a relatively long time (assumed to be 22 h in this study). The weathering processes and movement of spilled oil and contamination distribution in the surrounding waters and coastlines were affected by wind and currents.

In general, the oil spill model integrating the OSCAR and FVCOMs has effectively simulated the offshore and onshore distributions of the M/V Marathassa oil spill. To our best knowledge, this is the first study that modeled the oil spill in the English Bay and Vancouver Harbour by using the OSCAR model.