Horizontal transport in oil-spill modeling
HHorizontal transport in oil-spill modeling
Rodrigo Duran a,b, ∗ , Tor Nordam c,d , Mattia Serra e,f , Chris Barker g a National Energy Technology Laboratory, U.S. Department of Energy, USA b Theiss Research, USA c SINTEF Ocean, Trondheim, Norway d Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway e School of Engineering and Applied Sciences, Harvard University, USA f Department of Physics, University of California San Diego, USA g O ffi ce of Response and Restoration, Emergency Response Division, National Oceanic and Atmospheric Administration, USA Abstract
Simulating oil transport in the ocean can be done successfully provided that accurate ocean currents are availablethis is often toobig of a challenge. Deficient currents can sometimes be remediated by parameterizing missing physicsthis is often not enough.In this chapter, we focus on some of the main problems oil-spill modelers face, which is determining accurate trajectories whenthe velocity may be missing important physics, or when the velocity has localized errors that result in large trajectory errors.A foundation of physical mechanisms driving motion in the ocean may help identify currents lacking certain types of physics,and the remedy. Recent progress in our understanding of motion in the upper centimeters of the ocean supports unconventionalparameterizations; we present as an example the 2003 Point Wells oil spill which had remained unexplained until recently. Whenthe velocity realistically represents trajectory forcing mechanisms, advanced Lagrangian techniques that build on the theory ofLagrangian Coherent Structures can bypass localized velocity errors by identifying regions of attraction likely to dictate fluiddeformation. The usefulness of Objective Eulerian Coherent Structures is demonstrated to the oil-spill modeling community byrevisiting the 2010 Deepwater Horizon accident in the Gulf of Mexico, and predicting a prominent transport pattern from animperfect altimetry velocity eight days in advance.
Keywords:
Oil spill, Lagrangian Coherent Structures, chaotic trajectories, forecast, Objective Eulerian Coherent Structures,TRAPs, altimetry
Contents1 Introduction 22 The physics, the mathematics and the numerics 23 Overview of oil transport in the ocean 4 ∗ Corresponding author.
Email address: [email protected] (Rodrigo Duran)
Preprint submitted to Elsevier September 29, 2020 a r X i v : . [ phy s i c s . a o - ph ] S e p . Introduction Successfully forecasting the movement of oil during an oilspill, along with knowing oil’s current location, are the key in-gredients to solving one of the most important problems thatneeds to be solved during emergency response operations: whereis the oil heading? Oil’s current and future locations are criticalfor emergency planning and response, including recovery andcontainment. Reconstructing past spills through hindcasts isalso critical to understand an oil spill, for modeling improve-ments, for environmental impact studies and forensics. Thechallenge of successfully forecasting or hindcasting an oil spillis a considerable one. There is no guarantee that a simulationwill be successful, and deviations between simulated transportand observations are the norm rather than the exception. Thischapter describes some of the main reasons why accurately sim-ulating trajectories in the ocean is complicated and presents re-cent progress along di ff erent fronts that help remediate some ofthe problems.One of the problems is related to local deficiencies in the ve-locity that propagate during integration when computing trajec-tories; errors induced by inaccuracies often grow exponentiallydue to the unstable nature of ocean circulation. For this prob-lem we present modern techniques called Objective EulerianCoherent Structures (OECS; Serra & Haller, 2016) that can by-pass errors in the velocity, identifying instantaneous attractingregions that influence transport exceptionally, and that are com-puted without the need to integrate the velocity. This method al-lows trajectory forecasting without future information, but canalso be applied to ocean-model forecasts to limit the e ff ect ofdiscrepancies between the forecasted flow and reality.The other type of problem we examine is related to a veloc-ity assumed to simulate oil’s motion, but that is missing someof physical processes that drive motion, with a focus in theupper centimeters of the ocean where atmospheric influenceis strongest. The solution is to parametrize physics that drivemotion when necessary. Studies in the last few years have im-proved our understanding of the velocity within a fine surfacelayer, providing information that is somewhat at odds to com-mon practice in oil-spill modeling. As an example, we showhow progress in our understanding can explain the transport ofoil during the 2003 Point Wells spill in the Salish Sea whichhad remained unexplained for 15 years. Also discussed are re-cent studies showing that parametrized near-surface processesare often responsible for oil beaching, and finally, the potentialfor water subduction at the sea surface to indicate regions ofinterest for oil recovery.This chapter focuses on hindcasts and forecasts, a problemcharacterized by the need to replicate or anticipate observed oiltransport. A brief review of the basics of oil-spill modelingis given in section 2. Section 3 gives an overview of trans-port in the ocean, the unstable nature of ocean currents, and adescription of velocity products typically available to simulateoil transport. Section 4 is about motion in the upper layer ofthe ocean, describing recent progress due to numerical simula-tions and unique at-sea experiments. Section 5 describes someof the novel tools that may help overcome velocity errors that would result in erroneous trajectories, and exemplify their useby revisiting transport patterns observed during the DeepwaterHorizon accident in May 2010. We conclude in section 6 withconclusions and an outlook of what progress can be expected,and how the techniques presented here fit into that picture.
2. The physics, the mathematics and the numerics
Oil-spill modeling is often a multidisciplinary endeavor; itis not uncommon for biologists, chemists, physicists, geogra-phers, mathematicians, oceanographers, engineers, and com-puter scientists to work together. It is therefore helpful to beginclarifying, somewhat informally, the basic physics, mathemat-ics, and numerical solutions used to simulate the transport ina fluid that ultimately results in an oil-spill simulation. Thephysical mechanisms that drive motion within the ocean are de-scribed in later sections.Fortunately, the physics used to simulate transport forcedby a vast variety of oceanographic processes are well known.Consequently, the mathematical equations are also well known(e.g. studied in most introductory partial di ff erential equationcourses). The mathematics of oil transport boils down to solv-ing the advection-di ff usion equation, also known as the dispersion-di ff usion equation. In fluid dynamics and physical oceanogra-phy, the e ff ects of advection and di ff usion are often referred toas stirring and mixing, respectively. In oil-spill modeling, thenumerical solution of these equations is done in a Lagrangianframework, i.e. by computing trajectories. However, it is il-lustrative to introduce the equations of mathematical physics inEulerian terms, i.e. as a function of space and time, and returnto the Lagrangian formulation when discussing the numericalsolution.Stirring within a fluid causes a tracer to deform into streaksand swirls while the along-path concentration of the tracer doesnot change. E ff ectively, the tracer is redistributed with the ve-locity field while preserving its concentration. Such behavioris called conservative, assuming the velocity divergence is neg-ligible it satisfies the following evolution equation for a tracer C : ∂ C ∂ t + u · ∇ C = u is the two-dimensional velocity of the fluid. This equa-tion is known as the advection equation or dispersion equation.Some scientists may use the word convection instead of advec-tion, although oceanographers often reserve the term convec-tion for a di ff erent type of physics (vertical motion related tobuoyancy changes). The advection equation describes how aconcentration evolves as the velocity field acts upon the gradi-ent of the concentration, causing a redistribution of the tracer.Mixing refers to the di ff usion of a tracer with concentration C , typically simulated with the di ff usion equation: ∂ C ∂ t = ∇ · ( κ ∇ C ) (2)If κ is a molecular di ff usion coe ffi cient, then (2) representsthe mixing of a tracer due to molecular collisions. By itself,2his is an ine ffi cient method of mixing a fluid. In the ocean,however, an eddy di ff usion coe ffi cient is used for κ , several or-ders of magnitude larger than the molecular coe ffi cient that ischaracteristic of the fluid. From a physics point of view, thelarge coe ffi cient means that the di ff usion equation is model-ing the mixing of concentration due to the collision of waterparcels, not molecules. This is an ad hoc way of parameteriz-ing small-scale stirring and overturning of water parcels as theyundergo turbulent motions. Ironically, the physics simulated byequation (2) are well understood, yet the ocean physics that itparametrizes includes a variety of processes that are di ffi cult toeven measure (e.g. Moum & Rippeth, 2009). It is a fortunateturn of events that results from using a well-understood equa-tion such as (2) are satisfactory for many purposes, includingoil-spill modeling.To simulate transport in the ocean, both stirring and dif-fusion are often used simultaneously: the evolution equationfor the transport of a tracer is then the advection-di ff usion (ordispersion-di ff usion) equation: ∂ C ∂ t + u · ∇ C = ∇ · ( κ ∇ C ) (3)In the case of oil-spill modeling, an additional source termcan be added to represent the addition of oil as it is spilling intothe ocean; a sink term can also be included to represent the re-moval of oil. This chapter is concerned with the transport ofoil, and we will not consider sources or sinks. Salmon (1998)presents a discussion (his section 14) on stirring and mixing,describing their individual and simultaneous e ff ects on tracervariability. Stirring, di ff usion and their interplay is also dis-cussed in section 7.3 of Tennekes & Lumley (1972).The eddy di ff usion coe ffi cient, also known as the turbu-lent di ff usion coe ffi cient, is an unknown that needs to be deter-mined. The production of turbulence, and therefore the magni-tude of the eddy di ff usion coe ffi cient, depends on many factorsincluding the spatial structure of seawater’s density, the spa-tial structure of the velocity field, heating or cooling of waterparcels, and Earth’s rotation. Turbulence in the ocean is oftenproduced by instabilities that range in length from centimetersto hundreds of kilometersthe interested reader is referred to thefree book on ocean instabilities by Smyth & Carpenter (2019).Due to a large number of instability types, the large range ofspatial and temporal scales at which they occur and interact,and their often anisotropic nature, determining an adequate tur-bulent di ff usion coe ffi cient is a di ffi cult problem.Numerical ocean models require accurate mixing of mo-mentum, heat, and salinity to produce a good simulation; theyuse sophisticated turbulence closure models that are computa-tionally intensive and that require considerable more informa-tion than what is typically available during oil-spill simulations.Fortunately, the need for mixing parameterizations in oil-spillmodeling is fundamentally di ff erent than in numerical oceanmodeling and is not nearly as consequential. In oil-spill mod-eling, the di ff usion equation is used to simulate the small-scalespreading of oil caused by oceanographic processes that are notresolved by the velocity in equation (1). In oil-spill modeling, the most e ffi cient way to determinethe eddy di ff usion coe ffi cient is to choose a constant coe ffi -cient that matches the observed spread of oil. This is a strat-egy used for hindcasts and forecasts where the main objectivefor the simulation of di ff usion is matching the spread of oil asobserved through overflights, ships and satellites. However,during response forecasts, the di ff usion coe ffi cient is chosento err on the high side, so that the simulated impacts are un-likely to underestimate observed impacts. For example, stan-dard practice for NOAA’s O ffi ce of Response and Restorationis to use a random walk with a constant di ff usion coe ffi cientas described, depending on the type of simulation. Despite itsrather crude and ad hoc nature, it is often important to simulatedi ff usive processes in oil-spill modeling because 1) it is neededto match observed oil spreading about the trajectories result-ing from equation (1), 2) it provides a least-regret conservativeestimate for the spreading and impact of oil during responseforecasts, and 3) for simulations without wind, it provides oneway for oil to beach which is otherwise not available from mostocean-current velocity products (e.g. ocean models set to zerothe velocity normal to the coast near the coastline, although itmay cause beaching if there is a mismatch between the velocityproduct and the model coastline, or due to numerical instabili-ties). Beaching due to di ff usion is most appropriate in the surfzone, the e ff ects of which are not typically simulated in oceancirculation models. Other (more realistic) mechanisms that maydrive oil beaching are discussed in section 4.There are types of oil-spill simulations that may need an au-tomated method of determining an eddy di ff usion coe ffi cient.For example, ensemble-type modeling uses a large number ofoil-spill simulations to evaluate the environmental impact of anoil spill in a statistical sense; such simulations may choose toinclude di ff usive processes. The solution for oil-spill simula-tions that compute an eddy di ff usivity as part of the problem (in-stead of chosen to match observed spread) is described in Ap-pendix A. Further discussion on stochastic parametrizations ofdi ff usion and their implementation can be found in the techni-cal documentation for NOAA’s GNOME model (Zelenke et al.,2012) and Duran (2016).The easiest way to numerically solve equation (3) for oiltransport simulations, is to separately solve (1) and (2) andthen add the motion induced by each process to obtain the oil’smovement. Solutions to (1) and (2) are often found separatelyin Lagrangian terms, i.e. by computing oil-parcel trajectories,rather than in Eulerian terms where the equations are solvedover a fixed numerical grid. Because the advection part oftransport preserves concentrations, the simulation of advectionreduces to integrating the velocity field to obtain trajectoriesgiven an initial time and position. Thus, the typical approach tosimulate the advective part of oil’s trajectory is the solution x ( t )to the equation: d x d t = u ( x ( t ) , t ) x (0) = x . (4)Equation (4) can be solved with regular numerical methods forordinary di ff erential equations, although there are some specific3onsiderations due to the discrete nature of the velocity databeing integrated (e.g. Nordam et al., 2018; Nordam & Duran,2020).The di ff usion part is typically modeled as a random walk,numerically simulating random motions induced by the colli-sions being modeled with equation (2). The resulting motion isthen added at each time step to the solution of (4). Typically,the precision needed for u in (4) is more consequential than thespread of oil modeled with di ff usion. Because of this, in thischapter, we will focus on the advective part of oil’s transport.In this approach, once the advection part is satisfactory, furthersimulations will add di ff usion to the advective part.The Eulerian representation of fluid flow and its associatednumerical representation could, in principle, be used for oil-spill modeling. In practice, however, it is much more e ffi cientto simulate oil transport using the Lagrangian representation.Among the problems inherent to the Eulerian representation isthe need to set up a numerical grid for each domain. Also,the Eulerian representation is more computationally intensivebecause it requires solving the advection-di ff usion equation ateach point on the grid, regardless of whether there is oil there ornot, while in the Lagrangian approach trajectories are integratedonly for existing oil parcels. In the Eulerian approach, the non-linear term in (3) can be problematic for numerical methods(e.g. Durran, 2010, sections 3.3, 3.4, 3.5.1, 3.5.2, 5.10), whilethe Lagrangian approach only requires integrating an ordinarydi ff erential equation which is, for the most part, straightforwardand highly accurate. Finally, the Lagrangian approach workswell with a velocity field saved at a series of discrete times, thisis convenient because it readily allows additional experimentsusing the same velocity field. A comparison of advection re-sults from Eulerian and Lagrangian formulations can be foundin Wagner et al. (2019).Oil-spill models can account for other processes related tothe fate of oil separately, for example, oil droplets breaking intosmaller-size droplets, oil dissolution at depth, or oil evaporationat the surface. These processes can a ff ect the oil’s buoyancy.As oil-spill and blowout models increase in complexity, the ef-fects of vertical motion will be included. Oil’s vertical motionmay be induced by ambient conditions such as waves, or byoil’s buoyancy (some oils are denser than water, some are lessdense), droplet size, and weathering. Further details on some ofthe processes causing vertical motion can be found in Nordamet al. (2020).In this chapter, we assume that oil’s vertical location is known,whether varying or fixed. This is a valid approach because:1) For some spills, modeling horizontal trajectories alone maybe enough to get satisfactory results. 2) When vertical motionis important, vertical and horizontal components of a parcel’strajectory are computed separately. This is because the mecha-nisms forcing horizontal motion, tend to be di ff erent from mech-anisms forcing vertical motion, and the resulting trajectory com-ponents can be added to give an updated location for the next in-tegration step. Whenever the vertical dimension is included, theproblem of determining oil’s horizontal motion is more compli- cated. The vertical position of oil must first be determined to beable to sample the correct horizontal velocity, additionally, thevelocity u driving horizontal motion will now need to be accu-rate at di ff erent vertical levels.For the purpose of this chapter’s discussion, horizontal trans-port will be defined as the motion of oil at a constant depth,whether at the surface or deeper down within the water column.Horizontal motion in this chapter also means motion along sur-faces of constant density (lateral motion in oceanographic ter-minology) as long as it does not involve abrupt vertical dis-placements, e.g. where constant density surfaces rise to inter-sect the ocean surface. Some comments will be made regardinghorizontal transport of oil at density fronts in subsection 4.3.Ocean currents are driven by a variety of physical processesthat, for oil-spill modeling, can be conveniently classified bytheir vertical location within the water column. The range ofphysical processes driving motion can then be narrowed downaccording to the vertical location of oil.
3. Overview of oil transport in the ocean
In practice, transport of oil in the ocean is successfully sim-ulated as the movement of parcels moving with the velocity u ofocean currents, as in equation (4), for example as demonstratedin Jones et al. (2016). Near the surface, additional movementinduced by wind and waves may help drive motion. In general,the vertical location of oil is consequential because the e ff ectof wind and wave changes abruptly in the upper meters of theocean, as will be discussed in section 4, and because ocean cur-rents also tend to change with depth.Over the last several decades we have come to understandthe ocean as a turbulent fluid in perpetual motion, rich in tem-poral variability. There are di ff erent types of turbulence in theocean and a vast variety of instabilities—processes that can trig-ger oscillations with speeds considerably larger than the meanflow. Some of these oscillations (e.g. eddies) often result inhyperbolic points in the velocity field. Intuitively, hyperbolicpoints in the Eulerian velocity field suggest that initially-closetrajectories are likely to undergo exponential separation. Whilethis basic intuitive idea turns out to be correct, the relation be-tween the velocity at each point and the trajectories traversingthis time-dependent velocity is not as intuitive, requiring a care-ful mathematical treatment to uncover, as discussed in section5. As an example to illustrate chaotic behavior, we use HyCOM-GoM (Hybrid Coordinate Ocean Model-Gulf of Mexico), a state-of-science ocean model that is likely to be used by oil-spillmodelers in the Gulf of Mexico, to advect two groups of fourtrajectories initiated less than 5 km apart. Within each group,four trajectories are initiated within 300–500 m of each other.The two groups of trajectories undergo exponential separationas they move with realistic ocean currents, ending 200km apartafter 5 days (Fig. 1). These trajectories are initiated close to theTaylor Energy well that continued for many years to spill oilfrom the seafloor starting in 2004 (Sun et al., 2018), represent-ing a realistic example of the uncertainty that an oil-spill mod-4ler might face. Among the group of four trajectories initiatedto the southeast of Taylor well, one trajectory separates 100 kmfrom the other three trajectories, despite being initially 300500m away, further showing the chaotic nature of trajectories in theocean. This example shows that for a time-dependent flow, theinteraction of trajectories with a hyperbolic point in the velocityis complex.Exponential separation of initially close trajectories (i.e. asign of chaotic behavior) is an important part of why predictingtrajectories in the ocean is di ffi cult. Small errors in the veloc-ity field or in the location at which trajectories are initiated arelikely to result in large errors over short periods of hours todays. This is a problem inherent to ocean currents, i.e. it can-not be corrected with higher-order integration. We return to thistype of problem in section 5. Figure 1: Sea-surface velocity (black vectors) near the Mississippi delta inthe Gulf of Mexico from a HyCOM Gulf of Mexico operational simulation onMarch 14, 2016. Four trajectories (orange) are released at noon March 9, 2016,just southwest of Taylor well (red circle). Another four trajectories (blue) arereleased just southeast of Taylor well at the same time. A hyperbolic behaviorseparates initially-close ( < < A di ff erent type of problem occurs when a trajectory is de-ficient because the velocity is not representative of all the pro-cesses driving motion. In this case, it is sometimes possible toparameterize missing physics to complement the velocity. Oildoes not necessarily remain at the sea surface. The depth atwhich oil is located is particularly important when parameteriz-ing missing physics, as driving mechanisms change drasticallyeven within tens of centimeters in the upper ocean, as discussedin section 4. In this section, we provide an overview of some commonsources of ocean current data for oil spill modeling. Two prod-ucts are from remotely-sensed measurements and are thereforelimited to producing a sea-surface velocity up to current time. Ocean models produce ocean currents from numerical simula-tions but can provide velocity through the water column andcan be used to forecast ocean currents.
The equations governing geophysical fluid dynamicsfluiddynamics on a rotating spherecan be discretized and solvednumerically. The equations themselves are very complicated,and their numerical solution is further complicated because mo-tion at di ff erent spatial scales, from thousands of kilometers tocentimeters, interact with each other in fundamental ways, yetcomputers are not powerful enough to simulate all such scales.Also, simulations are necessarily initiated from imperfect ini-tial and boundary conditions, and geophysical flows tend to bechaotic. Notwithstanding, ocean models are surprisingly accu-rate in portraying a variety of physical processes in the ocean,and are often used as the source for the velocity u needed tosolve equation (4).When using an ocean model, the vertical resolution shouldbe a concern even when simulating oil transport exclusivelyat the sea surface. This is because the model’s output for asurface velocity will be in reality a representation of the ver-tically sheared currents over the height of the top grid cell ofthe model, not a representation of the velocity at the very sur-face. The model output that is closest to the surface is a depth-averaged value, where averaging takes place over the width ofthe model’s upper vertical level. Naturally, coarse resolutionsresult in greater smoothing, and therefore a less realistic repre-sentation of the surface velocity.Producing accurate velocity products with models is alsocomplicated because hyperbolic trajectories are often linked toocean instabilities that are not completely understood and aredi ffi cult to accurately simulate. For example, ocean eddies canhave a profound e ff ect on Lagrangian transport, and a numericalsimulation may develop an eddy that does not exist in the oceanbeing simulated. Even if a model accurately simulates an eddy,small displacements in the eddy location relative to the correctposition of the eddy in the ocean can result in large trajectoryerrors. Thus, ocean models often do not represent the oceanicstructures that are most influential on trajectories with enoughaccuracy. This is true even when the numerical model assimi-lates a variety of ocean measurements in an attempt to replicatethe real ocean. Observations used by data-assimilating modelsinclude satellite products such as sea-surface height, tempera-ture and salinity, and in situ observations from oceanographicbuoys, drifters, gliders, and other autonomous platforms. A re-cent overview of progress and challenges in ocean modelingcan be found in Fox-Kemper et al. (2019). High-frequency (HF) radars can measure sea-surface cur-rents remotely near the coastline ( <
200 km) with resolutionstypically 500 m to 6 km, and up to hourly in time. Velocity fromHF radar is an exponentially-weighted vertical average, with adecay scale that is proportional to the wavelength of the radarsignal (e.g. R¨ohrs et al., 2015). Due to the resolution, which is5nable to determine small-scale structures, and processing er-rors, trajectories from drifters designed to sample similar oceancurrents as those measured by HF radar, di ff er from trajectoriescomputed from HF radar velocity. Carefully calibrated radars atresolutions higher than about 1.5 km and 3 hours can replicatedrifter trajectories with a separation rate of about a few kilo-meters over a day (Rypina et al., 2014; Kirincich et al., 2012).However, the quality of HF radar processing and HF radar res-olution varies. Improving HF radar velocity products to betterrepresent coastal currents is an ongoing endeavor (e.g Kirincichet al., 2019).HF radars are an important part of operational models thatassimilate the surface velocity to minimize the model’s error.HF radar can measure currents near the coast, making theman excellent complement to altimetry that can only produce ageostrophic velocity (section 3.1.3) further from the coast.There are methods to improve HF radar data for Lagrangianpurposes. For example, to get trajectories from HF radar thatare closer to drifter trajectories, the Eulerian velocity may becorrected using trajectory data (e.g. Berta et al., 2014). A down-side of this approach is that it requires deploying drifters andallowing them to drift for some time before the corrections canbe applied.HF radar is increasingly available in the U.S. and around theworld (Roarty et al., 2019); a review on HF radar can be foundin Paduan & Washburn (2013). Satellites with altimeters are able to measure sea-surfaceheight (SSH) with enough accuracy that a geostrophic veloc-ity proportional to the SSH gradient can be computed: u g ( x , y , t ) = − gf ∂η ( x , y , t ) ∂ y (5) v g ( x , y , t ) = + gf ∂η ( x , y , t ) ∂ x (6)where ( u g , v g ) is the geostrophic surface velocity, η is theSSH, g is the acceleration of gravity and f = Ω sin θ is theCoriolis parameter, Ω = . × − s − is the rotation of theEarth, and θ is the latitude.Velocity from altimetry has been shown to give good resultsfor Lagrangian transport applications and may give superior re-sults than numerical models that assimilate the same altimetrydata. We cite a few studies as examples where using altimetricvelocity for Lagrangian transport applications has been shownto be a good choice.Ohlmann et al. (2001) compares surface velocity from driftersand altimetry and finds very good agreement in the Gulf ofMexico deeper than the 2000-m isobath and good agreementbetween the 200- and 2000-m isobaths. They find that thesecorrelations depend on the length scale over which the di ff er-entiation in equations (5) (6) is computed, with best resultsat ∂ x , ∂ y ∼ ff shorechlorophyll advection (see their figures 1 and 2); they proposea modification to the data assimilation scheme as a correction.NCOM and HyCOM are the former and current models usedby the U.S. Navy for their Global Ocean Forecast System, andthat assimilate a variety of observations through the Navy Cou-pled Ocean Data Assimilation (NCODA) system. The resultsof Olascoaga et al. (2013) and Jacobs et al. (2014) are for July2012, here we will analyze a similar transport pattern duringthe Deepwater Horizon accident in May 2010 in section 5.2.Liu et al. (2014) found that di ff erent altimetry velocity productsperform similarly and that trajectories simulated from altime-try perform better than from data-assimilative models. Bertaet al. (2015) compare satellite-tracked drifters to synthetic tra-jectories from altimetry finding “satisfactory average results”,they also show how blending drifter data into the altimetric ve-locity considerably improves trajectory hindcasts, and restoresmissing physics that cannot be explained by Ekman superpo-sition. Beron-Vera et al. (2013) and Beron-Vera et al. (2018)show the relevance of Lagrangian coherence associated witheddies detected objectively from altimetry. Essink et al. (2019)found that trajectories advected with altimetry do well in repli-cating the main transport patterns observed with drifters, andeven though they find that a variety of statistics from altimetrytrajectories do not closely resemble those from observed trajec-tories, we note that the concern for oil-spill modelers is identi-fying prevailing oil movement.Work towards remotely-sensed velocity products of higherresolution are currently underway (e.g. Chelton et al., 2019).
4. Transport in the upper layer of the ocean.
Among the challenges oil-spill modelers face is that someof the physical processes driving oil’s motion near the oceansurface may not be represented in available velocity products.In this section, we discuss how motion in the upper layer ofthe ocean may be strongly influenced by wind drift (windage)6nd Stokes drift from waves. As we will see, motion in the up-per centimeter of the ocean can be significantly di ff erent thanin the upper meter. Velocity from HF radar and ocean mod-els do not typically include Stokes drift or windage. See R¨ohrset al. (2015) for a discussion on Stokes drift in HF radar oceancurrents. Velocity from altimetry does not include Stokes drift,windage or Ekman transport, although the latter is sometimesadded from additional satellite measurements. Because it is dif-ficult to obtain a velocity that is representative of the uppercentimeter, it may be necessary to parameterize certain typesof physics if the trajectories of interest are in the order one-centimeter upper layer of the ocean. We note that near theocean’s surface, the vertical location of oil makes a big di ff er-ence; even when oil has a surface expression, large amounts ofoil may remain beneath the upper centimeter. The vertical loca-tion of oil may be determined by oil’s density (it can be heavierthan water) or due to a dynamic balance between entrainment,vertical mixing, and rise due to buoyancy (Nordam et al., 2020).Only recently have adequate observations resulted in insightinto the movement within the upper layer of the ocean. It istherefore timely to review how this information is relevant foroil spill modelers, as it suggests possibilities that were not typ-ically considered in the past. We include the Point Wells spillin the Salish Sea as a recent example where this type of physicswas used to explain the oil’s trajectory after remaining a mys-tery since 2003.Cross-shelf transport is crucial for oil-spill modeling be-cause it is needed for oil to beach, and beaching is one of themore consequential events during oil spills. A large fraction ofocean currents are in approximate geostrophic balance, stronglyconstraining flow in its ability to cross isobaths (Brink, 2016).Consequently, Lagrangian transport near the coast tends to moveparallel to the coastline (more precisely along geostrophic con-tours) and is limited in its ability to move perpendicular to thecoast (LaCasce, 2008). Among the processes capable of caus-ing cross-shelf transport there is eddying activity, and ageostrophicprocesses such as Ekman transport, Stokes drift, and windage.We also present recent evidence that windage and Stokes driftare important because they are e ff ective in driving large-scalebeaching. For some time now, it has been noted that as the wind in-creased in magnitude, the e ff ect it had on motion near the seasurface increased. Recently, Lodise et al. (2019) used data fromone of the largest Lagrangian experiments to date, to show thatundrogued drifters sampling the upper 5 cm of the ocean movewith a velocity that is 3.46.0 % of the wind under strong wind(1220 m / s) conditions, with a deflection to the right of wind di-rection increasing from 5 to 55 ◦ , as wind increased from 12to 20 m / s. For drogued drifters sampling the upper 60 cmof the ocean, the angle of deflection increased with increasingwind (12–20 m / s) from 30 to 85 ◦ , and windage ranged between2.34.1% of wind speed. In those experiments, an additionalvelocity component from Stokes drift was found to be about1.21.6% of the wind for undrogued drifters and about 0.51.2% of the wind for drogued drifters, with a deflection to the leftfrom wind’s direction of about 5 ◦ . Overall, windage and Stokesdrift accounted for about 70% of the total velocity of drogueddrifters, and about 80% of the velocity of undrogued drifters.Laxague et al. (2018) use a variety of instruments includingdi ff erent drifters to measure currents in the upper meters of theocean with an emphasis in the upper centimeters. They find thatthe velocity in the upper centimeter, about 60 cm / s, is twice thevelocity averaged over the upper meter, and four times the ve-locity averaged over the upper 10 meters. R¨ohrs & Christensen(2015) use two types of drifters, an undrogued drifter to sam-ple the surface layer and a drogued drifter to sample the upper70 cm. They find that the upper layer is influenced by wind,while subsurface motion has a stronger link to ocean dynam-ics, the result being that the surface response to wind forcingis distinct from the response 70 cm below. Androulidakis et al.(2018) also shows examples of how drogued drifters have sig-nificantly di ff erent trajectories than undrogued drifters, with thelatter heavily influenced by wind.Often the direction of motion induced by wind is not thesame as the wind direction. However, the angle of deflectioncan vary with a number of factors including the ocean’s strat-ification, the buoyancy of the particle, latitude, and the mag-nitude of the wind, making it di ffi cult to predict the directionof windage. Similar to the strategy where the eddy di ff usioncoe ffi cient is adjusted to match the observed spread of oil, itmight be necessary to use observations when possible, to ad-just the windage coe ffi cient and the angle of deflection to matchobserved motion. Some discussion on the diverse range of de-flection angles that have been measured at sea can be found inDuran (2016); windage for di ff erent objects can be found inBreivik et al. (2011); Maximenko et al. (2018). On the midnight of December 30, 2003, almost 5,000 gal-lons (about 110 barrels) of fuel spilled into the Puget Sound af-ter a tank barge accidentally overtopped near Richmond Beachin Shoreline, Washington. Helicopter overflights early in themorning observed that the surface expression of the spill driftedsouth, yet oil-spill models forecasted northward movement. Thetemporal and spatial extent of the spill was small enough (abouta day and a 15 km trajectory) for overflights to su ffi ce responsee ff orts. However, the reason why typical oil-spill model forc-ing, such as ocean currents and 1–3% of the wind, could notexplain the surface-oil trajectory remained a mystery. This wasparticularly puzzling in an enclosed sea where predictable tidesare an important component of the circulation.Published work where windage reaches 6% of the wind isunusual in the oil spill modeling community, although not un-precedented. In Duran et al. (2018b) it was conjectured thatoil’s motion during the Point Wells oil spill was driven by acombination of 6% of the wind with a 9 ◦ deflection to the rightof wind’s direction and ocean currents. This hypothesis ex-plained the trajectory of the spill that had previously been amystery, although the use of such a high windage coe ffi cientwas unusual. It wasn’t until a year later that measurements of7otion in a fine upper layer of the ocean were published byLodise et al. (2019), documenting motion dominated by windageat 6% of the wind speed.Hindcasts of the Point Wells spill advecting oil only withwind from a meteorological station in the vicinity of the spillwere suggestive for two reasons: 1) excellent agreement withthe observed trajectory in the first six hours and 2) it forced oiltowards the south, an observed-trajectory feature that had beendi ffi cult to emulate with ocean currents and typical windage(Fig. 2).When 6% of the wind from a meteorological station, andocean currents from an ocean model that replicated tides withhigh skill, were combined to force the trajectory computation,the resulting trajectory matched the correct locations at the cor-rect times throughout the spill, finally beaching at the correct lo-cation in the afternoon of December 30, 2003. A full account ofthe numerical experiments can be found in Duran et al. (2018b). Figure 2: Trajectory (orange, red circles mark locations at hourly intervals)initiated at the time and location of the Point Wells 2003 spill in the Salish sea(white star) forced only with 6% of the wind measured by the NOAA windstation (white circle with black cross near bottom). Locations and times whereoil was observed are marked with white squares and text.
Stokes drift is a net drift in the direction of wave propaga-tion caused by the asymmetrical orbital motion of particles nearthe surface induced by passing waves. Some authors convey theidea that Stokes drift is canceled in the mean due to the Coriolise ff ect. However, there is a large body of evidence suggestingcancellation in the near-surface is negligible in the presence ofturbulence induced by wind stress, which is the typical condi-tion in the ocean (e.g. Clarke & Van Gorder, 2018). Stoke’sdrift was an important driver of oil during the Deepwater Hori-zon: it was responsible for the observed beaching patterns and it is believed to have avoided oil being entrained by the LoopCurrent (Carratelli et al., 2011; Le H´ena ff et al., 2012; Weisberget al., 2017). This is consistent with other studies reporting thatStokes drift can exceed the Eulerian mean in the cross-shelf di-rection (e.g. Monismith & Fong, 2004).Stokes drift is mainly driven by high-frequency waves, thatis, waves forced by local wind rather than remote swell (D’Asaro,2014; Clarke & Van Gorder, 2018). Using many years of hourlyconcurrent wind and directional wave spectra from buoys in theGulf of Mexico and the Pacific, Clarke & Van Gorder (2018) de-rived a simple formula with which Stokes drift can be parametrizeddirectly from local wind, with good accuracy (within about 1cm / s from the average Stokes drift) for common wind speeds(between 1 to 50 m / s): u Stokes = . u ∗ ln (0 . U / u ∗ ) (7)where u Stokes is the magnitude of Stokes drift, U is the windspeed 10 meters above sea level and u ∗ = (cid:112) τ/ρ is the frictionalvelocity, the square root of wind-stress magnitude divided bya reference seawater density. The direction of Stokes drift isgiven by the unit vector in the direction of the wind. This isgood news because local wind data is often available, whetherfrom meteorological stations or operational models. Furthergood news is that this result holds even in the presence of swell.Both Clarke & Van Gorder (2018) and Onink et al. (2019)note that an additional contribution to transport at the surface,also in the direction of the wind, may be necessary due to wavebreaking. It is also possible that swell may induce Stokes driftnear the coastline, as waves become increasingly nonlinear dueto interactions with the bottom. Deeper down within the wa-ter column, Stokes drift from internal waves at the pycnoclinehas also been observed to be an e ff ective driver of oil transport(Shanks, 1987). How tracers respond in the upper ocean when they samplevelocity structures with influential vertical motion along theirpath, is an active topic of research that is receiving consider-able attention as new observational tools and experiments allowmeasuring smaller-scale processes, while higher numerical-modelresolutions become accessible. Considerable progress was madewhen such technological advances coincided with funding thatbecame available following the Deepwater Horizon accident.Identification of submesoscale structures, such as fronts andLangmuir cells, can be important because actionable counter-measures while responding to an oil spill require oil to reacha certain thickness. Intense convergence of oil along water-mass boundaries may, therefore, create an ideal location formitigation strategies, when conditions are right; for example,it seems that water subduction should persist for long enough(Gula et al., 2014). Water-mass subduction detection is there-fore suggested as an aid to identifying regions of thick oil. Itshould also be cautioned that using divergence as a diagnosticto identify regions of accumulation can lead to false positives8nd negatives (Serra et al., 2020). Clustering may happen in aregion of positive Eulerian-velocity divergence, we present anexample in section 5.1. Another potential caveat is the e ff ect ofstrong wind acting directly on buoyant oil. The experiment inRomero et al. (2019) suggests that wind order 10 m / s does notimpede strong vertical motion at fronts, although the tracer intheir experiment was neutrally buoyant, so how does the behav-ior of a positively-buoyant tracer di ff er?Types of vertical motion that are known to a ff ect the hori-zontal distribution of oil are related to ocean fronts, filaments,and Langmuir circulation. When two water masses meet, afront is formed along the boundary between the two. Whetherone water type sinks under the other because it is heavier, or be-cause of cabbeling, fronts in the ocean tend to be accompaniedby water subduction. A vertical circulation due to similar rea-sons also forms along the boundaries of filaments (McWilliams,2017). Thus, a downward vertical velocity is induced at theboundaries of water masses, which implies convergence in thehorizontal plane. Frontal regions are characterized by relativevorticity and negative divergence that are several times greaterthan planetary vorticity. This can have a profound local e ff ecton transport, collapsing floating material to essentially a point.An example of drifters originally spanning a width of about 10km, collapsing to 60 meters, can be found in D’Asaro et al.(2018). A comparison between two- and three-dimensional cir-culation at a scale of about 100 m, illustrating the distribution ofa neutrally buoyant tracer due to vertical motion, can be foundin Romero et al. (2019), showing the tracer sinking relativelyrapidly. If the tracer were buoyant, as is oil, an agglomera-tion of tracer could be expected at the surface. Androulidakiset al. (2018) found that the front-induced circulation dominatedthe trajectories of undrogued drifters, even under considerablewind, although wind may modulate their speed along the front,Also, the wind is one of the factors determining the location ofthe front. Fronts are also of interest because they may serve astransport barriers (Androulidakis et al., 2018).Langmuir circulation in its most basic form results from theinteraction of Stokes drift induced by surface waves, and thevertical shear induced by the turbulent transfer of momentumfrom wind to the upper ocean (Thorpe, 2004; Sullivan et al.,2012). The book by B¨uhler (2014) describes how a mean flow isinduced by an instability (Craik-Leibovich instability), a mech-anism that turns out to be robust and therefore explains whyLangmuir cells are ubiquitous in the ocean. As with the cir-culation associated with fronts, Langmuir circulation also hasan important vertical component, that likewise concentrates oilinto bands within minutes to hours, typically at scales of me-ters, to hundreds of meters (Chang et al., 2019; Simecek-Beatty& Lehr, 2017; D’Asaro, 2000). Convergence due to frontal cir-culation may dominate convergence due to Langmuir circula-tion (Romero et al., 2019). Changes in the vertical locationof oil droplets (e.g. McWilliams & Sullivan, 2000) induced byLangmuir circulation enhances the dispersion of oil by sub-jecting droplets to di ff erent ocean currents, as determined byvertical shear (Thorpe, 2004). Also, at least sometimes, Lang- muir circulation may be a more important part of ocean dy-namics than previously thought (D’Asaro, 2014). For exam-ple, the newest ocean climate modelsdesigned to study climatechangenow parametrize the e ff ect of Langmuir mixing at theocean’s surface. Evidence that Langmuir circulation may in-duce a large-scale coastal circulation can be found in Kukulkaet al. (2012).In the larger picture, the spatial scales of intense subduc-tion structures (few kilometers to hundreds of meters) implythat they are likely to be embedded within larger ( >
5. Modern Lagrangian tools
The regions where the separation of initially-close trajecto-ries and the attraction of initially-separate trajectories occur areregions of special interest when studying horizontal motion. Itis these regions that have an exceptional influence on the move-ment of nearby parcels, and thereby play a leading role in or-ganizing the flow into identifiable and predictable patterns. Ina two-dimensional flow, these regions are hyperbolic lines, in athree-dimensional flow they are surfaces. A rigorous approachto detecting these regions has been developed by identifying re-gions with maximal normal attraction, typically referred to ashyperbolic Lagrangian Coherent Structures (LCS; Farazmand& Haller, 2012; Haller, 2015).LCS theory builds on the concept of a flow map, a func-tion that maps every initial ( t = t ) position within a domainof interest x ∈ U , to its current ( t = t ) position x ( t ), that is:9 t t ( x ) (cid:66) x ( t ; x , t ). The Jacobian of the flow map D F t t isgiven by D F t t ( x ) (cid:66) ∂ x ∂ x ∂ x ∂ y ∂ y ∂ x ∂ y ∂ y (8)Informally, the Jacobian of the flow map (8) can be used to maptrajectory perturbations from one time to another, and this lin-ear approximation can then be used to optimize quantities ofinterest. For example, normal attraction of nearby fluid parcelsalong a trajectory over a time interval can be maximized withrespect to perturbations of the initial-time normal vector. Thisis the strategy used to find hyperbolic LCS, trajectories char-acterized by maximal normal attraction or repulsion. Workingout the math for this optimization problema formal account ofwhich can be found in Haller (2015) and references thereintheCauchy-Green (CG) strain tensor arises naturally, defined as C t t ( x ) (cid:66) (cid:104) D F t t ( x ) (cid:105) (cid:62) D F t t ( x ) . (9)In particular, the eigenvalues 0 < λ ( x ) < λ ( x ) and nor-malized eigenvectors ˆ ξ ( x ) ⊥ ˆ ξ ( x ) of (9), are used to set upordinary di ff erential equations from which hyperbolic, ellipticand parabolic LCS can be found (Haller, 2015). Thus, the CGtensor is central in LCS theory. Note that to obtain the CG ten-sor one must integrate the velocity; we return to the CG tensorin section 5.1.The mathematical formality behind LCS has proven a ver-satile approach to understanding Lagrangian motion. Hyper-bolic LCS will accurately identify how fluid will deform (i.e.along attracting hyperbolic LCS), anticipating the more influ-ential transport patterns. However, the final results ultimatelydepend on the accuracy of the velocity field, which is whatinduces Lagrangian transport in the first place. If a velocityfield is relatively accurate while having localized errors, trajec-tories traversing the time and location of errors in the velocityare likely to give wrong results relative to observed trajecto-ries. LCS will be negatively a ff ected by those velocity errors aswell, correctly identifying transport induced by the velocity, yetremaining incorrect relative to observed transport. The compu-tation of trajectories propagates localized velocity errors, oftenresulting in trajectories that are incorrect relative to observa-tions when hindcasting or forecasting transport of oil. The needto integrate the velocity field limits the suitability of observa-tional velocity data sets and of ocean models that assimilatesuch data, for Lagrangian transport purposes.Velocity products in our time can be relatively accurate thanksto relatively accurate measurements over wide areas, with satel-lite altimetry being particularly relevant because of good globalcoverage, and because it captures what is often an importantpart of the velocity at the sea surface (section 3.1.3). However,the coarse temporal and spatial resolution of altimetry meansthat the resulting velocity will almost certainly have deficien-cies. Ocean models assimilating data inherit these deficienciesand have shortcomings of their own. Given the chaotic sensitiv-ity of trajectories, it has therefore been a natural development to try to bypass the sensitivity resulting from localized velocityerrors.Another complication with hyperbolic LCS from an appliedpoint of view is that there is a timescale T involved in the com-putation. When computing (8), a choice must be made for theinitial time t and the final time t = t + T . The choice for T , sometimes even the choice for t , are often not clear a pri-ori, forcing subjective choices. Since LCS are material linesmoving with turbulent flow, these choices can result in big dif-ferences. In subsection 5.1 we describe a way to bypass sensi-tivity to the velocity field, and the need to choose T . The fundamental equation translating from the Eulerian andLagrangian characterizations of fluid flow is (4), an equivalencebetween the velocity of a parcel traversing a trajectory x ( t ), andthe Eulerian velocity u , at the parcel’s time and location. It istherefore a brilliant idea to exploit this correspondence by ex-ploring Eulerian counterparts for LCS. Serra & Haller (2016)developed Objective Eulerian Coherent Structures (OECS), in-cluding attracting hyperbolic OECS. As mentioned in section 5,the CG tensor is central to finding hyperbolic LCS, trajectoriesthat maximize normal attraction, thus maximizing the influenceon nearby water parcels and organizing flow. As shown in Serra& Haller (2016), the Taylor expansion of the CG tensor with re-spect to time is given in terms of the strain-rate tensor S : C t t ( x ) = I + S ( x , t ) ( t − t ) + O (cid:16) | t − t | (cid:17) This means that for time close enough to t , Lagrangian de-formation is given by the Eulerian strain-rate tensor, the i j -thentry of which is given by (cid:16) ∂ u i /∂ x j + ∂ u j /∂ x i (cid:17) . The strain-rate tensor has eigenvalues s , s with corresponding eigen-vectors e and e . Attracting hyperbolic OECS are tangent to e , their cores given by minima in the eigenvalue s , which isthe rate of change of the length of the normal eigenvector e due to the deformation induced by the flow; equivalently, s isthe strength of attraction normal to e . Negative values of s mean that there is attraction normal to e , the more negative thestronger the attraction. Thus, Serra & Haller (2016) extendedthe theory of LCS from finite time to their instantaneous limit interms of an Eulerian quantity, where there is no longer a need tointegrate the velocity field, e ff ectively bypassing the attendantsensitivity.The strain-rate tensor is objective, i.e. the results from com-puting Eulerian Coherent Structures are frame invariant (Haller,2015; Serra & Haller, 2016). This means that changes of ref-erence frames characterized by time-dependent rotations andtranslations will not a ff ect the results. This is important be-cause non-objective methods might give di ff erent results undercoordinate transformations, e.g. Eulerian-velocity hyperbolicpoints are not Galilean invariant (Serra & Haller, 2016).Serra et al. (2020) use attracting hyperbolic OECS, whichthey call TRAPs (TRansient Attracting Profiles), to demonstratein a series of experiments with satellite-tracked drifters, that10RAPs organize flow and perform better than trajectory com-putations in predicting drifter locations. They use a carefullycalibrated HF radar velocity and a data-assimilating model sim-ilar to what the U.S. Coast Guard would use for search and res-cue operations. Similar to hyperbolic LCS, TRAPs are linesin a two-dimensional flow; they have a core which is wherenormal attraction is maximal, with attraction strength decayingalong the rest of the TRAP. We present an example of satellite-tracked drifters converging to TRAPs computed from HF radarvelocity in Martha’s Vineyard in Massachusetts (Fig. 3). In thisexample, there is confluence of drifters at a TRAP where theEulerian horizontal velocity divergence is positive. A descrip-tion of these data, how TRAPs organize Lagrangian motion,and how TRAPs can be used for search and rescue operations,can be found in Serra et al. (2020). In the next section (sec-tion 5.2), we present an example where TRAPs can predict themovement of oil at least 8 days in advance, while trajectoriesdiverge from the observed transport due to a likely-erroneoushyperbolic point. Figure 3: TRAPs (red lines) and their cores (red circles) computed from HFradar velocity o ff Martha’s Vineyard, MA, plotted over the velocity divergence(color contours; day − ) with satellite-tracked drifters (white dots) convergingto TRAPs. Black lines are streamlines. TRAP A is in a region of positivehorizontal velocity divergence. TRAPs remain invisible to divergence fieldsand streamlines. The di ffi culty of simulating Lagrangian transport can beeasily experienced by trying to replicate observed trajectories.During the Deepwater Horizon, at least six di ff erent ocean mod-els were used in an attempt to forecast the location of oil, to pro-vide critical information for response and planning (Liu et al.,2011). However, there was enough intermodel variability thatensemble averaging was recommended to produce a forecastthat was more likely to occur. Even then, forecasts were limitedto two days due to forecast error growth. We note that these were ideal conditions as often there aren’t that many oceanmodels available for ensemble averaging.Here we present a di ff erent method (section 5.1) that seeksto bypass the sensitivity that causes error growth, using TRAPswhich are computed from instantaneous snapshots of an Eule-rian velocity. We show that an analysis combining TRAPs andLCS is enough to 1) accurately forecast the observed movementof oil at least 8 days in advance, and 2) understand why the sim-ulated Lagrangian transport does not conform to observations.In this example, forecasts only depend on a previous-day sin-gle velocity snapshot, and the LCS are computed from onlypast information to complement the information obtained fromTRAPs and pinpoint the source of error in the velocity field,and its Lagrangian manifestation.A variety of velocity products from altimetry are available,some of them including an Ekman component. The product thatwe use here is a daily velocity by GEKCO2 (Geostrophic andEKman Current Observatory; Sudre et al., 2013), from satellitealtimetry and wind. We confirm our results using a daily instan-taneous velocity from HyCOM Global at about 9 km resolutionin the Gulf of Mexico, the current U.S. Navy Operational model(Burnett et al., 2014) that assimilates a variety of observationsusing the Navy Coupled Ocean Data Assimilation (NCODA)and that is forced with the NAVy Global Environmental Model(NAVGEM). Transport simulated with HyCOM Gulf of Mex-ico, which has a similar setup as HyCOM Global, but at a 4 kmresolution, produces similar simulated transport as the other ve-locity products used in this experiment.We analyze daily forecasts of the movement of oil duringthe Deepwater Horizon accident by comparing TRAPS to theobserved outline of oil and the advection of oil obtained by in-tegrating the velocity field, i.e. computing trajectories betweenMay 1117, 2010, that are initiated at the location of oil observedon May 10, 2010. TRAPs are computed from snapshots of thevelocity on previous days, thus providing forecasts within thishindcast exercise.The forecast on May 10 shows a weak TRAP (near 28.2N,88.3W) suggesting slight oil movement towards the southeast,coinciding with the outline of oil observed on May 11, and withtrajectories computed between May 10–11, although the simu-lated oil and the TRAP are slightly o ff set from the observedoil (Fig. 4). The strength of attraction of the TRAP is low(about 0.3 day − ), accurately forecasting slight oil movement.The only other TRAP core in contact with the observed oil (near29.1N, 88W), is the core of a TRAP almost entirely containedwithin the observed oil on May 11, and therefore cannot be ex-pected to cause a significant rearrangement of oil outside of theobserved oil outline. TRAPs in the southeast section of figure 4have higher strengths of attraction of about 1 day − , but are stillrelatively far from the oil.11 igure 4: The blue line is the outline of oil as observed from satellites on May10, 2010, used as initial conditions for trajectory computations. In orange arethe final positions (May 11) of trajectories initiated within the blue outline onMay 10, computed by integrating GEKCO2 velocity. The black line is theoutline of oil as observed on May 11, 2010. Black vectors are the velocity fromGEKCO2 on May 10, and the red lines are TRAPs computed with the velocityon May 10, TRAP cores (colored circles) are colored according to attractionstrength (day − ; colorscale on the right). By May 13, oil trajectories continued along the path to-wards the southwest, the weak TRAP computed from the May12 velocity accurately forecasting that path (Fig. 5). This willbe the last day this weakly-attracting TRAP appears near 28.2N,88.3W. The same TRAPs with strong attraction computed withthe May 10 velocity remain when TRAPs are computed withthe May 12 velocity (Figs. 4 and 5). These TRAPs in thesoutheast of our domain continue to have the strongest attrac-tion, three to four times stronger than TRAPs directly interact-ing with oil.
Figure 5: Same as in figure 4 but on May 13, 2010.
By May 15, observed oil has aligned with one of the strongest-attraction TRAPs, the one that has remained at the same lo-cation since May 10 near 27.4N and 87.25W. Meanwhile, thesimulated oil trajectory continues its original path towards thesouthwest, by now clearly diverging from the observed path(Fig. 6). The weak TRAP originally indicating the path to-wards the southwest (near 28.2N, 88.3W in Figs. 4 and 5) is no longer present in the May 14 velocity, and will not be seenagain during the rest of our analysis.
Figure 6: Same as in figure 4 but on May 15, 2010.
By May 17, observed oil has deformed towards the souththen east, while the simulated oil trajectory has deformed to-wards the south then west, thus the observed and simulated tra-jectories are heading in opposite directions (Fig. 7). LCS il-lustrate a hyperbolic point near 28.4N and 87.7W where trans-port splits, the western part heading towards the south then west(simulated oil follows this LCS) and the eastern part movingsouth then east (observed oil follows this LCS). Thus, LCSshow that the simulated tracer just barely missed the observedtransport pattern that is accurately depicted by the LCS on theeastern side of the hyperbolic point. Note there is a TRAPabove the LCS on the eastern side of the hyperbolic point near28.4N and 87.7W, accurately selecting the altimetric LCS thatagrees with the path of observed oil transport. Further confir-mation comes from the strong TRAP near 27.4N and 87.25Wforecasting elongation of oil along that direction since at leastMay 10. This shows how TRAPs and LCS provide complemen-tary information, together explaining the discrepancy betweensimulated and observed trajectories. With such a detailed char-acterization of the available information, an oil-spill modelercan then identify which patterns are most likely to occur andwhich patterns are likely spurious. In this example, a TRAP isable to predict the observed movement of oil at least 8 days inadvance, starting from the velocity snapshot on May 10, 2010,while the simulated trajectory for oil initiated from the observedoil on May 10, 2010, is caught on the wrong side of a hyperbolicpoint and ends up moving in the opposite direction relative toobserved transport.12 igure 7: Same as in figure 4 but on May 17, 2010; purple lines are attractingLCSs computed back in time between May 17 and May 10.
Although initially the velocity is correct in inducing trans-port towards the southwest, the TRAP that accurately identifiessouthwest motion is weak and it disappears after a few days(Figs. 4, 5 and 6). The hyperbolic point causing the diver-gence of simulated transport relative to observed transport (Fig.7) therefore seems to originate from a disparity in the veloc-ity arising from coarse temporal resolution and resulting in er-ror accumulation during Lagrangian integration. The problemmay be related to the period between passes of altimetry satel-lites being too long to capture changes in the ocean velocity ontimescales of a few days, and the attendant influence on trajec-tories when integrating the velocity. Fortunately, velocity fromaltimetry accurately captures the features that result in the maintransport patterns, it is just that velocity integration is not theadequate tool to extract this information.The above results are from the GEKCO2 velocity; very sim-ilar results are obtained using HyCOM Global (not shown).Simulated transport is also very similar when using HyCOMGoM (not shown). The computation of TRAPs for flows athigh resolutions (about 4 km or less in the ocean) may requirefiltering and is a topic of current research. Global GEKCO2velocity is available since 1993 to date minus two days.
6. Conclusion and Outlook
Recent advances suggest that better results for oil-spill mod-elers is a reachable goal. In this chapter, we have shown exam-ples of how a basic understanding of the physics driving mo-tion in the sea and the use of novel Lagrangian and EulerianCoherent Structures techniques can result in improved oil-spillmodeling.Recent progress in Coherent Structures techniques— com-puting attracting structures that shape material transport froman Eulerian snapshot—is a promising development for oil-spillmodeling e ff orts. We have shown how TRAPs (or AttractingOECSs) are able to bypass errors in the velocity that producelarge errors in simulated trajectories.Olascoaga & Haller (2012) explored similar ideas by search-ing for the most persistent hyperbolic cores using 15-day inte- grations. Lagrangian integration over such a period acts as afilter, removing short-term variabilities and focusing on finite-time mesoscale features. Our results are consistent with theirs:a highly attractive hyperbolic core persists for over a week, ac-curately anticipating prominent fluid deformation. The advan-tages of TRAPs are that they do not require velocity integra-tion, they can be computed from a single velocity snapshot, andthey predict hyperbolic attraction cores whether persistent orephemeral.As described in section 5.2, TRAPs were able to identifythe correct transport patterns, while LCSs were influenced bya deficient velocity, as were the simulated trajectories. The er-roneous simulated transport is initially correct as evidenced bythe observed movement of oil, aptly identified in the velocity bya weak, ephemeral TRAP. However, simulated trajectories be-come erroneous as integration causes velocity-error accumula-tion, while a strong persistent TRAP marked the correct regionof oil confluence well in advance of observed deformation.Higher resolution observations of sea-surface velocity andsurface processes are expected to advance considerably our un-derstanding, ultimately resulting in improved velocity products.For example, using high-resolution observations, theoretical mod-eling, and coupled ocean-atmosphere-wave numerical modelsto study the coupled e ff ects of wind, upper ocean processes, andwind waves can be expected to significantly improve our under-standing of the ocean’s surface (Villas Bˆoas et al., 2019). Im-provements in observations and understanding should translateto improved oil transport forecasts. As velocity products im-prove by including more of the physics relevant to simulatingoil’s movement, the techniques highlighted here will becomemore relevant. A basic understanding of ocean physics willcontinue to be needed to supplement velocity products lack-ing certain types of forcing, but also to understand new velocityproducts that will incorporate more types of physics than previ-ously available.Despite improvements in velocity products, the sensitivityof trajectory computations to small errors will likely continueto produce erroneous trajectories. This suggests that the useof novel techniques that seek to bypass the sensitivity inher-ent to trajectory computations are likely to become importanttools for the oil-spill modeler. Thus, as an e ff ort parallel to im-proving velocity products, progress in techniques bypassing theproblems inherent to the unstable nature of ocean currents canbe expected. Examples include Objective Eulerian CoherentStructures for instantaneous transport patterns and climatolog-ical Lagrangian Coherent Structures for climatological trans-port patterns. The latter is an empirical approach developedin Duran et al. (2018a) where it was found that the filtering ofthe velocity by computing a climatology is surprisingly accu-rate for identifying recurrent Lagrangian transport patterns ifthe proper Lagrangian tools are used to extract the transportpatterns. Among their results, the transport pattern studied insection 5.2 turns out to be a recurrent pattern, and therefore apattern that is likely to be seen in May through August of anygiven year. A climatological approach should not replace fore-13asts, yet it does provide a valuable general understanding ofpersistent transport barriers, trajectories, regions of persistentattraction, and persistent isolation. Thus, Lagrangian climatolo-gies in the sense of Duran et al. (2018a) complement the inter-pretation of forecasts while providing a general understandingof Lagrangian motion in a region of interest. The climatologicalapproach suggests that progress can be made by understandingthe connection between the inherently time-dependent trajecto-ries of an instantaneous ocean velocity, and a low-pass filteredclimatological velocity. An alternative approach to understand-ing uncertainty in oil-spill modeling is the use of ensemble sim-ulations to create a surrogate model (Zhang et al., 2020). It ispossible that future work might be able to bridge the Lagrangianclimatology strategy of Duran et al. (2018a) with the ensemble-based surrogate model of Zhang et al. (2020).As ocean observations and ocean models improve with newsatellite products, a greater number of HF radars, drifters andautonomous vehicles, and advances in data assimilation and nu-merical modeling, oil-spill modelers should be able to capital-ize from the material presented here, achieving a higher rate ofsuccess in forecasting and hindcasting the movement of oil.Accurately simulating Lagrangian transport in the ocean isof considerable societal interest for a variety of reasons includ-ing oil spills, the fate of other contaminants, fisheries, oceanecology, search and rescue, tracing accidents or crimes back intime (forensic work), climate change and weather predictions,among others. Many countries have conducted oceanographicresearch for several decades now. Consequently, enough progresshas been made to where simulating trajectories in the ocean of-ten produces valuable information. For the needed progress tocontinue, we must understand the ocean’s importance for soci-ety at large, and that the relevance of oceanographic endeavoris increasing due to pressing issues such as climate change,coastal development, population growth, and globalization.
7. Acknowledgments
The GEKCO2 product used in this study was developed andextracted by Jol Sudre et al. at LEGOS, France. GEKCO2data can be requested at . Funding for the developmentof HyCOM has been provided by the National Ocean Partner-ship Program and the O ffi ce of Naval Research. Data assimila-tive products using HyCOM are funded by the U.S. Navy. Com-puter time was made available by the DoD High PerformanceComputing Modernization Program. The output is publiclyavailable at https: // hycom.org. Deepwater Horizon oil productsfrom NOAA are available at . RD would like to thank M. J. Olasco-aga for helpful conversations.The work of RD was performed in support of the US De-partment of Energys Fossil Energy, Oil and Natural Gas Re-search Program. It was executed by NETLs Research and In-novation Center, including work performed by Leidos Research Support Team sta ff under the RSS contract 89243318CFE000003.M.S. acknowledges support from the Schmidt Science Fellow-ship and the Postdoc Mobility Fellowship from the Swiss Na-tional Foundation.This work was funded by the Department of Energy, Na-tional Energy Technology Laboratory, an agency of the UnitedStates Government, through a support contract with Leidos Re-search Support Team (LRST). Neither the United States Gov-ernment nor any agency thereof, nor any of their employees,nor LRST, nor any of their employees, makes any warranty,expressed or implied, or assumes any legal liability or respon-sibility for the accuracy, completeness, or usefulness of anyinformation, apparatus, product, or process disclosed, or rep-resents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process,or service by trade name, trademark, manufacturer, or other-wise, does not necessarily constitute or imply its endorsement,recommendation, or favoring by the United States Governmentor any agency thereof. The views and opinions of authors ex-pressed herein do not necessarily state or reflect those of theUnited States Government or any agency thereof. Bibliography
Androulidakis, Y., Kourafalou, V., ¨Ozg¨okmen, T., Garcia-Pineda, O., Lund,B., Le H´ena ff , M., Hu, C., Haus, B. K., Novelli, G., Guigand, C., Kang,H. S., Hole, L., & Horstmann, J. (2018). Influence of river-induced frontson hydrocarbon transport: A multiplatform observational study. Journal ofGeophysical Research: Oceans , . doi: .Beron-Vera, F. J., & LaCasce, J. H. (2016). Statistics of simulated and observedpair separations in the Gulf of Mexico.
Journal of Physical Oceanography ,. doi: . arXiv:1505.03475 .Beron-Vera, F. J., Olascoaga, M. J., Wang, Y., Trianes, J., & P´erez-Brunius, P. (2018). Enduring Lagrangian coherence of a Loop Cur-rent ring assessed using independent observations. Scientific Reports , , 11275. URL: https://doi.org/10.1038/s41598-018-29582-5 .doi: .Beron-Vera, F. J., Wang, Y., Olascoaga, M. J., Goni, G. J., & Haller, G. (2013).Objective detection of oceanic eddies and the agulhas leakage. Journal ofPhysical Oceanography , . doi: .Berta, M., Bellomo, L., Magaldi, M. G., Gri ff a, A., Molcard, A., Marmain,J., Borghini, M., & Taillandier, V. (2014). Estimating Lagrangian transportblending drifters with HF radar data and models: Results from the TOSCAexperiment in the Ligurian Current (North Western Mediterranean Sea). Progress in Oceanography , . doi: .Berta, M., Gri ff a, A., Magaldi, M. G., ¨Ozg¨okmen, T. M., Poje, A. C., Haza,A. C., & Josefina Olascoaga, M. (2015). Improved surface velocity andtrajectory estimates in the Gulf of Mexico from Blended satellite altimetryand drifter data. Journal of Atmospheric and Oceanic Technology , . doi: .Breivik, Ø., Allen, A. A., Maisondieu, C., & Roth, J. C. (2011). Wind-induceddrift of objects at sea: The leeway field method.
Applied Ocean Research , , 100–109.Brink, K. (2016). Cross-Shelf Exchange. Annual Review of Marine Science , .doi: .B¨uhler, O. (2014).
Waves and Mean Flows . Cambridge Monographson Mechanics (2nd ed.). Cambridge University Press. doi: .Burnett, W., Harper, S., Preller, R., Jacobs, G., & LaCroix, K. (2014). Overviewof operational ocean forecasting in the us navy: Past, present, and future.
Oceanography , . URL: https://doi.org/10.5670/oceanog.2014.65 .Carratelli, E. P., Dentale, F., & Reale, F. (2011). On the E ff ects of Wave-Induced Drift and Dispersion in the Deepwater Horizon Oil Spill. In Mon-itoring and Modeling the Deepwater Horizon Oil Spill: A Record BreakingEnterprise . Wiley. doi: . hang, H., Huntley, H. S., Kirwan, A. D., Carlson, D. F., Mensa, J. A., Mehta,S., Novelli, G., ¨Ozg¨okmen, T. M., Fox-Kemper, B., Pearson, B., Pearson,J., Harcourt, R. R., & Poje, A. C. (2019). Small-scale dispersion in thepresence of langmuir circulation. Journal of Physical Oceanography , .URL: https://doi.org/10.1175/JPO-D-19-0107.1/ . doi: .Chelton, D. B., Schlax, M. G., Samelson, R. M., Farrar, J. T., Molemaker,M. J., McWilliams, J. C., & Gula, J. (2019). Prospects for future satelliteestimation of small-scale variability of ocean surface velocity and vorticity.
Progress in Oceanography , . doi: .Clarke, A. J., & Van Gorder, S. (2018). The Relationship of Near-SurfaceFlow, Stokes Drift and the Wind Stress.
Journal of Geophysical Research:Oceans , . doi: .Csanady, G. (1973).
Turbulent di ff usion in the environment . Dordrecht, Hol-land: D. Reidel Publishing Company.Dagestad, K.-F., & R¨ohrs, J. (2019). Prediction of ocean surface trajectoriesusing satellite derived vs. modeled ocean currents. Remote sensing of envi-ronment , , 130–142.D’Asaro, E. (2000). Simple Suggestions for Including Vertical Physics inOil Spill Models. Spill Science & Technology Bulletin , . doi: .D’Asaro, E. A. (2014). Turbulence in the Upper-Ocean MixedLayer.
Annual Review of Marine Science , . doi: .D’Asaro, E. A., Shcherbina, A. Y., Klymak, J. M., Molemaker, J., Novelli, G.,Guigand, C. M., Haza, A. C., Haus, B. K., Ryan, E. H., Jacobs, G. A., Hunt-ley, H. S., Laxague, N. J., Chen, S., Judt, F., McWilliams, J. C., Barkan, R.,Kirwan, A. D., Poje, A. C., & ¨Ozg¨okmen, T. M. (2018). Ocean convergenceand the dispersion of flotsam.
Proceedings of the National Academy of Sci-ences of the United States of America , . doi: .Davidson, P. A. (2015).
Turbulence: An introduction for scientists and engi-neers . (2nd ed.). Oxford, UK: Oxford University Press.Duran, R. (2016).
Sub-Grid Parameterizations for Oceanic Oil-Spill Simulations . Technical Report EPAct Technical ReportSeries; U.S. Department of Energy, National Energy Technol-ogy Laboratory. URL: https://edx.netl.doe.gov/dataset/sub-grid-parameterizations-for-oceanic-oil-spill-simulations .doi: .Duran, R., Beron-Vera, F. J., & Olascoaga, M. J. (2018a). Extracting quasi-steady Lagrangian transport patterns from the ocean circulation: An ap-plication to the Gulf of Mexico.
Scientific Reports , . doi: .Duran, R., Romeo, L., Whiting, J., Vielma, J., Rose, K., Bunn, A., & Bauer,J. (2018b). Simulation of the 2003 Foss Barge - Point Wells Oil Spill: AComparison between BLOSOM and GNOME Oil Spill Models.
Journal ofMarine Science and Engineering , . doi: .Durran, D. R. (2010).
Numerical Methods for Fluid Dynamics - With Applica-tions to Geophysics . Springer. doi: .Essink, S., Hormann, V., Centurioni, L. R., & Mahadevan, A. (2019). CanWe Detect Submesoscale Motions in Drifter Pair Dispersion?
Journal ofPhysical Oceanography , . doi: .Farazmand, M., & Haller, G. (2012). Computing lagrangian coherent struc-tures from their variational theory.
Chaos: An Interdisciplinary Journal ofNonlinear Science , , 013128.Fox-Kemper, B., Adcroft, A., B¨oning, C. W., Chassignet, E. P., Curchitser, E.,Danabasoglu, G., Eden, C., England, M. H., Gerdes, R., Greatbatch, R. J.,Gri ffi es, S. M., Hallberg, R. W., Hanert, E., Heimbach, P., Hewitt, H. T.,Hill, C. N., Komuro, Y., Legg, S., Sommer, J. L., Masina, S., Marsland, S. J.,Penny, S. G., Qiao, F., Ringler, T. D., Treguier, A. M., Tsujino, H., Uotila,P., & Yeager, S. G. (2019). Challenges and prospects in ocean circulationmodels. doi: .Gula, J., Molemaker, M. J., & McWilliams, J. C. (2014). SubmesoscaleCold Filaments in the Gulf Stream. Journal of Physical Oceanography , , 2617–2643. URL: https://doi.org/10.1175/JPO-D-14-0029.1 .doi: .Haller, G. (2015). Langrangian Coherent Structures. AnnualReview of Fluid Mechanics , , 140906185740003. URL: . doi: .Jacobs, G. A., Bartels, B. P., Bogucki, D. J., Beron-Vera, F. J., Chen, S. S., Coelho, E. F., Curcic, M., Gri ff a, A., Gough, M., Haus, B. K., Haza, A. C.,Helber, R. W., Hogan, P. J., Huntley, H. S., Iskandarani, M., Judt, F., Kirwan,A. D., Laxague, N., Valle-Levinson, A., Lipphardt, B. L., J. Mariano, A.,Ngodock, H. E., Novelli, G., Olascoaga, M. J., ¨Ozg¨okmen, T. M., Poje,A. C., Reniers, A. J., Rowley, C. D., Ryan, E. H., Smith, S. R., Spence,P. L., Thoppil, P. G., & Wei, M. (2014). Data assimilation considerations forimproved ocean predictability during the Gulf of Mexico Grand LagrangianDeployment (GLAD). Ocean Modelling , . doi: .Jones, C. E., Dagestad, K. F., Breivik, Ø., Holt, B., R¨ohrs, J., Christensen,K. H., Espeseth, M., Brekke, C., & Skrunes, S. (2016). Measurement andmodeling of oil slick transport.
Journal of Geophysical Research: Oceans , .doi: .Kim, S., Samelson, R. M., & Snyder, C. (2011). Toward an uncertainty bud-get for a coastal ocean model.
Monthly Weather Review , . doi: .Kirincich, A., Emery, B., Washburn, L., & Flament, P. (2019). Improving sur-face current resolution using direction finding algorithms for multiantennahigh-frequency radars.
Journal of Atmospheric and Oceanic Technology , .doi: .Kirincich, A. R., De Paolo, T., & Terrill, E. (2012). Improving HF radar esti-mates of surface currents using signal quality metrics, with application to theMVCO high-resolution radar system.
Journal of Atmospheric and OceanicTechnology , . doi: .Kloeden, P. E., & Platen, E. (1992).
Numerical Solution of Stochastic Di ff eren-tial Equations . Springer-Verlag Berlin Heidelberg.Kukulka, T., Plueddemann, A. J., & Sullivan, P. P. (2012). Nonlocal transportdue to Langmuir circulation in a coastal ocean. Journal of Geophysical Re-search: Oceans , . URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2012JC008340 . doi: .LaCasce, J. H. (2008). Statistics from Lagrangian observations. doi: .Lawrence, G. A., Ashley, K. I., Yonemitsu, N., & Ellis, J. R. (1995). Naturaldispersion in a small lake. Limnology and Oceanography , , 1519–1526.Laxague, N. J., ¨Ozg¨okmen, T. M., Haus, B. K., Novelli, G., Shcherbina, A.,Sutherland, P., Guigand, C. M., Lund, B., Mehta, S., Alday, M., & Mole-maker, J. (2018). Observations of Near-Surface Current Shear Help De-scribe Oceanic Oil and Plastic Transport. Geophysical Research Letters , .doi: .Le H´ena ff , M., Kourafalou, V. H., Paris, C. B., Helgers, J., Aman, Z. M., Hogan,P. J., & Srinivasan, A. (2012). Surface evolution of the deepwater horizonoil spill patch: Combined e ff ects of circulation and wind-induced drift. En-vironmental Science and Technology , . doi: .Liu, Y., Weisberg, R. H., Hu, C., & Zheng, L. (2011). Trajectory Forecast asa Rapid Response to the Deepwater Horizon Oil Spill. In
Monitoring andModeling the Deepwater Horizon Oil Spill: A Record Breaking Enterprise (pp. 153–165). Wiley. doi: .Liu, Y., Weisberg, R. H., Vignudelli, S., & Mitchum, G. T. (2014). Evaluationof altimetry-derived surface current products using Lagrangian drifter tra-jectories in the eastern Gulf of Mexico.
Journal of Geophysical Research:Oceans , . doi: .Lodise, J., ¨Ozg¨okmen, T., Gri ff a, A., & Berta, M. (2019). Vertical structureof ocean surface currents under high winds from massive arrays of drifters. Ocean Science , . doi: .Lund, B., Haus, B. K., Horstmann, J., Graber, H. C., Carrasco, R., Lax-ague, N. J., Novelli, G., Guigand, C. M., & ¨Ozg¨okmen, T. M. (2018).Near-surface current mapping by shipboard marine X-band radar: A vali-dation.
Journal of Atmospheric and Oceanic Technology , . doi: .Lynch, D. R., Greenberg, D. A., Bilgili, A., McGillicuddy Jr, D. J., Manning,J. P., & Aretxabaleta, A. L. (2014).
Particles in the coastal ocean: Theoryand applications . Cambridge University Press.Maximenko, N., Hafner, J., Kamachi, M., & MacFadyen, A. (2018). Numericalsimulations of debris drift from the great japan tsunami of 2011 and theirverification with observational reports.
Marine pollution bulletin , , 5–25.McWilliams, J. C. (2017). Submesoscale surface fronts and filaments: sec-ondary circulation, buoyancy flux, and frontogenesis. Journal of Fluid Me-chanics , , 391–432. doi: .McWilliams, J. C. (2019). A survey of submesoscale currents. GeoscienceLetters , . doi: . cWilliams, J. C., & Sullivan, P. P. (2000). Vertical Mixing by Lang-muir Circulations. Spill Science and Technology Bulletin , . doi: .Monismith, S. G., & Fong, D. A. (2004). A note on the potential transport ofscalars and organisms by surface waves.
Limnology and Oceanography , .doi: .Moum, J. N., & Rippeth, T. P. (2009). Do observations adequately resolvethe natural variability of oceanic turbulence?
Journal of Marine Systems , .doi: .Murthy, C. (1976). Horizontal di ff usion characteristics in lake ontario. Journalof physical oceanography , , 76–84.Nordam, T., Br¨onner, U., Skancke, J., Nepstad, R., Rønningen, P., & Alver, M.(2018). Numerical integration and interpolation in marine pollutant transportmodelling. In Proceedings of the Forty-first AMOP Technical Seminar .Nordam, T., & Duran, R. (2020). Numerical integrators for la-grangian oceanography.
Geoscientific Model Development Discus-sions , , 1–30. URL: https://gmd.copernicus.org/preprints/gmd-2020-154/ . doi: .Nordam, T., Skancke, J., Duran, R., & Barker, C. (2020). Vertical transport inoil-spill modeling. In Marine Hydrocarbon Spill Assessments . Elsevier.Ohlmann, J. C., Niiler, P. P., Fox, C. A., & Leben, R. R. (2001). Eddy en-ergy and shelf interactions in the Gulf of Mexico.
Journal of GeophysicalResearch: Oceans , . doi: .Okubo, A. (1971). Oceanic di ff usion diagrams. Deep sea research and oceano-graphic abstracts , , 789–802.Okubo, A., & Levin, S. A. (2013). Di ff usion and ecological problems: Modernperspectives . New York Berlin Heidelberg: Springer-Verlag.Olascoaga, M. J., Beron-Vera, F. J., Haller, G., Tri ∼ nanes, J., Iskandarani, M.,Coelho, E. F., Haus, B. K., Huntley, H. S., Jacobs, G., Kirwan, A. D., Lip-phardt, B. L., ¨Ozg¨okmen, T. M., H. M. Reniers, A. J., & Valle-Levinson,A. (2013). Drifter motion in the Gulf of Mexico constrained by alti-metric Lagrangian coherent structures. Geophysical Research Letters , .doi: .Olascoaga, M. J., & Haller, G. (2012). Forecasting sudden changes in environ-mental pollution patterns.
Proceedings of the National Academy of Sciencesof the United States of America , . doi: .Onink, V., Wichmann, D., Delandmeter, P., & van Sebille, E. (2019). The Roleof Ekman Currents, Geostrophy, and Stokes Drift in the Accumulation ofFloating Microplastic.
Journal of Geophysical Research: Oceans , . doi: .Paduan, J. D., & Washburn, L. (2013). High-Frequency RadarObservations of Ocean Surface Currents.
Annual Reviewof Marine Science , , 115–136. URL: https://doi.org/10.1146/annurev-marine-121211-172315 . doi: .Richardson, L. F. (1926). Atmospheric di ff usion shown on a distance-neighbourgraph. Proceedings of the Royal Society of London, Series A , , 709–737.Richardson, L. F., & Stommel, H. (1948). Note on eddy di ff usion in the sea. Journal of Meteorology , , 238–240.Roarty, H., Cook, T., Hazard, L., George, D., Harlan, J., Cosoli, S., Wy-att, L., Alvarez Fanjul, E., Terrill, E., Otero, M., Largier, J., Glenn,S., Ebuchi, N., Whitehouse, B., Bartlett, K., Mader, J., Rubio, A.,Corgnati, L., Mantovani, C., Gri ff a, A., Reyes, E., Lorente, P., Flores-Vidal, X., Saavedra-Matta, K. J., Rogowski, P., Prukpitikul, S., Lee,S.-H., Lai, J.-W., Guerin, C.-A., Sanchez, J., Hansen, B., & Grilli, S.(2019). The Global High Frequency Radar Network. Frontiers in MarineScience , , 164. URL: . doi: .R¨ohrs, J., & Christensen, K. H. (2015). Drift in the uppermost part of the ocean. Geophysical Research Letters , . doi: .R¨ohrs, J., Sperrevik, A. K., Christensen, K. H., Brostr¨om, G., & Breivik,Ø. (2015). Comparison of HF radar measurements with Eulerianand Lagrangian surface currents.
Ocean Dynamics , . doi: .Romero, L., Ohlmann, J. C., Pall`as-Sanz, E., Statom, N. M., P´erez-Brunius,P., & Maritorena, S. (2019). Coincident observations of dye and drifterrelative dispersion over the inner shelf.
Journal of Physical Oceanography ,. doi: .Rypina, I. I., Kirincich, A. R., Limeburner, R., & Udovydchenkov, I. A.(2014). Eulerian and Lagrangian correspondence of high-frequency radarand surface drifter data: E ff ects of radar resolution and flow compo- nents. Journal of Atmospheric and Oceanic Technology , . doi: .Salmon, R. (1998).
Lectures on Geophysical Fluid Dynamics . Oxford Univer-sity.Serra, M., & Haller, G. (2016). Objective eulerian coherent structures.
Chaos ,. doi: . arXiv:1512.02112 .Serra, M., Sathe, P., Rypina, I., Kirincich, A., Ross, S. D., Lermusiaux, P.,Allen, A., Peacock, T., & Haller, G. (2020). Search and rescue at seaaided by hidden flow structures. Nature Communications , . doi: . arXiv:1909.07828 .Shanks, A. L. (1987). The onshore transport of an oil spill by internal waves. Science , , 1198–1200. URL: https://science.sciencemag.org/content/235/4793/1198 . doi: .Simecek-Beatty, D., & Lehr, W. J. (2017). Extended oil spill spreadingwith Langmuir circulation. Marine Pollution Bulletin , . doi: .Smagorinsky, J. (1963). General circulation experiments with the primitiveequations: I. the basic experiment.
Monthly weather review , , 99–164.Smyth, W. D., & Carpenter, J. R. (2019). Instability in Geophysical Flows .Cambridge University Press.Spivakovskaya, D., Heemink, A. W., & Deleersnijder, E. (2007). Lagrangianmodelling of multi-dimensional advection-di ff usion with space-varying dif-fusivities: theory and idealized test cases. Ocean Dynamics , , 189–203.Sudre, J., Maes, C., & Garc¸on, V. (2013). On the global estimates ofgeostrophic and Ekman surface currents. Limnology and Oceanography:Fluids and Environments , . doi: .Sullivan, P. P., Romero, L., Mcwilliams, J. C., & Kendall Melville, W. (2012).Transient evolution of langmuir turbulence in ocean boundary layers drivenby hurricane winds and waves.
Journal of Physical Oceanography , . doi: .Sun, S., Hu, C., Garcia-Pineda, O., Kourafalou, V., Le H´ena ff , M., &Androulidakis, Y. (2018). Remote sensing assessment of oil spillsnear a damaged platform in the Gulf of Mexico. Marine Pollu-tion Bulletin , , 141–151. URL: . doi: https://doi.org/10.1016/j.marpolbul.2018.09.004 .Tennekes, H., & Lumley, J. L. (1972). A first course in turbulence . MIT press.Thorpe, S. A. (2004). LANGMUIR CIRCULATION.
Annual Review ofFluid Mechanics , , 55–79. doi: .Villas Bˆoas, A. B., Ardhuin, F., Ayet, A., Bourassa, M. A., Brandt, P., Chapron,B., Cornuelle, B. D., Farrar, J. T., Fewings, M. R., Fox-Kemper, B., Gille,S. T., Gommenginger, C., Heimbach, P., Hell, M. C., Li, Q., Mazlo ff , M. R.,Merrifield, S. T., Mouche, A., Rio, M. H., Rodriguez, E., Shutler, J. D., Sub-ramanian, A. C., Terrill, E. J., Tsamados, M., Ubelmann, C., & van Sebille,E. (2019). Integrated Observations of Global Surface Winds, Currents, andWaves: Requirements and Challenges for the Next Decade. Frontiers inMarine Science , . doi: .Wagner, P., R¨uhs, S., Schwarzkopf, F. U., Koszalka, I. M., & Biastoch, A.(2019). Can Lagrangian tracking simulate tracer spreading in a high-resolution ocean general circulation model?
Journal of Physical Oceanog-raphy , . doi: .Weisberg, R. H., Lianyuan, Z., & Liu, Y. (2017). On the movement ofDeepwater Horizon Oil to northern Gulf beaches.
Ocean Modelling , .doi: .Zelenke, B., O’Connor, C., Barker, C., Beegle-Krause, C., & Eclipse, L. E.(2012).
General NOAA Operational Modeling Environment (GNOME),Technical Documentation, NOAA Technical Memorandum NOS OR & R 40 .Technical Report October NOAA.Zhang, R., Wingo, P., Duran, R., Rose, K., Bauer, J., & Ghanem,R. (2020). Environmental economics and uncertainty: Reviewand a machine learning outlook. URL: https://oxfordre.com/environmentalscience/view/10.1093/acrefore/9780199389414.001.0001/acrefore-9780199389414-e-572 .doi: . Appendix A. Automated oil-spill simulations
To simulate an oil spill with advection and di ff usion butwithout the need to choose an eddy di ff usion coe ffi cient, the16dvection-di ff usion equation is solved in Lagrangian terms in-cluding an automated method to determine an eddy di ff usioncoe ffi cient. Mathematically, and considering two-dimensionalhorizontal transport, this amounts to solving a stochastic di ff er-ential equation (SDE), given byd X = ( u + ∇ κ ) d t + √ κ d W ( t ) , (A.1)where κ is assumed to be a scalar function of space and time,and where the random variable W ( t ) is a two-dimensional Wienerprocess (see, e.g., Kloeden & Platen, 1992, p. 28, 70). Here, wehave assumed that the di ff usivity is isotropic, i.e., it is the samein both horizontal directions. For details of the anisotropic case,the interested reader is referred to Spivakovskaya et al. (2007).If we solve this equation for a su ffi ciently large number of par-ticles, the density of particles will evolve in time in the sameway as the concentration, C , described by Eq. (3) (Lynch et al.,2014, p. 122–126).The di ff usion part is typically modeled as a random walk,by numerically solving Eq. (A.1), with u = u =
0) is X n + = X n + ( ∇ · κ ) ∆ t + √ κ ∆ W n . (A.2)Here, X n is the position of a particle at time t n , ∆ t is the timestep, and ∆ W n are the increments of the two-dimensional Wienerprocess. That is, ∆ W n is a vector with two independent, iden-tically distributed random components, with zero mean, andvariance ∆ t . If the di ff usivity is spatially variable, accountingfor its gradient in Eq. (A.2) avoids nonphysical transport awayfrom regions of high di ff usivity (Lynch et al., 2014, p. 125).However, this problem is usually more important for verticaltransport, as di ff usivity gradients are usually both sharper andmore persistent in the vertical (Nordam et al., 2020).The di ff usivity can be estimated by di ff erent means, and issometimes provided by an ocean model, but it will usually in-clude uncertainty and errors. It is also important to remem-ber that the eddy di ff usivity does not directly correspond to anyphysical, measurable quantity in nature. Rather, it is a param-eterization of the combined e ff ect of unresolved eddy motion(subgrid stirring), and molecular di ff usivity. Note that since theeddy di ff usivity is intended to compensate for unresolved fea-tures in the ocean model, the eddy di ff usivity will be higher forlow-resolution models, and smaller for high-resolution mod-els. A simple scheme suggested by Smagorinsky (1963) scalesthe eddy di ff usivity with the square of the model grid cell size,which may be a useful rule-of-thumb.Another option is to use a time-dependent di ff usivity, whichincreases with the time since the release. The rationale for thisapproach is found in observations. In the ocean, eddies exist at awide range of spatial scales, from the largest basin-scale gyres,down to the Kolmogorov length scale of millimeters or less.The e ff ect of these eddies on a patch of dissolved tracer dependson the size of the eddy, relative to the size of the patch. Eddiesthat are much larger than the patch will only advect it, with littleor no change to its shape. Eddies that are much smaller than the patch will only serve to deform its surface, without changing itsoverall shape. Eddies that are of the same size as the patch, onthe other hand, will significantly change its shape, by stretchingout filaments in di ff erent directions, thus increasing the overallsize of the patch.A small patch will initially be most a ff ected by small ed-dies, but as it grows in size, it will be a ff ected by increasinglylarge eddies. By constructing an argument based on the typicalturnover time of eddies of di ff erent sizes, it is possible to arriveat an expression for how fast the size of the patch will grow intime. Following the argument of (Davidson, 2015, pp. 257–258), we let R be the mean radius of an initially small andspherical patch, and let the typical speed of an eddy of size r be v r ∼ ( ε r ) / , where ε is the turbulent kinetic energy dissipa-tion rate per unit mass. Since the patch is mainly a ff ected byeddies of its own size, we getd R d t ∼ v R ∼ ( ε R ) / . (A.3)Rewriting this expression by using that dd t R = R d R d t , we getd R d t ∼ ε / R / , (A.4)which is known as Richardson’s four-thirds law Richardson (1926).This expression is only valid for η (cid:28) R (cid:28) (cid:96) , where η is thescale of the smallest eddies (the Kolmogorov scale), and (cid:96) isthe scale of the largest eddies (Davidson, 2015, p. 258). A fur-ther limitation in our case is that on large scales, the ocean isessentially two-dimensional. We will return to this point.Since R is proportional to the variance of a patch of tracer,we see that the rate of increase of the variance is size-dependent,and thus time-dependent, when a patch is subject to turbulentmixing. This is contrary to the case in Fickian di ff usion de-scribed by Eq. (2), where the variance grows linearly with time,proportional to the di ff usivity:d R d t ∼ K . (A.5)From the above, we can derive a time-dependent “e ff ective dif-fusivity”, K e ff ( t ), for a patch subject to turbulent mixing. Inte-grating Eq. (A.3), we find that r ∼ ε / t / , and inserting thisinto Eq. (A.4) we find K e ff ( t ) ∼ ε t . (A.6)Hence, we find that the variance of a patch subject to turbulentmixing is proportional to t , since it grows at a rate proportionalto t .Early experimental investigation of the above results includeobservations of balloons released into the atmosphere Richard-son (1926), and bits of parsnip thrown into a loch by Richard-son’s cabin in Scotland (Richardson & Stommel, 1948). Okubo(1971) published a summary of several dye release experiments,covering spatial scales from 100 m to 10 km and time scalesfrom hours to several weeks. When plotting variance as a func-tion of time (Fig. 1 in Okubo (1971)), he found R ∼ t . , and17hen plotting e ff ective di ff usivity as a function of spatial scale,he found K e ff ∼ R . . These results were later expanded withmore observational data, by, e.g. , Murthy (1976) and Lawrenceet al. (1995), still showing approximately the same trends.If we for the moment accept sloppy notation with respectto units, an explicit expression for the time-dependent appar-ent horizontal di ff usivity, K a , may be obtained from Eq. (3) inOkubo (1971), σ rc = . · t . , (A.7)where σ rc is measured in cm and t is in seconds. Combiningthis with the relation K a = σ rc / t , we get K a = . · t . , (A.8)where K a is given in units cm / s . These observation-basedresults do not agree with the theoretical considerations sum-marised in Eq. (A.4). However, it is clear that the ocean cannotbe considered to be three-dimensional when considering a patchof size 100 m or above, released in the mixed layer. Hence, thetheoretical results cannot be expected to hold exactly.Based on the discussion above, it might seem reasonable touse a time-dependent di ff usivity in an oil spill model. However,it is important to remember that the e ff ective di ff usivity is in-tended to mimic the mixing due to eddies in the ocean currents.If high-resolution current data is used as input to the modeling,more of those eddies will already be represented in the data,and need not be accounted for in the time-dependent di ff usiv-ity. Hence, the di ff usivity should in some sense be matched tothe resolution of the ocean current data.If the horizontal resolution of the current data is ∆ x , thenany patch of tracer with R (cid:28) ∆ x will only be advected alongthe currents, without changing its shape significantly. Hence,it makes sense to apply a time-dependent di ff usivity to smallpatches. However, once the patch grows in size such that R >∆ x , di ff erential advection by eddies represented in the currentdata will lead the patch to spread out further. Applying an addi-tional time-dependent di ff usivity to such a patch will then leadto too much di ff usion.In practice, it is easier to truncate the e ff ective di ff usivitybased on time, rather than spatial scales. It is a simple mat-ter to keep track of the “age” of numerical particles, and usea time-dependent di ff usivity in the random walk for each par-ticle. Future work will be need to determine when to truncatethe time-dependent di ff usivity, and to quantify the di ff erence oftime-dependent di ffff