Evaluation of Oil Spill Modeling in Ice Against In Situ Drifter Data from the Beaufort Sea

The International Association of Oil and Gas Producers (IOGP), in support of the Arctic Oil Spill Response Technology – Joint Industry Programme (JIP), sought to advance and expand oil spill modelling capabilities within the Arctic. For simulation of oil spill trajectory and fate in Arctic waters, oil spill models depend upon coupled ice-ocean (hydrodynamic) and meteorological models to provide wind, current, and ice drift rates, coverage, and thickness. However, available coupled ice-ocean model systems (described below) were designed primarily for evaluating Arctic-scale climate-related issues. They have been of lower spatial/temporal resolution and accuracy than desirable for operational oil spill model forecasts in ice-infested waters. Recognizing this limitation, Phase 1 of the JIP initiative focused on improving ice modeling algorithms and the resolution of ice dynamics. In Phase 2, resulting improvements in oil spill modeling using the updated ice models were examined.

Available coupled ice-ocean models include the Nucleus for European Modelling for the Ocean (NEMO) system, coupled with the thermodynamic-dynamic sea ice model Louvain-la-Neuve Sea Ice Model (LIM; http://www.elic.ucl.ac.be/repomodx/lim/) and the (Towards) an Operational Prediction system for the North Atlantic European coastal Zones (TOPAZ; http://topaz.nersc.no/) developed by NERSC. NEMO provides ocean currents in the Arctic regions at daily time intervals and 1/12 degree resolution (available via http://marine.copernicus.eu/). LIM is used in several global climate models contributing to the assessment reports of the Intergovernmental Panel on Climate Change (IPCC) and in the operational oceanography system MERCATOR (http://www.mercator-ocean.fr). Two versions of the model are available: LIM2 (Fichefet et al., 1997; Bouillon et al., 2009) and LIM3 (Vancoppenolle et al., 2009). The LIM model resolves sea ice concentration, sea ice/snow thicknesses, sea ice drift and sea ice thermal content at the same resolution as NEMO.

TOPAZ4 incorporates the HYCOM hydrodynamic model (Bleck, 2002) coupled with a seaice model using an elastic-viscous-plastic (EVP) rheology (Hunke and Dukowicz, 1997), and a 100-member ensemble Kalman filter (EnKF; Evensen, 1994), a multi-data and multi-variate algorithm for data assimilation. The EnKF assimilates remotely-sensed sea level anomalies, sea surface temperature, sea ice concentration, and Lagrangian sea ice velocities (winter only), as well as temperature and salinity profiles from Argo floats. Wind stress for the TOPAZ4 model is from the ERA-40 (European Center for Medium-range Weather Forecast, ECMWF REANALYSIS) wind model (http://www.ecmwf.inf). The model predictions have been compared to data in broad scale; Sakov et al. (2012) found that TOPAZ4 produced realistic estimates of the mesoscale circulation in the North Atlantic and sea ice variability within the Arctic. The horizontal resolution of TOPAZ4 is approximately 12–16 km in the Arctic Ocean.

Recently NERSC updated the ice modeling approach and specifically the rheology for pack ice, which is coupled with the TOPAZ ocean model. The resulting Lagrangian-based Elasto-Brittle (EB) model (named neXtSIM, Rampal et al. 2009a), which uses a more accurate rheology than classic models (which use EVP), showed significant improvement in performance over the present (Elastic-Viscous-Plastic, EVP) modeling approach used in their operational TOPAZ4 coupled ice-ocean model as well as other existing ice-ocean models. The EB-based neXtSIM model better preserves dynamic features like cracks, ridges, and leads in the ice than traditional Eulerian models (Rampal et al., 2016a). NERSC focused their validation (Rampal et al., 2016b) on examination of transport (as compared to drifters and interpretations of satellite imagery), statistical analysis of ice coverage, and on the degree of dispersion (comparing model predictions to drifters). In the Beaufort Sea, in particular, NERSC noted that the predicted ice movements (generally westward) were too fast in TOPAZ coupled with the standard EVP rheology, but more similar to observations using the neXtSIM model.

This study made use of three refined ice modeling products provided by NERSC, as well as publically-available TOPAZ4 data products, to evaluate potential improvements in modeling oil transport within the Beaufort Sea:

1. TOPAZ-EVP: TOPAZ hydrodynamics and ice model using updated EVP rheology, free-run for 10-years (2001–2010), provided as 6-hourly averages.

2. TOPAZ-EVP-WIM: TOPAZ (free run) hydrodynamics and ice model using updated EVP and MIZ rheology and WIM-generated wave data for hindcasts of 5 selected months (March, June, September, and December of 2008 and May 2009) in a portion of the Beaufort Sea, provided as hourly data.

3. neXtSIM: Updated TOPAZ reanalysis hydrodynamics and ice model with EB rheology in pack ice, run for 10 winter periods (i.e., November 1 – May 15, 2000–2010), provided as 6-hourly averages.

4. TOPAZ4: publically available product from the operational model, in reanalysis mode, with ice model using the standard EVP rheology.

Using trajectory simulations, we have compared transport using these ice-ocean model datasets against available in situ drifter data in the Beaufort Sea from the IABP.

The model algorithms in OILMAP and SIMAP (French McCay, 2003, 2004) have been developed over the past three decades to simulate oil transport and fate under a variety of environmental conditions. The SIMAP model quantifies trajectory, areas swept by floating oil of varying thicknesses, and fates and concentrations of subsurface oil components (dissolved and particulate). Processes simulated include spreading (gravitational and by shear), evaporation of volatiles from surface oil, transport on the surface (in ice and open water) and in the water column, randomized dispersion from small-scale motions (turbulent mixing), emulsification, entrainment and resurfacing of oil as droplets of varying sizes into the water (natural and facilitated by dispersant) or returned to the surface, dissolution of soluble components, volatilization of dissolved hydrocarbons from the surface water, adherence of oil droplets to suspended sediments, adsorption of soluble and semi-soluble aromatics to suspended sediments, sedimentation, stranding on shorelines and landfast ice edges, and degradation. Lower-molecular-weight aromatic hydrocarbons, i.e., mono-aromatics (MAHs) and polynuclear aromatic hydrocarbons (PAHs), are the soluble and semi-soluble components that are most bioavailable to aquatic biota, inducing most of the effects (French McCay, 2002). These and other “pseudo-components” representing volatile aliphatic hydrocarbons are tracked separately from residual oil in the model. Whole oil (containing non-volatile residual oil and volatile components not yet volatilized or dissolved from the oil) is simulated as floating slicks, emulsions and/or tarballs, oil in ice, or as dispersed oil droplets of varying diameter (some of which may resurface). Sublots of the spilled oil are represented by Lagrangian elements (“spillets”), each characterized by mass of hydrocarbon components and water, location, thickness, diameter, density, and viscosity. A separate set of Lagrangian elements is used to track movement of the dissolved hydrocarbons.

OILMAP is a simplified version of SIMAP, designed for operational use and for contingency planning. It uses a reduced number of pseudo-components to represent the oil, and so requires less data to define the oil composition. It tracks the whole oil and volatiles through evaporation, but does not track dissolution or the transport of the dissolved hydrocarbons.

Oil interactions with mobile sea ice or immobile landfast ice involve several processes that affect oil transport and fate. Oil released at or above the water surface, may spill into water and/or onto the surface of the ice. Oil deposited on ice may absorb into surface snow, run off and become trapped between cracks or in open water between ice floes, and/or become encapsulated in the ice. Oil released into and under water may become trapped under the ice in ridges and keels, or build up along and become trapped in sea or landfast ice edges (Drozdowski et al., 2011). Many of these interactions and processes are at a finer scale than can be captured in oil spill models using inputs that are currently available from large scale meteorological, hydrodynamic and coupled ice-ocean models. However, the influence of ice on net transport and fate processes is simulated by considering the location of the oil (on, under or between ice) and potential reduction in oil surface area (i.e., reduced spreading), which changes the wave environment, entrainment rate, oil movement, dissolution, volatilization and mixing in the water column.

When oil interacts with mobile sea ice, some fraction will become contained (either on top, in, or underneath the ice) and will travel with the ice floe (Drozdowski et al., 2011). This oil can drift rapidly and over great distances in the Arctic (Peterson et al., 2008). The fraction of oil moving with the ice versus that in open water depends on conditions and specifics of the release. In some cases, all of the oil becomes frozen in the ice and remains there until the ice melts. This scenario is readily modeled (i.e., 100% of oil drifts with ice). However, in most cases since sea ice can be patchy, only limited amounts of oil may become either encapsulated or trapped (e.g., between ice fragments or under ice in cavities; Drozdowski et al., 2011), depending on ice coverage, subsurface ice roughness, winds and currents, and ice formation/melting dynamics.

Ice coverage information available in coupled hydrodynamics and ice models is typically resolved at relatively large scales (>1 km). While detailed information regarding ice coverage and conditions are not available from these models, the information provided can be used as an indicator of whether oil would move predominantly with the surface water currents or with the ice. A rule of thumb followed by past modeling studies is that oil will generally drift with ice when ice coverage is greater than 30% (Drozdowski et al., 2011; Venkatesh et al., 1990).

The OILMAP and SIMAP models use the ice coverage data (at the available resolution) to determine whether floating (or ice-trapped) oil is transported by the surface water currents or the ice. If the ice coverage is less than a minimum threshold defining the MIZ (MIZ_min, default assumption is 30%), the oil is assumed to move with surface water currents. If ice coverage exceeds MIZ_min, the ice coverage limits the spatial coverage of the oil, and oil is transported with the ice using the ice velocities from the ice-ocean model. Horizontal diffusion coefficients for floating oil are reduced in the ice, proportionately to fraction of ice cover.

In areas and at times where ice cover < MIZ_min, floating oil is transported with surface water currents and a wind drift algorithm to account for wind-induced drift current not resolved by the hydrodynamic model plus Stokes drift caused by wave motions. Wind drift is predicted in the oil spill model based on the modeling analysis of Stokes drift and Ekman flow by Youssef (1993) and Youssef and Spaulding (1993, 1994). According to this algorithm, in the northern hemisphere at moderate wind speeds, floating oil drifts 20° to the right of downwind at about 3.5% of wind speed. In areas where ice exceeds MIZ_min, the ice drift model explicitly accounts for wind drift, and so no additional wind drift is added to the oil transport.

In the model, when floating oil encounters landfast ice it is assumed to trap at or move along the ice edge (depending on the current and wind directions at the location and time). If oil becomes entrapped within landfast ice (by surfacing there or as landfast ice extends over the area), it remains immobile until the ice retreats. When landfast ice is no longer present at the location of trapped oil, the oil is released back into the water as floating oil.

Drifter buoys have been deployed and tracked in the Arctic Ocean since 1978 as part of the IABP to provide met-ocean data in real time. The number of buoys and the coverage varies each year, but an average of 25 buoys are in service at any time (Colony and Munoz, 1986; Rigor and Ortmeyer, 2001). The data are processed by the University of Washington’s Polar Science Center, and are interpolated to produce gridded fields of surface temperature and velocity from drifter positions at 12-hour intervals (0000 UTM and 1200 UTM) daily. The interpolated position, and interpolated ice velocity fields, which are estimated and computed from the buoy positions, are available online from 1979 to present (Rigor, 2002).

Drifter trajectories for buoys in the Beaufort Sea during years and areas where NERSC ice modeling products were available were extracted from IABP records. Oil spill model trajectories, using the NERSC data as input, were compared with the observed drifter trajectories. Simulations for locations in pack ice during winter-spring months were used to test the advancements of the neXtSIM model’s EB ice rheology over the recently updated EVP rheology in TOPAZ-EVP and the publically-available TOPAZ4 model. TOPAZ-EVP-WIM includes the updated EVP rheology and a MIZ rheology that incorporates results of a wave model (Bergh and Ólason, 2015), the performance of which was compared to the other models. Since TOPAZEVP-WIM model results were available for March, June, September and December of 2008 and May of 2009, the comparisons were performed for all IABP observed trajectories available during those five months in the Beaufort Sea area covered by the NERSC products. As the neXtSIM model products were only provided for March, December and May, i.e., months where ice coverage was >90%, we focused the analysis on these “winter “ months.

The differences in the model-predicted versus observed tracks over time were measured to evaluate the model predictive skill. First, the trajectories of ice model and drifter data were plotted as progressive vector diagrams from the location of each drifter buoy at the selected initialization time, i.e., the first time step of the IABP trajectory segment examined. The predicted and observed ice drifter trajectory were compared using a simple ratio of the model predicted path length divided by the observed path length. A dimensionless skill score was also calculated by using the Lagrangian separation distance between endpoints of simulated and observed drifters divided by the observed path length (Liu and Weisburg, 2011; Ivichev et al., 2012). The separation index S is defined as:

where sepDistance_i is the separation distance and ObservedPath_i is the cumulative length of the observed path at time step i. N is the number of (e.g., 12-hour) time steps since the beginning of the simulation. Thus, the separation index S is a measure of the cumulative error in predicted distance traveled, normalized by the cumulative length of the observed path. This approach accounts for the differences in the speeds of movement over the path. Differences in speed were evaluated using the ratio of the modeled to the observed path length (LenRatio).

The skill score SS is defined as:

T is a user-selected tolerance threshold that is typically set to 1. At T = 1, the tolerance for error is equal to the normalized cumulative separation distance, i.e., error should not exceed the magnitude of the cumulative movement. The skill score SS increases as model skill increases, but is zero if the model has “no skill”, where S is greater than the tolerance threshold T.

Multiple drifter experiments were extracted from the IABP data set, trajectories simulated, and skill scores determined. Drifter experiments were taken from continuous observed drifter trajectories as all sequential, non-overlapping, observed paths of individual buoys over 5-day, 10-day and 15-day intervals that were within the ice model data domain and the temporal periods noted above. Five days is an appropriate and desirable forecast period for ice-ocean and oil spill modeling; the 10-day and 15-day intervals were examined to determine model skill farther into the future. The use of non-overlapping time periods ensures that the results are unbiased by auto correlations in the trajectory movements.

The trajectory example shown in Figure 1 is drawn from modeling analyses of 61 drifters that moved through the TOPAZ-EVP-WIM study domain in the Beaufort Sea where we have data from all ice model data products. The predicted trajectory patterns are similar, reflecting the underlying current and wind field inherent in all models, and performance is reasonably good for both TOPAZ-EVP and neXtSIM based simulations. For these and all drifters studied, the transport paths are characterized by periods of relatively constant directional movement, interrupted by rapid changes in direction. There are typically three to five such events during the one month simulation period, although occasionally there are periods where motion is slow and buoys move in seemingly random directions. One can see divergences accumulate over time, but most of the displacements occur at the events. In the example in Figure 1, the neXtSIM EB model prediction is very close to the observed drifter path for the first 10 days as it changed direction 5 times, whereas the predictions using the TOPAZ models are close to the observed path for the first 7 days before diverging (primarily by moving too fast). This example shows the limits of the ice-ocean models’ abilities to forecast as being about 7–10 days, although in other examples accuracy degrades after 3–5 days of forecasting. While the accuracy of individual oil model trajectories projected weeks to months into the future would be expected to be low, in the event of a spill, forecasts could be updated frequently (on a time scale of hours to days) with satellite information, aircraft observations, drifter data, and other observations to improve reliability. The right panel of Figure 1 demonstrates the model trajectories when reinitialized at the observed drifter position on 15 March 2008. In the second 15-day period, all the models correctly predict the turn to the east, but all (the neXtSIM EB and TOPAZ-EVP-WIM less so than TOPAZ4 and TOPAZ-EVP) move faster than the observed drifter.

The ratios of model:observed path length (LenRatio), separation index (S) and skill scores (SS) calculated for each of the 5-day, 10-day and 15-day simulations of IABP drifter segments using the different ice-ocean models are summarized for the winter months in Tables 1–3, respectively. The LenRatio and S are lowest using neXtSIM, highest using TOPAZ4 and intermediate using TOPAZ-EVP and TOPAZ-EVP-WIM as trajectory model inputs. The skill scores SS show the opposite pattern. The average skill scores are generally below 0.7. In segments where the model skill is such that S<1, SS is calculated as >0, i.e., the model has some skill. The differences between skill scores for neXtSIM, TOPAZ-EVP and TOPAZ4 are significant at p<0.05 for the 5-day and 10-day trajectories, and are significant at p<0.10 for the 15-day trajectories, based on a two-tailed t-test for unpaired observations. The differences between LenRatio and S between the models were of less significance (because of higher variance than SS, which treats all poor fitting model trajectories equally as no skill). TOPAZEVP-WIM simulations were intermediate in performance between TOPAZ-EVP and neXtSIM, but the differences were not significant. Thus, the SS results indicate that the neXtSIM model has better skill than both TOPAZ models with standard EVP rheology, but results for neXtSIM do not show significantly better performance than TOPAZ-EVP-WIM simulations with EVP and MIZ rheology.

The LenRatio results for the winter months (Tables 1–3) indicate that the model-predicted ice drift speeds are faster than observed, more so for the TOPAZ models than for neXtSIM. However, examination of results by month showed that the speeds were too high (speeds averaged 5–8 times the observed for the 5-day trajectories, Table 4) for all models in March 2008, but neXtSIM, TOPAZ-EVP and TOPAZ-EVP-WIM were not significantly different from the observed in December 2008 and May 2009, whereas TOPAZ4 speeds averaged 1.5 times the observed (Table 5). Model skill (SS) was much lower for March 2008 than for December 2008 and May 2009 (Tables 4–5). In December 2008 and May 2009, model performance was relatively good for simulations 10 days and 15 days in length, and while neXtSIM performed significantly better than TOPAZ4, TOPAZ-EVP-WIM performed as well as neXtSIM by all three metrics. For the summer months (June and September 2008), the LenRatio, S and SS scores calculated for each of the 5-day, 10-day and 15-day drifter segments are not significantly different using TOPAZ4, TOPAZ-EVP and TOPAZ-EVP-WIM as trajectory model inputs.

The Separation distance (km) between the model trajectory and the observed path provides a simple metric of error for a 5-day, 10-day or 15-day forecast. Separation distances after 5 days are on average 20–22 km using the TOPAZ models (based on results for all 5 months) and 14 days using neXtSIM (in the winter months). The expected separation distances increase to about 45 km days for TOPAZ and 26 km for neXtSIM after 15 days. The coefficient of variation is 70–80% of the mean for TOPAZ and 62–64% of the mean for neXtSIM.

Figures 2 to 5 demonstrate for a single example drifter (#66276), tracked in March and April 2008, the TOPAZ4, TOPAZ-EVP, TOPAZ-EVP-WIM and neXtSIM model-predicted trajectories compared to the observed trajectory, respectively. The model-predicted trajectories are reinitialized to the observed locations every 5 days. The simulations with re-initialization show improved matches to the observations. All of the models perform well in this operational mode where the modeled locations are updated regularly.

Rampal et al. (2016b) compared modeled trajectories using their updated ice model products with IABP drifter paths, examining speeds and diffusivities, finding the EB-based neXtSIM model more closely agreed with the IABP observations than TOPAZ-EVP. They found that in general, trajectories were longer when modeled using the TOPAZ data than observed in the IABP data. Our findings based on the modeled trajectories herein, as measured by the ratio of modeled to observed path lengths (indicative of relative speed) and skill scores, are consistent with their findings related to relative performance. A comprehensive test series involving statistical comparisons of the modeled trajectories to ~400 5, 10, 15 and 30-day intervals of IABP drifter trajectories demonstrated that the new neXtSIM model predictions agree more closely with the observations than the operational status quo models represented by TOPAZ4 using the older (standard) EVP rheology. There was a significant improvement in the forecast performance (as measured by cumulative divergence from the observed buoy path) in TOPAZ accomplished with the TOPAZ-EVP model, and the neXtSIM EB model trajectories were significantly better than both the TOPAZ-EVP and TOPAZ4 models in speed and direction of oil-in-ice movements. On average, the neXtSIM EB model trajectories were statistically similar in length to the observed buoys (i.e., indicating similar speed of ice movements), whereas oil trajectories using the TOPAZ models using standard EVP rheology had significantly longer lengths 1.5 up to 5 times longer than the paths of the observational buoys.

Based on the present analysis, use of the EB-based neXtSIM model would improve the accuracy of oil spill trajectory modeling in areas of fast ice when compared to the use of results from TOPAZ with EVP rheology. From the sample of drifters examined here, it appears that the updated TOPAZ-EVP and TOPAZ-EVP-WIM models perform better than the publically-available TOPAZ4 reanalysis version, but this may be due in part to differences in the wind forcing used (TOPAZ4 used the interim ERA, while TOPAZ-EVP used the final reanalysis ERA) and/or to the daily averaging used to deliver (on the web) the TOPAZ4 product as compared to 6-hourly or hourly data delivered in the TOPAZ-EVP and TOPAZ-EVP-WIM products, respectively. Provision of TOPAZ-EVP model data at smaller time steps (than daily, as on the web server, or even than 6-hourly) and without time averaging would provide improved performance in oil spill modeling.

Further analysis would provide more insight as to the locations, times and conditions under which these noted improvements in model performance should be expected. In using these ice model products for oil spill forecast modeling, one needs to consider the increasing divergence between model predicted and observed trajectories of drifters in ice that accumulate over time using any of the ice model products. Separation distances after 5 days are on average 20–22 km using the TOPAZ models and 14 km using neXtSIM (in the winter months). The accuracy of individual oil model trajectories projected more than a few days into the future would be expected to be low. In practice, operational oil spill trajectory forecasts should be frequently reinitialized to reduce the forecast errors that are accumulated from the initial locations (Liu and Weisberg, 2011; Ivichev et al., 2012). However, the overall transport patterns and results of an ensemble of trajectories would provide useful information for planning and risk assessments based on typical current and ice movement patterns.

This research is part of the International Association of Oil and Gas Producers (IOGP), Arctic Oil Spill Response Technology – Joint Industry Programme (JIP) supported by nine oil and gas companies – BP, Chevron, ConocoPhillips, Eni, ExxonMobil, North Caspian Operating Company, Shell, Statoil, and Total.