Simulating the Behaviour and Fate of an Oil Spill Using a Coupled Three-Dimensional Hydrodynamic Model

Oil spill modelling is a crucial tool in the permitting process for assessing the potential consequences of an oil spill on the environment. Mitigation response organizations have also been using oil spill modelling as a way to assess and improve their spill response strategies.

Coastal fjords and estuarine regions such as the Strait of Georgia and Juan de Fuca Strait, B.C. are a complex system to model: temperatures and salinities are primarily governed by estuarine processes with, in this case, the Fraser River entering as surface waters and oceanic waters from the Pacific entering as a deep inflow. Nonetheless, these regions are also areas with high hydrocarbon shipping traffic proposed, serving the international market, particularly Asia. The understanding of the behaviour and fate of an oil spill in this area has become a sensitive topic to be addressed by several oil shipping companies.

Although the dominant currents affecting an oil spill are the surface currents, the best way to obtain realistic currents is to use a three-dimensional model. In this way, processes such as wind-driven currents, river plumes and large-scale estuarine circulation are correctly included in the calculation of surface currents. The use of an oil spill tracking model which simulates the weathering of the oil and the advection based on three-dimensional currents is therefore the most appropriate tool to conduct a realistic assessment.

For this paper, the main models used were:

1. A three-dimensional hydrodynamic model, H3D;

2. A wave model, SWAN; and

3. An oil spill tracking model, SPILLCALC.

Tetra Tech's proprietary oil spill model, SPILLCALC, was used for the simulations described here. SPILLCALC can be configured as either a stand-alone model, that relies on other models and observational data bases, or can be embedded within H3D, in which case a full three-dimensional simulation is possible.

A hypothetical spill site, located in central Juan de Fuca Strait, between Vancouver Island, B.C. and Washington State, was selected for the spill simulations described in this paper, as shown in Figure 1. The modelling domain covers the Strait of Georgia, Haro Strait and Juan de Fuca Strait, as well as the adjacent continental shelf and extends to Johnstone Strait in the north. First, the models used for this paper are briefly described. Then, a basic analysis of the local oceanography is conducted. Finally, numerical modelling results are presented. Stochastic modelling is conducted for winter and summer seasons to show the differences in the spill behaviour based on the season. A specific, or deterministic, simulation is then conducted within the same timeframe, where the impact of three-dimensional currents and stratification are shown on the behaviour of dissolved oil.

Surface currents for this paper were hindcast using Tetra Tech's proprietary three-dimensional hydrodynamic model, H3D. This model is derived from GF8 (Stronach et al., 1993) developed for Fisheries and Oceans Canada. H3D has been used on several studies along the B.C. coast, including the Enbridge Northern Gateway Project and the Trans Mountain Pipeline Expansion Project. An extensive application of an operational version of this model to the St. Lawrence Estuary is described in Saucier and Chassée (2000).

H3D has been extensively validated in the area of study (Strait of Georgia, Juan de Fuca Strait, Fraser River). The following key points provide further information on the hydrodynamic characteristics of the model.

• Tidal constituents from the Canadian Hydrographic Services (CHS) and global tidal models (Schrama and Ray, 1994) were used to provide water level data at the oceanic boundary of H3D;

• Wind forcing is derived from coastal Meteorological Service of Canada stations, NOAA stations and moored buoys. More details are provided in the following section “Wind Forcing for H3D and SPILLCALC”;

• River inputs from 50 rivers and creeks throughout the model domain were incorporated);

• Meteorological data are also needed to compute heat flux into the waterbody and thus its temperature structure were derived from MSC station data;

• Turbulence modelling based on the Smagorinsky shear-dependent formulation in the horizontal (Smagorinsky, 1963) and Mellor and Yamada Level 2 in the vertical (Mellor and Yamada, 1982);

• Oceanic boundary conditions for salinity and temperature were available via models maintained by the Alaska Ocean Observing System (AOOS): Temperature and salinity data from the AOOS Pacific Northeast ROMS implementation (12 km resolution) was obtained at 6-hour intervals and used to force the open boundary during inflowing currents; and

• H3D was implemented on a 1 km x 1 km grid, rotated to align with the major axis of the Strait of Georgia.

The oil spill model, SPILLCALC, requires wave conditions as an input to its weathering processes. Wave conditions for the simulation period were hindcast using SWAN version 40.72 (Booij et al., 2006). For consistency with the hydrodynamic inputs, wave conditions were simulated on the same set of computational grids as were used for the hydrodynamic modelling.

SPILLCALC is a time-stepping model that computes the motion and weathering of liquid hydrocarbon spills. SPILLCALC uses currents from the H3D model to advect the spill. The model has been recently expanded to simulate the behaviour and fate of diluted bitumen (dilbit). Because of the wide range of chemical and physical properties of the various constituents of dilbit, pseudo-components were defined based on the Canada Wide Standard, but with a greater number of pseudo-component bins, and with the addition of light fractions such as volatiles (C₁ – C₅) and heavy components such as resins and asphaltenes, for a total of 17 pseudo-components. This pseudo-components definition was supported by a detailed chemical analysis of a diluted bitumen sample. The pseudo-component approach significantly reduces the number of individual constituents that must be considered for evaporation and dissolution, compared to considering each individual chemical species in the mixture.

Oil released on the water surface is represented as a large number of independent floating particles, referred to as slicklets. A total of 50,000 slicklets was used in this study. Individual slicklets are not intended to be physically meaningful. Instead, the cloud of slicklets as a whole is the area covered by the spill, and its progress is the spill's dispersion and trajectory. Each slicklet knows its location (longitude, latitude), its volume and the volume fraction of each pseudo-component, its age, and whether or not the oil is in the form of a tar ball. SPILLCALC includes several weathering processes, as shown in Figure 2. A focus has recently been on quantifying the key weathering processes for diluted bitumen.

In SPILLCALC, there are two processes that determine the rate of evaporation: first, the standard bulk aerodynamic process, which calculates the mass flux from the surface slick based on wind speed, temperature, equilibrium pressure for each the pseudo-component and molar concentration of the pseudo-component in the total product. This method is used in NOAA ADIOS2 model, for instance.

However, SPILLCALC includes an additional mechanism to account for the effect of the slow rate of molecular diffusion within diluted bitumen. Molecular diffusion is responsible for bringing the lighter fractions to the evaporating surface, to replace the losses due to evaporation. In general the rate of molecular diffusion through the vertical extent of the slick is considerably slower than the rate of evaporation from the surface, so that in fact the controlling mechanism is the internal diffusion process. A ten-layer vertical diffusion equation is implemented for each SPILLCALC cell, for each pseudo-component, to calculate the flux of volatiles to the evaporating surface.

Breaking waves drive small droplets of the oil into the water column. Depending on the natural turbulence in the water and the size and density of the droplets, the dispersed oil will generally stay suspended in the water column and will be prevented from resurfacing as long as the dispersing mechanism, breaking surface waves, remain active. When wind and waves die down, the dispersed oil will generally rise to the surface. The process of vertical dispersion was implemented in SPILLCALC using equations developed by Delvigne and Sweeney (1988). The opposing process of resurfacing was implemented in SPILLCALC using the equations developed by Tkalich and Chan (2002). The characteristics of oil droplets for vertical dispersion and resurfacing such as diameter were kept constant over the simulation. A unique feature of SPILLCALC is that the wave field is generated by a reliable and widely used wave model SWAN, whereas many spill models estimate waves from wind speed and fetch. The use of SWAN provides much more realistic wave energy for computing vertical dispersion.

The shoreline is based on British Columbia and Washington State data bases, and includes shore location, coastline type, and a value for oil retention. Oil retention was calculated based on shore types and the known properties of diluted bitumen, especially its relatively high viscosity (Coastal & Ocean Resources, 2013). Refloating of beached oil was not considered in these simulations.

A potential issue of concern is the extent to which oil would come into contact with intertidal sand and mud flats and adversely affect benthic invertebrates and bio-films. In addition to entering beach and mud flat sediment via the shore contact process, SPILLCALC contains an algorithm to simulate stranding of oil as water levels fell below the level of the beach or sand flat.

SPILLCALC keeps track at each time-step of the density of each slicklet based on its composition. Depending on the initial density of the product (927 kg/m³ for diluted bitumen), its composition and the receiving environment (i.e. marine, estuaries or freshwater), the potential for submergence varies. In the three-dimensional simulations, if a slicklet does become submerged, it is then advected horizontally and vertically based on water currents at the depth where the water density and slicklet density are the same. If the slicklet density is greater than any part of the water column, the slicklet is considered sunk, and is transferred to the sea bed.

The formation of oil-mineral aggregates is another process that can affect the behaviour of an oil slick. In river and estuary areas, where the fine sediment load is usually higher than in the ocean, the interaction between oil and fine sediment is crucial in assessing the impact of a spill on the environment. The method used in the SPILLCALC model follows the approach proposed by J.R. Payne (1987) and incorporates the effect of water turbulence. The model was calibrated with experimental data report by Khalifa et al. (2008). Sediment concentration over the model domain was obtained from the H3D model. Monthly concentrations of sediment were used as in input into H3D, based on samples collected at different times and locations in the Fraser River.

Emulsification is a process whereby oil and water co-mingle and form an emulsion, usually requiring wave energy to mix the two liquids. Formulas for the water uptake and the emulsion stability were proposed by Mackay et al. (1980) and Mackay and Zagorsky (1982) respectively. Amongst others, the emulsification has a strong impact on the evaporation process. The inhibition of evaporation increases with increasing water content and slick thickness. With respect to impact on evaporation, SPILLCALC follows the method developed by Ross and Buist (1995): hydrocarbon evaporation is assumed to have a linear relationship with the water content, decreasing as the water content rises.

Some of the lighter hydrocarbon fractions, such as benzene, are soluble in water; they will dissolve into the underlying water column. The potential for dissolution is a function of the pure component solubility, the mole fraction of the hydrocarbon in the source (the spill) and the receptor (surface waters of Juan de Fuca Strait) and the mass transfer coefficient. The rate of dissolution is computed according to the equation published by MacKay and Leinonen (1977) and uses their value for a mass transfer coefficient: 2.36 10⁻⁶ m/s.

A biodegradation module was also incorporated in SPILLCALC. Since the initial bacteria population is rarely well known, SPILLCALC uses a first order bacterial decay process in which the rate of oil biodegraded is proportional to the initial mass of oil and an empirical decay coefficient.

Wind forcing causes both currents and water level differences. Consideration of wind forcing is also important because wind energy has a notable effect on vertical mixing, and therefore scalar distributions. For the modelling described in this paper, wind stresses acting at the water surface are derived from wind records collected from coastal Meteorological Service of Canada stations and moored buoys. Offshore winds from the North American Regional Reanalysis (NARR) were also used due to the scarcity of offshore observations. Winds from NARR matched the buoy data well, but did not accurately represent winds inland of Vancouver Island so were only used offshore. The raw or reanalysis data were processed into hourly time-series of over-water winds at the observation points, and then spatially interpolated inside the three dimensional hydrodynamic model, H3D. A simple inverse distance weighting scheme in H3D is used for the spatial interpolation of hourly wind data from the available locations. The same interpolation method was used by the oil spill model, SPILLCALC.

The prevailing winds at the study site exhibits strong seasonality, typically originating from the east during winter and from the west during summer (as recorded at Race Rocks, B.C.). The winds blowing along the Strait are frequently up to 10 m/s and occur almost continuously in spring and summer but only intermittently in fall and winter. The winds coming off the land, however, are typically less than 5 m/s and dominate the fall and winter periods. Currents in Juan de Fuca Strait during summer are dominated by the large-scale two-layer estuarine flow in these areas, whereas in Haro Strait, currents tend to be more tidal. The seasonality in currents at the study site reflects the seasonality in river inflows into the Georgia-Fuca basin.

Figure 3 shows the temperature and salinity profiles at the study site in Winter (January 30, 2012) and in Summer (August 15, 2012). The stratification can be well observed in the first 50 m depth. As a result, the two periods of interest are winter and summer: stochastic and deterministic spill modelling will be conducted during both of these periods. For the purposes of this study, the winter period was defined as January, February and March and the summer period as July, August and September.

The entrance to Juan de Fuca Strait from the Pacific Ocean shows a strong seasonal variability due to the seasonal changes in the shelf circulation: the Vancouver Island Coastal Current and Davidson Current are dominant in winter and head northward along the west coast of Vancouver Island, whereas the Shelf Break Current heads southward and is dominant in summer.

A stochastic simulation is essentially a collection of separate specific spill simulations. In each separate and independent simulation, most aspects of the spill, such as location and volume, are held constant but the start time and date of the spill are varied so as to capture the full range of environmental conditions likely to be encountered. The resulting set of spills can be analyzed to obtain statistical information on the extent of the spill, the probability to reach a certain area, and the amount of shoreline (minimum, average, median, maximum length) that can be potentially contacted. Stochastic modelling is widely used to develop an understanding of the likely behaviour of an oil slick. Typically, the major driving force for slick motion is wind-driven currents, and it is fairly common to randomly select a number of scenarios, i.e., random sampling of a wind dataset should produce a smaller number of wind events to be modelled, but with the same statistics (means, max, etc. ) as the original series. Because of the strong seasonal signal in the estuarine waters under study, it was thought more important to incorporate sufficient simulations to adequately describe a complete year in a statistical sense, rather than random parts of several years. As well, provision of adequate high quality boundary data, such as that from Alaska Fisheries, was limited to the period 2011–2012. A meteorological analysis was conducted and confirmed that the period 2011–2012 is a representative year, i.e. neither overly stormy or overly calm compared to other recent years.

For the winter period, 364 independent simulations were conducted; for the summer period, 368 independent simulations were conducted at the study site, representing three complete months in each case. Each independent simulation released 16,500 m³ of diluted bitumen over 13 hours at the study site, corresponding to the loss of one tank of a partly loaded Aframax tanker. Each spill was tracked for 15 days.

The behaviour of the slick is very different from one season to another, due to the prevailing winds and surface currents. Figure 4 shows the probability for an area to be contacted by the oil during the winter period. The oil is mainly confined to Juan de Fuca Strait. Less than 20% enters Puget Sound, Haro Strait or exits Juan de Fuca Strait into the Pacific Ocean. There is a 10% probability for the oil to travel along the west coast of Vancouver Island. This transport is associated with the Vancouver Island Coastal Current, heading northward in winter.

A completely different slick behaviour can be observed during the summer season in Figure 5. Because of the large-scale estuarine flow, most of the oil is pushed westward of the study site towards the Pacific Ocean. About 40% exits Juan de Fuca Strait. When entering the Pacific Ocean, the oil moves then southward towards Washington, as a result of prevailing southward shelf flow (summer shelf break current). There is a vanishingly low probability that Haro Strait would be affected during the summer season.

The mass balance after 15 days of tracking time for each of the simulations is similar for the winter and the summer period. Most of the oil ends up on shore, about 67%, since there is no mitigation consideration in this study. About 20% of the oil evaporates; the difference between summer and winter primarily lies in the area covered by the spill and wind speed. The OMA represents the percentage of oil that formed oil-mineral aggregates. No dispersed oil is observed at the end of the simulation. The dispersion process is active during parts of the simulation, based on wind and wave energy; however the resurfacing process follows shortly after the reduction in wave energy. The dispersion process should be considered as an intermediate, but not final, weathering process.

Other outputs from the stochastic model such as travel time are available: first time for a specific segment of shoreline to be contacted, as well as probability maps of oil presence on the water after 6 hrs, 12 hrs, 24 hrs and 48 hrs, provide a good understanding of the travel time and probabilistic evolution of the slick.

The deterministic modelling simulates the behaviour and fate on a three-dimensional computational grid: this approach allows simulation of the fate of dissolved oil in the water column combined with the movements of the slick on the surface. Two simulations were conducted: one in winter, spill released on January 30, 2012, and a second in summer, spill released on August 15, 2012. The release of 16,500 m³ of diluted bitumen oil was assumed to occur over 13 hours. First, the behaviour of the spill on the surface is discussed, followed by a study of the dissolved components of oil.

The movements of the slick on the water surface are followed on a variable time-step based on current drift velocity. The extent of the slick and its thickness are shown after 13 hrs, 24 hrs, 2 days and 7 days as snapshots for both seasons. Figure 6 shows these results for the winter season and Figure 7 for the summer season.

Weak winds in winter combined with the Vancouver Island long-shore current carry the oil out of Juan de Fuca Strait onto the Pacific Ocean along Vancouver Island, as shown in Figure 6. This feature was observed in the stochastic modelling results. About 3% of the oil is left on the water after 7 days.

In summer, as shown in Figure 7, currents are dominated by the combined effects of a large-scale estuarine flow and local winds. The estuarine flow initially brings the oil westward, but, after two days, strong winds blowing from the northwest push the slick towards the Washington coast, where the slick is driven onto the shore and retained there, resulting in less than 1% of the oil being left on water after 7 days. These differing patterns of slick behaviour, and there day-to-day variability, underline the importance of using observed winds and currents from a three-dimensional model which can resolve such complexity between tide, estuarine flow and wind forcing.

The second section of this analysis focuses on the dissolved oil. Since the ecological implications are much greater with toxic components, this paper focuses on dissolved benzene, although more complete analyses have been done for all toxic components such as dissolved TEX (Toluene, Ethylbenzene and Xylene). Benzene is an important constituent to consider because of its relatively high solubility.

A comparison is conducted between dissolved benzene resulting from a winter spill and from a summer spill. For each simulation, three transects were selected and the modeled benzene concentration plotted. All figures have the same scale, with a cut-off at 0.001 mg/L or 1 μg/L. For comparison, the benchmark from the Canadian Council of Ministers of the Environment (CCME) for the protection of marine aquatic life for benzene is 110 μg/L.

Figure 8 shows the benzene concentration in the surface layer at the end of the release, i.e., after 13 hours, during the winter period. At this time, about 1% of the oil dissolved, i.e. 107 m³, (9 m³ of benzene) although a further 533 m³ of oil will dissolve over the duration of the simulation. Three transect line, identified as the west, central and east transects, are shown in Figure 8. Vertical benzene concentrations along the west, central and east transects, are shown in Figure 9, from the point of view of an observer looking towards the east, and for the top 30 m of the water column. The contour lines in the transect figures are temperature. The benzene concentration mixes in the water column between 20 m to 30 m depth depending on transect.

Similarly, Figure 10 and 11 show the surface benzene concentration after 13 hours during the summer period: 126 m³ of oil (9 m³ of benzene) is dissolved at this time. An additional 495 m³ of oil will dissolve over the duration of the simulation. All west, central and east transects show that the concentration is higher near the surface; and the dissolved benzene does not mix in significant quantities below a depth of 5 m, due to the summer stratification. The west transect shows some down-welling occurring, carrying the benzene down to a depth of 20 m along the shoreline.

Complex oceanographic systems such as coastal fjords and estuarine regions required a three dimensional model to adequately simulate current dynamics and oil fate. An extensively validated model for the study area, H3D, was used and provided currents to the oil spill tracking model SPILLCALC. Weathering processes such as evaporation, dissolution, shore retention, OMA formation and biodegradation were simulated. For this paper, an accidental release of 16,500 m³ of diluted bitumen was modeled as a demonstration of model capability.

Stochastic modelling showed a significant seasonal difference in spill trajectories. Oil released at the study site in winter mainly stays confined inside Juan de Fuca Strait with a potential of being carried along the west coast of Vancouver Island. On the other hand, the summer period shows that the oil has a large potential to exit Juan de Fuca Strait and travel along the west coast of Washington State, but is also strongly affected by inflow winds. The differences in the oil trajectory when exiting Juan de Fuca Strait and entering the Pacific Ocean are due to the seasonal changes in the shelf circulation.

Amongst the 17 pseudo-components resolved by the model, dissolved benzene was selected for discussion. Results showed that dissolved benzene concentration is higher near the surface in summer; and does not mix in significant quantities below a depth of 5 m, due to the summer stratification. The mixing of dissolved benzene goes much deeper during winter, down to a depth of 30 m.

These results confirm that an oil spill model coupled with a three-dimensional model is crucial if one wants to look at the behaviour and fate of the oil in complex oceanographic systems.

Two key applications resulting from such a study would be an environmental impact assessment with more accurate analysis, especially regarding toxicity in the water column and on intertidal areas. Spill response would be the second key application to this modelling approach. Indeed, simulating a hypothetical oil release supported by observed winds and modelled three-dimensional currents would greatly help optimizing the mitigation strategy plan. Historical cases pointed out that logistic were a limiting factor in the recovery of the oil. Understanding the oil behaviour and its potential spreading or thickening is crucial for an optimal spill recovery, as it provides time scales for response, and the on-side storage requirements for effective mitigation.

We would like to acknowledge a few persons, who participated to some degrees to this paper. First, we would like to thank Captain Bikramjit Kanjilal (Trans Mountain Expansion Project Marine Lead) and Kinder Morgan Canada. We would also like to acknowledge TERA and Stantec for providing us with the list of pseudo-components. Finally we would like to thank John Harper (Coastal and Ocean Resources) for providing a shoreline database and Elliott Taylor (Polaris Applied Sciences) for conducted lab testing about diluted bitumen.