UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simulating and visualizing marine oil spills van Blankenstein, David 1992

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1993_spring_van_blankenstein_david.pdf [ 7.83MB ]
JSON: 831-1.0051275.json
JSON-LD: 831-1.0051275-ld.json
RDF/XML (Pretty): 831-1.0051275-rdf.xml
RDF/JSON: 831-1.0051275-rdf.json
Turtle: 831-1.0051275-turtle.txt
N-Triples: 831-1.0051275-rdf-ntriples.txt
Original Record: 831-1.0051275-source.json
Full Text

Full Text

SIMULATING AND VISUALIZINGMARINE OIL SPILLSByDavid van BlankensteinB. Sc. (Computing Science) Simon Fraser UniversityA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTERS OF SCIENCEinTHE FACULTY OF GRADUATE STUDIESCOMPUTER SCIENCEWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1992© David van Blankenstein, 1992In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.Department of  Co tr, F si+ t se, 6 G I e The University of British ColumbiaVancouver, CanadaDate^7;') -- I.^Q2.-DE-6 (2/88)AbstractThe main goals of this thesis are to develop a system to simulate and visualize thedynamics of a marine oil spill. This necessitates the development of computational modelsof the properties of petroleum, the containing environment, and a melding of them intoa consistent simulation of the fate and transport of oil spills. Techniques of computergraphics will be applied to present the state of this system as a visual output. Througha large set of accessible parameters and data, the variation of the system in terms of spillsimulation, and model tuning, will be visualized over both time and space.The resulting system is a combination of research from a spectrum of scientific andtechnological areas. It provides a means of abstracting the problem of pollution throughan accessible and consistent collection of mathematical models, and presents its resultsfor interpretation through visualization. As a practical illustration of the system, it willbe tested and verified in the region of southwestern British Columbia.iiTable of ContentsAbstract^ iiTable of Contents^ iiiList of Tables viiiList of Figures^ ixAcknowledgement xi1 Introduction 11.1 Oil Spill Pollution ^ 11.2 Visualization and Simulation ^ 21.3 This Thesis ^ 31.3.1^The Database ^ 41.3.2^The Interactive Visualizer ^ 41.3.3^The Extended Visualizer 41.3.4^Testing the System ^ 52 Oil Spills 62.1 Oil Spills ^ 62.2 The Properties of Petroleum ^ 72.3 The Physics of an Oil Spill 92.3.1^Spreading ^ 9iii32.3.2^Drift or Transport ^2.3.3^Evaporation 2.3.4^Vertical Dispersion, Dissolution and Emulsions ^2.3.5^Grounding and Sinking ^2.3.6^Other Processes 2.4^Spatial and Time Scales of an Oil Spill ^2.5^Visual Characteristics of Petroleum on Water The Environment3.1^The Environment ^3.1.1^Study Area 3.2^Current Inducing Processes ^3.2.1^The Tides ^3.2.2^The Winds 3.2.3^Other Oceanographic Flows ^3.2.4^Waves and the Sea State 3.3^Processes Not Considered ^121313141516161919192021232525264 Visualization 294.1 The Role of Visualization ^ 294.2 Current Trends in Visualization 314.3 Interactive Simulation and Visualization ^ 344.4 Visualization of Model Problem ^ 355 Previous Work 375.1 Statistical Modelling ^ 375.2 Transport Fate Models 39iv5.3 Visualization of Previous Systems ^  426 Modelling^ 446.1 Modelling  446.1.1 Restriction to the Ocean's Surface^  446.2 Modelling the Environment ^  456.2.1 The Terrain ^  456.2.2 The Land Query  466.2.3 Grid Data Description ^  486.2.4 The Tide ^  496.2.5 The Wind  526.2.6 Empirical Wind Modelling ^  536.2.7 Numerical Wind Modelling  546.2.8 Interpolation of Wind Field Data ^  556.2.9 Interpolation of Wind and Tide  566.2.10 Other Advecting Flows ^  586.3 Other Environmental Factors  586.4 Modelling of the Oil Spill Process ^  596.4.1 Oil Release ^  616.4.2 Transportation  616.4.3 Turbulent Diffusion and Spreading ^  626.4.4 Evaporation and Dissolution ^  666.4.5 Grounding ^  676.4.6 Weathering Processes not Modelled ^  686.5 Determining the Time Step ^  697 The System 717.1 Interfaces ^ 717.2 The Database 727.2.1^Data Storage and Extraction^ 737.3 The Interactive Visualizer ^ 747.3.1^Running a Simulation 757.3.2^Displaying the Simulation ^ 777.3.3^Sampling from a Simulation 797.3.4^Flow Visualization ^ 807.4 The Extended Visualizer 837.4.1^Determining the Thickness Map ^ 847.4.2^Interpolation of the Slick Between Frames ^ 867.4.3^Other Rendering Types^ 897.4.4^Texture Mapping an Oil Spill ^ 928 Experiments and Examples 948.1 Standard Slick Formulation ^ 948.2 Testing for Sensitivity to Space and Time^ 948.3 Testing for Spreading and Advection 978.4 Testing the Oil Source ^ 1008.5 Interpolation Between Frames ^ 1018.6 Altering Parameters Based on Tests 1028.7 Texture Mapping a Slick into the 3-Dimensional Model ^ 1039 Conclusion 1089.1 Future Work^ 1089.2 Conclusion 110viAppendices^ 111A Appendix 112A.0.1 Determining the Wind Drift ^  112A.0.2 System Specifics^  114Bibliography^ 117viiList of Tables2.1 Various Oils and Their Properties^ 93.1 Semi-Diurnal Components of the Tide  228.1 Default Slick Model ^  94viiiList of Figures2.1 Schematic of Processes at Work on Oil Slick ^ 102.2 Time Scale of the Various Processes at Work on a Large Oil Spill ^. . . 173.1 Tidal Forces Resulting from the Moon ^ 203.2 Vancouver and Surrounds ^ 286.1 Bilinear Interpolating of Tide Flow ^ 576.2 Spreading of Model Slick as Concentric Ellipses ^ 656.3 Particle-Shoreline Avoidance ^ 687.1 A Model Slick as Viewed in the Visualizer ^ 787.2 Tidal Flow, Time 0900:17-Jul-92 ^ 817.3 Drift Method, Time 0900:17-Jul-92 817.4 Wind Flow, Time 0900:17-Jul-92 ^ 827.5 Viewpoint and Sampling Region 867.6 Sampling of p x q Grid ^ 867.7 Interpolation of an Oil Parcel Over Time ^ 907.8 Current Flow in Howe Sound ^ 927.9 Current Flow in Vancouver Harbour ^ 928.1 The Standard Oil Spill ^ 958.2 Position 200m West 968.3 Position 200m East ^ 968.4 Time Set 15m Back 97ix8.5 Time Set 15m Ahead ^  978.6 No Turbulent Diffusion  988.7 No Physical Spreading ^  988.8 7h = 1.0m2s-1 ^  998.9 -yh = 25.0m 2 .9-1 ^  998.10 No Wind or Tide, Both Spreading Types On ^  1008.11 a =0.08 ^  1018.12 at = 0.10  1018.13 No Effects from Spreading or Diffusion ^  1028.14 Adaptive Time Step ^  1038.15 A i = 300 Seconds  1038.16 Density po = 0.75 ^  1048.17 Density po = 0.99  1048.18 Slick Thickness Map ^  1058.19 Resulting Edges  1058.20 Results of Interpolation ^  1058.21 Slick Produced Using Adaptive Spread Algorithm ^ 1068.22 An Oil Slick in Vancouver Harbour ^  107AcknowledgementThe work presented here would not have been possible without the help from a numberof people scattered across four Departments and two Universities.I would like to thank my supervisor Alain Fournier, for his frank and helpful discus-sions and guidance. Thanks too, to my second reader Brian Klinkenberg of the Depart-ment of Geography. Appreciation is extended to Guido Marionne of the Department ofOceanography for his help in acquiring the tidal model and for assisting me in manymatters revolving around oceanography. Thanks also to Douw G. Steyn in the Depart-ment of Atmospheric Sciences, and to Keith Ayotte of York University for providingmeteorological advice and the wind model.In the Department of Computer Science I would like to thank the many characterswho lurk about the graphics labs, in particular John Buchanan, Peter Cahoon, ChrisHealey, Pierre Poulin, Chris Romanzin, Rob Scharein, and Marcus Tessmann. This goesdouble for Scott Ralph, and for all his years of help and friendship.And finally, a special thank you to my patient and loving wife, Marlene for all hersupport and encouragement over my long years of schooling.xiChapter 1Introduction1.1 Oil Spill PollutionAll industrialized societies use vast quantities of petroleum for fuel and materials. Asthe point of origin and the point of consumption are rarely the same, petroleum mustbe transported both for refining, and for distribution to consumers. This thesis willbe concerned with accidents that occur during the transport and use of petroleum in amarine environment, and in particular with the pollution which results, termed a surfaceoil spill. The frequency of movement, and the great distances covered by shipping, resultsin the inevitable occurrence of marine accidents. Further, land-based human activitiesalong coastal waterways also contribute to the chronic problem of oil spill pollution. Thisthesis will model these seaborne events.Massive spills such as with the Exxon Valdez, where fifty million liters of PrudhoeBay crude was spilled into Prince William Sound, are a major stress on the regionalecosystems, over both time and space. Seemingly innocuous events, such as a few litersof diesel fuel spilling over a vessel's side are far less devastating, of short duration, and yetare far more common. In general, the frequency of oil spill events are inversely related tothe volume of the spill. In the end both contribute to the long and short term degradationof the local marine environment.An oil spill is a complex event over spatial and temporal ranges. The circumstance ofthe spill (in terms of its proximity to the local coastline, size and composition), and the1Chapter 1. Introduction^ 2state of the marine environment (described by the winds, wave, tide and other naturalprocesses), will dictate where the slick will transport, spread and eventually disperse intothe water column The consistent modelling of this transportation and fate was a majorgoal of this thesis.1.2 Visualization and SimulationComputer graphics has fostered the growth of visualization as a mainstream method ofinterpreting the results of experimentation. With the advent of the digital computer atool was provided which could both analyze, transform and display data in large quan-tities. Over the years, as processors and display devices have increased in power andnumbers, and the techniques and algorithms of computer graphics have improved andmultiplied, analysis of data by visual means has become more and more attractive. Thisthesis draws from this body of knowledge, and models, simulates, and visually presentsthe problem of oil spill pollution.The system is designed for consistency where possible, but due to the intrinsic com-plexity of the subject area many broad assumptions are made. The domain is restrictedto be a small subset of the earth's surface, where approximations are made of the majoroceanographic, meteorological and chemical processes which directly effect the petroleumslick's transport and degradation. Within this modelled region, and in context to theassumptions made, estimates of the fate of spills in this region, and their potential effecton the marine biota, can be realized.The methods involve a hybrid combination of both numerical, analytical and empiricalmodels, and empirical and estimated data. From the resulting scheme, an inferencetowards a course of action can be made from interpretation of the potential damagecaused in a simulated environment by a given case of an oil spill. This is judged primarilyChapter 1. Introduction^ 3by visual means, through representations of the oil spill and its interactions with itsenvironment. The primary flavours of visualization are in the animation of the transport,degradation and breakup of the spill, representations of the flows and processes thatform the marine environment, and the rendering of affected areas. These representationsrange from simplistic images of component models, to realistic animations of specific spillscenario. This allows not only for the inference of the potential damages as the result ofa given oil spill, but can also present the problem in a setting which is familiar to peoplein general. To view an animation of an oil spill spoiling a local beach can bring homethe possible result of petroleum use and misuse, without actually causing damage to theenvironment.Variation of the many parameters of the component models, and visualization of theresults have been facilitated. Viewing the model and the spill event not only over timebut over the variation of certain parameters provides information on both the oil spillprocess, and on the merit of the component models. This allows the user of the systemto adjust certain parameters of the component models for their specific needs.1.3 This ThesisThis thesis investigates the problem of oil spill pollution, and how it interacts and isaffected by the environment it is contained in, by providing a visual representation forthe user's interpretation. The problem and properties of marine oil spills will be discussed.A review of previous work in modelling and simulating oil spills, and in the current trendsin visualization of scientific and industrial phenomenon will be presented.Computational models of a subset of the driving meteorological, and oceanographicforces and processes, and the transport and dispersion of oil spill pollution are developed.These models and paradigms are integrated into three distinct environments.Chapter 1. Introduction^ 41.3.1 The DatabaseA database of state data for the various models is developed. A graphical front end willallow the selection of subsets of the spatial domain, and the portions of the various modelswhich apply to it. The ability to introduce data from other sources into this database,and extensions to other models are discussed. This primary duty of the database is toprovide consistent and reliable data sets for testing and composing oil spill scenarios.1.3.2 The Interactive VisualizerInteractive simulation and visualization environment are constructed. Their task is toallow the user to initialize, execute and visualize oil spill scenarios. It also allows theuser to save the state of these simulations over time, visualize the major components ofthe environmental models, and allow the straightforward variation of particular modelparameters. Flow visualization, and the slick distribution over time, will be the primaryfocus of the visual portion of this system. Sampling parameters, selection of variousspatial subsets, and altering of graphics settings is also supported.1.3.3 The Extended VisualizerThe final environment allows for extended sampling of the system and the post processingof simulation data. Some forms of visualization are too computationally intensive to in-clude in an interactive system. The visualization tasks allocated to the interactive systeminvolve the display of model output readily extracted or generated from the simulationin action. Transformation of simulation output into more complex visual representationsis the task of the extended visualizer. Other duties of this subsystem include the con-version of simulation data to other (graphic, visualization and animation) systems. Thesynchronization of time sampled output sets is formalized in a simple script language,Chapter 1. Introduction^ 5which also provides the mechanism for generating output from this subsystem.1.3.4 Testing the SystemA series of scenarios are generated from this system. These are intended to verify themodels, and to allow for a broad range of views of the outputs. These visualization andgraphic images range from the crude to the complex, and will cover not only oil spillsimulations, but also the major environmental models.The thesis concludes with a discussion of the oil spill problem and the knowledgegained from modelling and visualization, and directions for this type of multidisciplinaryapplication of computer graphics and computational modelling.Chapter 2Oil Spills2.1 Oil SpillsMarine oil spills are the intentional or accidental venting of some volume of petroleuminto a coastal or oceanic system. It is a chronic problem which plagues all nations whichrely heavily on petroleum, and as a result effects essentially all mankind.The specific actions that result in oil spills are many and varied. In many cases it isa result of slight and perennial seepage of oil into the environment. Offshore drilling, on-shore storage, pipelines, petroleum refineries and areas of intense marine activity usuallycause this problem to some degree. This form of petroleum pollution usually manifestsitself locally about the source or sources in the form of long term changes in bio-diversity,as some marine and tidal organisms perish under the prolonged effect of petroleum con-tamination [Pe186, Tip82].The more abrupt and certainly more recognizable form of oil pollution is the spilling ofa raw or refined petroleum in its liquid state from a vessel or processing station. Spills oflarge proportion are usually the result of a major accident, such as a grounding or sinking,or from inopportune neglect. The source may be fixed, as in the case of a drilling platformor grounded vessel, where the oil escapes from a rupture over time. In other cases, thevessel may drift with the currents, or be moving under its own power. On March 24,1989 the fully laden super tanker Exxon Valdez tore its bulkheads open on the rocks ofPrince William Sound in Alaska. Approximately 50,00,000 liters of Prudhoe Bay crude oil6Chapter 2. Oil Spills^ 7spilled into the sound, spoiling extensive stretches of the Alaskan coastline. Ecosystems,primarily the intertidal zone and various estuarine regions, were heavily contaminated bythe slick, which travelled over 500 kilometers through the western Alaskan islands beforefinally breaking up in the Pacific [Ven90]. The effect of this spill is still being felt today,as oil buried in shoreline sediments leaches. Other massive spills, such as the AmocoCadiz spill (which occurred of the coast of northern France in March 16, 1978 spilling180,000,000 liters of crude oil), resulted in similar short and long term environmentaldevastation of a large coastal region. The environmental effect of these spills is thedisruption of the inter-tidal and estuarine biota. Toxic level contamination of the watercolumn, and the total smothering of the intertidal zone are the most drastic environmentalshock realized. Further, many human activities such as commercial fisheries, and tourismsuffer as a result of these spills [Tip82]. These large spills tend to receive a large amountof attention, but are relatively infrequent. This statement is meant in no way to belittlethe magnitude of the damage caused by large spills, but only that they are not theentire problem. A far greater volume of oil finds its way into the marine water systemfrom smaller, more frequent spills [Tip82, Mac78]. Vessels flushing their bilges, carelessrefueling, and other small scale accidents introduce levels of petroleum contamination ona regular basis in heavily traveled coastal areas. While seemingly insignificant relativeto their larger counterparts, these spills still cause local damage to the environment.Further, if they are frequent enough, long term degradation of a shoreline will occur.2.2 The Properties of PetroleumThe term oil refers to a variety of organic and inorganic liquids. Petroleum, which islatin for rock oil, is distinct from oil in general in that it is formed in rock formationsover millions of years. For the remainder of this thesis the terms oil and petroleum willChapter 2. Oil Spills^ 8be interchanged, but the reference is always to petroleum.Petroleum is usually found in the form of a dark coloured, highly viscous liquid.However, it can also be found in solid, gaseous and composite forms (oil sands andshales). It is extracted and refined into various fuels, lubricants and petrochemicals.Petroleum is a colloidal dispersion of gas, liquid and solids, formed primarily of shortand long chain hydrocarbon molecules, with impurities in the form of sulphur, oxygen,nitrogen and miscellaneous metals [Jor80, Pur57]. The fraction of the various long andshort chain hydrocarbons largely determines the physical properties of a given petroleumcompound. Short chain molecules tend to be light, with low boiling points, while thelong chain molecules tend to be denser, with higher boiling points.Some important measures assigned to petroleum are the fractions of the various hydro-carbons and impurities, its density, viscosity and specific gravity. A given petroleum typeor product can be described by providing the percentage or fraction that each distinctcompound forms. This is usually stated as a mole fraction' The density and viscosity 2describe the general property of the petroleum. The density is usually reported as aspecific gravity3. Low viscosity oils, such as gasoline and kerosene have a viscosity lessthan 30 centistokes, with crudes and light bunker oils in the range of 30-2000, and highviscosity petroleum, such as heavy bunkers and weathered crudes, having values greaterthan 2000 centistokes. Table 2.1 lists the properties of petroleum crudes and products[Enc87].1 The mole fraction of a component of a substance is the ratio of the number of molecules of thecomponent per molecules of the substance.2Viscosity of a solid or liquid is internal resistance to spread.3Specific gravity is the mass per unit volume of a compound over the mass per unit volume of water.Chapter 2. Oil Spills^ 9PetroleumTypeSpecificGravityViscosity(Centistokes)Arabian Heavy Crude 0.886 35.8Arabian Light Crude 0.831 5.65Iranian Heavy Crude 0.966 840Kerosene 0.787 4.74Light Residual 0.924 146Heavy Residual 0.984 5400Table 2.1: Various Oils and Their Properties2.3 The Physics of an Oil SpillWhen petroleum is spilled into a large body of water, a complex series of competing pro-cess begin to act immediately. Occurring at both large and small scales, these processeslead to the eventual dispersal of the oil into the water body, atmosphere and the incidentshorelines. They are driven by the composition of the oil, and the state of the marineenvironment, as illustrated in figure SpreadingThe most important process which an oil spill undergoes in its early stages is dispersionunder the influence of turbulent and physical spreading. This is the lateral dilation ofthe spill under the influence of number of competing forces. This is one of the mostimportant process which an oil spill undergoes, as it will largely determine the area ofcontamination, and the rate of the various weathering processes.Initially oil spreads under the force of gravity. This force is the dominant spreadingforce for up to the first week of an oil spill's life [Pay85, Jor80, Hua83]. This force isrelated to the thickness of the oil film, and the difference in density of the oil and thewater. Gravitational spreading is counterbalanced in the first hours by inertia, whichis rapidly overwhelmed by friction. The frictional term becomes the dominant force inChapter 2. Oil SpillsFigure 2.1: Schematic of Processes at Work on Oil Slickrestricting the oils spread, and is in effect from the first hours up to a week after spilling[Jor80, Spa88]. Friction increases with surface area, and, as a result, with diminishingfilm thickness. The rate of this spreading is affected by the thickness gradient acrossthe slick body. Differential thickness results in different spreading rates across the oil.Spreading under the influence of surface tension dominates during the later stages ofan oil spill. This tends to occur only with slicks that last more than seven days. Thedifferences between air/water and air/oil surface tensions are the prime cause of thisforce, which is independent of slick thickness. It tends to affect only those oils which10Chapter 2. Oil Spills^ 11have low density and a high concentration of aromatic compounds.If we consider the special case of a calm sea with no wind or current, then the spreadingis uniform and radial about the source. Dependent on the density and viscosity of theoil, a minimum thickness will be reached, and spreading will stop. However, this is anextremely ideal case which would never occur outside of a laboratory testing tank. Inthe presence of wind and surface currents the body of the slick will distort. The velocitygradients of the currents will result in non-uniform spreading. In the case of point sources,the body tends to elongate into an elliptic shape. The long axis of the ellipse tends togrow at a far greater rate than the short axis, and is oriented in the direction of thesurface winds [A1R89]. Further, most of the volume tends to accumulate near the leadingedge (with respect to the local wind currents) of the slick. This process is called sheardiffusion [E1186, Spa88], and is strongly influenced at this scale by high frequency, shortterm wind patterns and turbulent surface currents.The oil's viscosity and density, and the ambient and water temperatures also areinvolved in determining the rate of the spreading. A decrease in viscosity and den-sity results in a decrease in spreading. An increase in temperature, either from theatmosphere or the ocean, results in increased spreading. The reciprocal also holds true.These effects are not observed to be of much significance when the spill involves heavierpetroleum compounds, such as unrefined crudes and residual oils, as their constituentsare less susceptible due to their higher boiling points. However, the effects are readilyobserved in lighter oils and refined fuels, as their constituent are affected by temperaturechanges present in the environment, and undergo greater changes in composition overtime [Mac77, Hua83, Spa88].Chapter 2. Oil Spills^ 122.3.2 Drift or TransportThe drift is the large scale movement of the center of mass of the slick body. It isprimarily driven by large scale forcing from wind, wave trains and surface currents. Thisprocess occurs independent of slick volume and the rate and extent of spreading.The largest influence on the trajectory of an oil spill is from low frequency, longperiod wind currents. Wind shear stresses at the air/water interface induce a surfacecurrent which will cause the slick body to drift. The magnitude of this drift has beenestimated in relation to the magnitude of the wind, and has been found to be in therange of one to five percent of the speed of the wind. Further, a drift angle of up to 45degrees counterclockwise of the wind direction' is associated with this current (This isexpanded in the Appendix). These drift currents extend into the water column, resultingin subsurface currents.The transport of the slick body is influenced by local tidal, geostrophic and waveinduced currents. In the open ocean the effect of tidal currents is of small consequence,as the oscillatory nature of these surface currents tends to lead to a small net trans-port. However in regions where shallows, islands, reefs and other forms of topographycause blocking and deflection of these currents, the effect on the transport is far morepronounced [Bow83, Cre88].The net advecting current is the superposition of these flow fields, although someresearch suggests that this simple linear sum may not be correct, due to non-lineareffects in the system [A1h75, E1186]. The transport of the slick by this current is theprimary cause of shoreline contamination, one of the most serious consequences of theslick's entry into the marine environment.41n the northern hemisphere.Chapter 2. Oil Spills^ 132.3.3 EvaporationOil is composed of various hydrocarbon compounds, which when considered as separatecomponents, have different properties such as boiling points, specific gravity, vapourpressure and densities. The differences in composition of these components results indifferential evaporation rates throughout the petroleum body. This thesis is primarilyconcerned with the lighter fractions of an oil body. Evaporation from light crudes andrefined fuels' can result in the removal of as much as 75 percent of the slick mass. Withheavier crudes, and residual compounds, the loss of mass can be as low as five percent.The majority of lighter components will be lost in the first 24-48 hours after exposure.Over time, the removal of these materials from the slick body results in increased viscosityand density. Evaporative rates are governed by the wind speed, the area of the slick,incident solar radiation, air and sea temperature and the turbulent state of the sea surface[Mac77, Hua83, Spa88).2.3.4 Vertical Dispersion, Dissolution and EmulsionsDispersion and dissolution occur when oil from the spill enters the water column in adissolved or dispersed form.In the dissolution process the lighter components are lost to the water column inparticle form, which is a result of surface agitation and wind speed. The mass lost tothis process tends to be small, as it is competing with an higher evaporation rate for thesame fraction of the oil body. The amount of loss due to dissolution has been reportedto be two to three orders of magnitude smaller than that due to evaporation [Hua83].Vertical mechanical dispersion is a result of the turbulent nature of the water surface.This mechanical action, in the form of wave and wind action, drives some amount of the'Such as gasoline,diesel or aviation.Chapter 2. Oil Spills^ 14oil body under the surface in the form of oil droplets. The majority of these droplets riseand rejoin the main body, while the smaller particles diffuse downwards into the watercolumn. This process is extremely important, in that it will determine the long termlifetime of the oil spill. On a calm sea, this process is very weak, and the slick body, afterthe removal of the volatile compounds by evaporation, will be denser and more stable,and can possibly sink. However if it does not sink, in the absence of vertical mechanicaldispersion, it can persist, increasing the possibility of the slick grounding. In a roughsea, a large amount of the slick volume can be dispersed into the water column. Thisprocess is a function of oil slick thickness, the interfacial tension between the oil andthe water, the sea state, and the fraction of the sea that is covered by breaking waves[Hua83, Spa88, Mac77].In some cases a spill of oil will spread out to monomolecular thickness, and thendisperse under the influence of evaporation and mechanical dispersion. In other casesoil-in-water and water-in-oil emulsions will form as a result of mechanical actions fromsurface turbulence and wave actions. The resulting compounds do not appear to dispersereadily. It is termed a water-in-oil emulsion if the water content is less than 80 percent.It has been found that when more than 50 percent water is present the slicks spreadingis retarded, and a stable emulsion will form [Nae79]. The emulsion will take much longerto dissipate, and will require an applied force in the form of high energy wave actionsto disperse the slick body. Further, the formation of these emulsions can increase theeffective volume of the spill up to five fold [Hua83].2.3.5 Grounding and SinkingSinking occurs when the oil slick has achieved a specific gravity greater than that of thesurrounding water body. This results partially from weathering processes, but also fromthe binding of sediments to the oil body.Chapter 2. Oil Spills^ 15Grounding is caused by the deposition of the oil spill onto the inter-tidal zone. Therate and amount of oil which will collect on the shoreline is dependent on the ebbing andflow of the tides, the nature of the shoreline, and the properties of the spilled petroleum.Each type of shoreline will have a finite capacity for oil accumulations. The importantfeatures associated with a shoreline are its slope, and its surface characteristics. Flat orlow sloping areas tend tend to allow the oil to collect upon them. Tidal flats, marsh areasand gently sloping beaches have this property. If the grade is is more severe, then the oilwill not as readily adhere to it. The surface properties of the shoreline can be describedin terms of its roughness and the porous nature of the shoreline. Rocky, gravelly, andsandy coastlines allow oil to adhere in differing amounts. Further, the mixing of oiland shoreline sediments will also determine the amount of oil which is deposited on ashoreline. Dependent on these features, and the amount of oil which has accumulated,shoreline oil will be partially removed by the receding tide. Further information on thistopic can be found in [Gun87].2.3.6 Other ProcessesPhotoxidation of hydrocarbons, caused largely by solar radiation, results in the produc-tion of water soluble compounds from a slick body. The physical processes at work arepoorly understood [Hua83, Spa88]. It has been observed to be a process which takesweeks to months to occur, and affects only a small portion of the slick volume.Biodegradation is the slow process which results from marine organisms decomposingthe oil slick. It is also a poorly understood process, with most knowledge derived fromlaboratory observations. It is of little consequence in the first 15-30 days of an oil spill,and becomes marginally more important after this period.Chapter 2. Oil Spills^ 162.4 Spatial and Time Scales of an Oil SpillThe spatial extent of an oil spill and the range of its transport, and the time scale underwhich an oil spill will exist, is clearly determined by the size and composition of the oilsource, and the state of the environment. See figure 2.2. This must be tempered in thecase of a large-scale accidental spill where the oil source is introducing petroleum intothe environment continually. Hence, the time that will elapse from the entry of the oilinto the marine environment, up until the spill body has completely dispersed, is highlyvariable.2.5 Visual Characteristics of Petroleum on WaterThe visual characteristics of oil slicks on water are a function of the oil type, sea state,slick thickness and the state of weathering that the slick body has undergone.In general the radiative properties of oil on water tends to be similar to the propertiesof the underlying water. If we assume a calm sea, a high solar altitude, and no bottomreflectance, then in the visible spectrum the radiance of an oil body is slightly higherthan the surrounding ocean. The radiant properties of clean sea water and oil contrasthighly in ultraviolet and the infrared bands of the spectrum. The reflectance propertiesof oil and water are almost identical in the green and red bands. This of course mustbe tempered by the composition of the water body. Turbid waters which are carryinglarge sediment loads reflect more light than clean water, and this will affect the contrastbetween the oil and water. The solar altitude largely determines the contrast between oiland water. With a low solar altitude oil is all but indistinguishable from the surroundingwater [Viz73] To be able to reliably distinguish between oil and water readily otherproperties must be observed.Oil has a dampening effect on the ocean's surface. Small scale surface waves, orChapter 2. Oil Spills^ 17Figure 2.2: Time Scale of the Various Processes at Work on a Large Oil Spillcapillary waves are removed by the affect of the riding oil body. The range of capillarywaves which will be dampened depends on the thickness of the oil slick. Large oil spillshave been observed to dampen out even wind driven waves. It is this property whichmakes radar the preferred method of remotely sensing an oil spill [ASP83].If an oil slick is at an advanced stage of weathering, and emulsion formation occurs,then the radiative properties will also change. The formation of water-in-oil and oil-in-water emulsions yields a compound which ranges from a muddy brown to an almostwhite appearance. In the event of emulsion formation, the contrasts between the oil andChapter 2. Oil Spills^ 18water body is far more pronounced.If the oil forms sheets which are sufficiently thin, then an irridescent lustre appears.This is the result of the difference in refractive indices between air, oil and water. Anytwo rays of light with the same wavelength emanating from a distant light source willbe in phase. The difference in the refractive indices of the oil, air and water resultsin a phase shift between rays that were parallel on entry. Destructive and constructiveinterference occurs. Destructive interference tends to remove certain wavelengths, whileconstructive interference reinforces others. In studying this effect in a laboratory setting,[Hor72] found that oil on water with a thickness less than 150 nm appeared as a colourlessfilm, but a thickness in the range of 300 - 1000 nm produced distinct iridescent patterns.The patterns dull from the 1000 nm thickness on, with the effect essentially removed inoils of thickness greater than 2000 nm.In this chapter the general properties of a marine oil spill were discussed. Thisincluded the properties, sources, and effects of oil spill pollution, and the affect of theregional environment on the oil spill.Chapter 3The Environment3.1 The EnvironmentThe rates governing the spread, degradation and eventual breakup of the oil spill are aresult of the state of the marine environment, and the resulting proximities of the slick toland and other features. The mechanics and dynamics of the atmosphere and ocean areextremely complex and variable. Those processes which are of concern to the model willbe identified. Further refinement will specify those activities that occur in the domainwithin which the model will be developed. The scale of the intended research will notinclude events of global scale, nor processes that are described at the level of the molecule.The middle, or mesoscale view of the environment will be discussed.3.1.1 Study AreaThe study area or model domain will consist of North America from the 49th to the 50thlatitude, and from the 122nd to 124th longitude (Figure 3.2). This includes the westernportion of Greater Vancouver, the mouth of the Fraser River, Howe Sound, the southeastern section of the Strait of Georgia, and the Northern Gulf islands. This is a regionof complex topography, formed primarily by the uplift of the Coast mountain range, andthe Fraser River valley. It is a strongly tidal estuarine' system.1 The term "Estuary" is given a broader definition than its standard as a tidal region at the mouthof a river; rather, it is a semi-enclosed coastal region with a free connection to an ocean within whichthe sea water is considerably diluted with fresh water19Chapter 3. The Environment^ 20 dMoonEarthFigure 3.1: Tidal Forces Resulting from the MoonHuman activity in the region, centered around Vancouver, a rapidly growing cityof approximately 1.5 million people, results in a high marine traffic density in bothcommercial and pleasure vessels, in addition to a large inshore fishery and forest industry.A large fresh water discharge from the Fraser River results in a highly stratified watercolumn in the Georgia basin. Many other smaller rivers and streams contribute to thefreshwater flow into the system. The region has a semi-diurnal unequal 2 tide. Winds inthe region are dominated by two large atmospheric pressure systems, the North Pacifichigh and the Aleutian low [Cre88]. Winds from the south east occur in the eastern portionof the domain when the low dominates, with winds from the north west occurring in thewestern half. The high pressure system induces westward blowing winds along the westernportion. In general, the low dominates in the winter, and the high in the summer.3.2 Current Inducing ProcessesAn oil body is transported and dispersed in an estuarine region by the combined effect ofsurface-driven wind currents, tidal currents, geostrophic and residual currents, and wind2 1n a semi-diurnal unequal tide, there are two maxima and two minima over a diurnal cycle. Further,there is a marked difference in the heights of the succesive maxima and minimaChapter 3. The Environment^ 21induced waves.The processes at work in the atmosphere and ocean are extremely complex, occurringat planetary and molecular scales. The laws of thermodynamics, fluid dynamics andplanetary geophysics describe these systems. This thesis cannot hope to completely de-scribe these systems, and so the interested reader is directed to the works of [Pon83] for ageneral discussion of dynamic oceanography, and [Ho179] for similar work on atmosphericprocesses.3.2.1 The TidesThe tide is visually perceived as the rise and fall of the elevation of the ocean's surface.However, the primary phenomenon is the horizontal flow of the ocean in response toforces imparted on the earth by the moon and sun. The rise and fall is merely theconsequence of the inward and outward flow of water. This tide generating force wasshown by Newton to arise from the law of gravitation. If a system comprised of the earthand moon is considered, neglecting the rotation of both bodies around the sun, andassuming they are fixed spheres, then the two bodies would move as a result of mutualgravitational attraction. The force imparted by the moon on the earth is F1 at point Pon the surface, and is F2 at the center of the earth C. This is shown in diagram 3.1. Theterm F is the tidal generating force, which is the horizontal component of F3 = F1 - F2.It can determined by the system3 ME aF = –2 —mm ( ci )3 g sin 20where ME and Mn, are the masses of the moon and the earth, a is the radius of theearth, d is the distance between their centers, and 0 is the latitude of the point P. Therelative positions of the moon and sun with respect to the earth vary with the orbitsand rotations of these three bodies, yielding a number of constituents, each describedChapter 3. The Environment^ 22by its period, phase and amplitude. These constituents are classified by their periods,with a period around 24 hours being diurnal, 12 hours a semi-diurnal, and greater than24 hours a long period constituent. Since the motions of the sun and moon are known,the resultant tide-generating potential can be determined by summing the individualharmonic constituents. Some of these constituents are listed in table 3.1.Constituent Symbol Period(hours)RelativeMagnitudePrincipal Lunar M2 12.42 100Principal Solar S2 12.00 47Large Lunar Elliptic N2 12.66 19Luni-Solar K2 11.97 13Table 3.1: Semi-Diurnal Components of the TideThere are two major theories on the ocean's response to tidal forces, the equilibriumtheory and the dynamic theory. In the equilibrium theory, proposed by Newton, it isassumed that the earth is covered by a sea of fixed depth and density. The ocean respondsto the applied tidal force by converging on the applied force, lifting the surface of theocean in one location, and depressing it in others. A new oceanic equilibrium is attainedwhen the hydrostatic pressure forces from the sloping ocean surface balances the geo-potential forces. This theory underestimates the magnitudes of the oceanic tides, andcannot account for phase lags in constituents.The dynamic theory, first derived by Laplace, also assumes a world girdled by a uni-form ocean. However, the tidal forces are seen to generate waves with the same period asthe constituent. The movement of these waves yields surface currents, while the intersec-tion of the waves with coastal features results in the rising and falling of the sea level weperceive. The M2 tidal force generates a forced wave with a period of 12.42 hours, overthe ocean surface of the model planet. In principal, this theory should provide a meansfor calculating the tides at any given point on the ocean's surface. In reality oceanicChapter 3. The Environment^ 23boundaries and subsurface topography complicate the constituent equations. These fea-tures affect both the phase and magnitude of the many forced waves. In constrictedchannels a tidal wave will respond by accelerating, a result of the changing width anddepth it encounters. If the channel narrows and decreases its depth, then the wave willbuild up on itself, increasing its amplitude. This phenomenon is called funneling. Fur-ther, boundary effects can be found in the vicinity of subsurface features, and points ofland, which cause the deflection of gross tidal currents. In bays, tidal resonance can causeabnormal fluctuations in the tidal elevation, which can be markedly different than thenearby open ocean. Further, tidal waves are strongly influenced by friction in shallow andconstricted waterways. These topographic affects on the tides abound in the Vancouverregion.Clearly the horizontal flow of the tide has a large effect on the transport and spreadingof a slick in a marine coastal region. Further, the rise and fall of the sea level willalternately expose and protect the shoreline, allowing a slick to affect the entire intertidalzone.The effect of the tide on petroleum slicks is small on the open ocean. The oscillationsof the tide over the diurnal cycle tend to result in a small net displacement, and themagnitude of the tidal currents tend to be smaller than that in enclosed coastal regions.However, in regions which are formed of islands, channels and bays, the effect on thetidal current is more pronounced, resulting in complex surface currents.For further reading on tidal phenomenon see [Bow83, Cra86, 0br85]3.2.2 The WindsThe planet earth is surrounded by a layer of air called the atmosphere. Kept in place bygravity, it is a low density mixture of gases such as nitrogen, oxygen, argon, water vapourand particulate matter. Although the earth and atmosphere receive as much energy asChapter 3. The Environment^ 24they radiate, this heat balance is only true on a global scale. Locally there are variationsin the exchange of heat, and these differences place the atmosphere into motion. Thisthesis will be concerned with the net processes that occur in the circulation of the loweratmosphere or troposphere, and ultimately in the horizontal surface winds. The primaryinterest is in diagnosing the winds in the study domain.Wind on a planetary scale is generated primarily by the pressure force and the Cori-olis force. The balancing of these forces is termed the geostrophic flow. The pressureforce results when air in regions of high pressure flows towards regions of low pressure.Regions of the atmosphere that differ in pressure arise primarily from local changes inheat exchanges. The Coriolis force is induced by the rotation of the earth about its axis.This results in the counterclockwise rotation of flow in the northern hemisphere. Thegeostrophic flow determines the global circulation patterns, and flows from the west inthe Northern Pacific. Local wind conditions are largely determined by this wind [Ho179].This thesis will be primarily concerned with wind patterns which have a scale in theorder of one to one hundred kilometers, and are of a duration of between 4 to 10 hours.These include convection currents in unstable air, gravity waves, oscillations in unstableair, and wind patterns induced by friction, topography, and heat exchanges between theland and atmosphere, at what is called the planetary boundary layer. These mesoscalewinds are characteristic of the region they influence.The primary wind flow that must be investigated is induced by boundary effectscaused when the geostrophic wind comes into contact with the terrain of a region. In thePacific Northwest, the mountainous terrain greatly modifies the weather patterns in thelower troposphere [Mas85]. Channeling, blocking and deflection induce patterns whichin no way resemble the synoptic patterns which drive them. These patterns are alsoinduced by differential heating and friction between the land and ocean surface [Pie74].In a statically stable atmosphere, the turbulent mixing in the boundary layer is generatedChapter 3. The Environment^ 25primarily by vertical shear at the ground. The turbulent mixing which transfers heat andenergy away, and momentum towards the surface during the diurnal cycle of heating andcooling, is mechanically, not thermally, driven.Further discussion of this topic will be deferred to chapter 6 where modelling of themesoscale terrain induced wind currents will be covered. See [Ste88, Smi90, Eli77a] forfurther readings on the nature and properties of weather systems and in particular windsystems.3.2.3 Other Oceanographic FlowsGeostrophic or oceanic currents are primarily the result of the rotation of the earth aboutits axis, which is manifested in the form of the Coriolis force. This force results in thecounterclockwise flow of large oceanic currents in the Northern Hemisphere. Further,variations caused by heat exchange also contribute to these large scale flows. However,in the Vancouver and Strait of Georgia region these currents have little or no effect, asthe semi-enclosed nature of this marine body provides a topographic buffer, and so thesecurrents will not be considered.3.2.4 Waves and the Sea StateWaves in the marine environment range over a wide spectrum. The concern in this thesisis primarily with surface waves, which are the result of a sustained wind stress over apatch of the ocean's surface. Larger waves, such as internal and planetary waves will notbe considered. Tidal waves are discussed in a previous section of this chapter. Surfacewaves have relatively short periods, ranging from 0.1 to 30 seconds, with wavelengthsfrom one centimeter to hundreds of meters. The vertical accelerations of these wavescan approach that of gravity with both horizontal and vertical accelerations overcominggeostrophic forces. These waves largely determine the sea conditions of a coastal region.Chapter 3. The Environment^ 26Surface waves can be found far from the area of influence which initially generated themin the form of swell. A swell will decay over time as it travels beyond its origins. If sucha wave enters shallow water, then the bottom will effectively rejuvenate it, causing theswell to peak. The swell will eventually disperse, either within the oceanic body, or onthe shoreline as surf.Wave analysis is a complex and fascinating field, but for the purposes of this researchwe are interested only in the general state of the sea surface, and not with individualwaves. The sea state is a rather qualitative term used to describe the wave state of agiven sea surface. Surfaces waves have both horizontal and vertical kinetic, and potentialenergy as a result of vertical displacement. A measure of the density of energy on theocean's surface is given by E = pgH2 /8 , where p is the density of the water, and H isthe wave height 3 and g is the acceleration due to gravity at sea level. This energy willactively work on promoting the degradation processes which are weathering an oil spill.The practical field of wave and sea state prediction will be turned to and expanded inchapter 6.3.3 Processes Not ConsideredThere are many components of the environment not considered. The effect of fresh andsalt water mixing is an important factor in an estuarine setting. At the mouth of a riverthere will be a distinct layering as the two liquids mix. The rise and fall of the tide helpsin determining if the fresh water forms the top or bottom layer. The mixing is primarilyvertical. Its effect on the net horizontal surface currents away from the mouth of theriver is slight, and as such it will not be considered. The effect of spatial and temporal3Wave height, which is measured from peak to trough, is meant in the statistical sense. One methodsof estimating wave height is to determine the significant wave height, or the mean height of the highestone third of the waves. In many cases this statistic is based purely on visual estimatesChapter 3. The Environment^ 27variation of the ambient and ocean temperatures is important to processes such as largescale weather patterns, and the exchange of heat in the atmosphere. At the scale of themodel problem they are not as important.The major processes which effect the movement and degradation of oil spills in amarine environment have been discussed. These processes include the terrain, wind andtidal currents, and are specific to the test region.Chapter 3. The Environment^ 28Figure 3.2: Vancouver and SurroundsChapter 4VisualizationVisualization in the context of computer graphics and science can be loosely defined asthe transformation of data from some initial state into a visual representation via theprocessing of a digital system. A more frank definition is that visualization is merelythe display of data [Tuf83]. From this representation, interpretation of the data anddetermination of the underlying processes at work is possible. The transformation canrepresent both the mathematical conversion of the data from one state into anotherthrough a model, or to the intermediate filtering and composing of the data into a visualform. This chapter will focus on the role of visualization or data display in scientificand industrial applications. A review of previous visualization research in related areaswill be covered. The function of interaction in these systems will be emphasized. Adiscussion of the application of visualization to the model problem of oil spill pollutionwill complete the chapter.4.1 The Role of VisualizationVisualization in its many forms has always be a major component of the scientific process.In its most unpretentious form, the two dimensional graph of a simple relationship oftencan do more to describe the actual workings of the root process to a human than either theraw data, or the observation of the event. However, in the presence of higher dimensions,more complex processes and volumes of data, this time honored visual technique becomeswoefully inadequate.29Chapter 4. Visualization^ 30Digital technology has both advanced and strained the field of scientific visualization.On one hand it has enhanced our ability to transform and visually present data fromboth observation and mathematical modelling. On the other hand it is producing devicesand systems which create or collect ever increasing volumes of data. A method mustbe applied to clearing the backlog of data produced by these systems. Thus scientificvisualization is presented with a number of major tasks that it must fulfil.The first task that visualization encounters is the interpreting of large volumes ofdata. High bandwidth data sources such as satellites and radio telescopes, and medicalimaging systems, to name but a few, produce voluminous quantities of data which must besummarized rapidly. Further, computational models of physical and chemical processesof ever increasing complexity are being developed. These models produce a boundlesssupply of data. The output from these sources must in some way be refined to allow forinterpretation without removing the underlying patterns of the data.Complex and multidimensional data needs to be presented in a form which allowsfor the determination of trends and order within these sets. We exist in a multidimen-sional world, but our visual system is based on a two-dimensional image of this world.Visualization must take these constraints, and present information, extracted from manydimensions, in a manner with which understanding of the processes in action is achieved.The role of visualization in this context is to rapidly communicate information ona given data set in the form of graphic presentations. It must provide a base of func-tionality, which can be either general or application specific, and that can be appliedto transforming a set of potentially multidimensional inputs into a visual output. Thisprovides a window into local and global properties of the data set.Chapter 4. Visualization^ 314.2 Current Trends in VisualizationScientific visualization has branched out over the years to encompass many fields in thecourse of transforming data into visual image. The methods of computer graphics, vision,computer aided design, signal processing, and computer human interactions to name buta few, are applied to producing visualizations of data from a wide variety of applications.Researchers reiterate the problem associated with attempting to present too muchdata to a user [Hib89, Nie89, Be187, Hab90]. When this occurs, the user's ability tointerpret the resulting image decreases and the effort required to generate the imageincreases. Each image should be derived to show no more than the effect of two or threeparameters. If more information is shown, the usefulness of the visualization begins todegrade. While in many cases exotic graphics can be visually appealing, they quite oftendo not provide any assistance in understanding the processes that they were generatedfrom.The application of animations to the visualization process has largely been realizedthrough the rapid advancements of processing and display hardware. The frames ofthese animations are usually but not necessarily dependent on time. In the case of themodel problem, an animation is formed of the spatial movements of the slick body overtime. However, what variable the frames of the animation are based upon will vary fromapplication to application. Researchers who emphasize the use of animation as a tool forvisualization include [Rap89, Bri91, Be187, Bra92, Hib89, Nie89].Many researchers discuss the linking of visualization to what can be termed as acomputational experiment. Rather than deal with static data generated from an externalmodel or observation, these systems have models coupled directly to the system. Thestate of the system is initialized, and the model executed, with the progress and resultspresented visually.Chapter 4. Visualization^ 32A computational experiment begins with a loosely stated problem. For this thesisthe statement is "What would the effect of a oil spill be on the local environment." Thenext step is to generalize and transform this statement into a well-posed mathematicalproblem, from which a reliable and adequate approximation can be developed. Thesimulation is executed and the output of the model is interpreted. The user can thenadjust various quantities, both for the purposes of verifying models and for consideringdifferent setups in these models. The evaluation of these models by visual means greatlyaccelerates the process.An example of a visualization system developed to be run as a virtual lab was thesystem developed by Briscolini and Sanangelo in [Bri91]. A visualization system wasdeveloped to study the properties of two dimensional turbulence and three dimensionalflows. The main form of visualization was the animation of flows over time, with thevarious physical parameters being altered over each set. The primary value measuredby visual means was the vorticity. Rather than render lines of isovorticity, the effortwas placed in determining the isovorticity of the surface. The blue and red componentsdenoted a clockwise and counterclockwise flow respectively, with areas of low vorticityenhanced by the inclusion of white and green respectively. When extending their workto three dimensional flows, they reported no conclusive method could be found for visu-alizing these flows. Three dimensional vector representations of flows are discussed, andthe authors conclude that this method was relatively ineffective if the flow is not smooth.Rather, the flow was rendered using a colour-table approach, where blue and red areused to represent low and high velocity flow. Intermediate flows grade from the blue intogreen and then red. They concluded that the two dimensional animations allowed forthe observation and understanding of the mechanisms that control turbulent flows.Hibbard and Santek developed a visualization system specific to the data and prob-lems associated with the earth sciences [Hib89]. The range of data sources includedChapter 4. Visualization^ 33meteorological data, image data, and the results of numerical simulations. The emphasisof the system was on the processing of extremely large sets of data, some containing inthe order of one billion observations. The numerical models, which were a prominentsource of data for the system, were far too complex for direct coupling to be considered.Rather the simulations were executed, and the resulting output presented to the system.As the sources were many and varied, an emphasis was placed in providing reliablemechanisms for bringing data into the system. The data, once in the system, could beconverted into various data structures could be used in the visualizations. These struc-tures included grids, images, paths, and structures of irregular data. Tools to interpolatethese structures over time and space were included. The resulting system provided a gen-eral arena for the visualization of meteorological data (both observational and numeric),tectonic and oceanographic data.The spatial data was placed in a setting which composed of a box with height andhorizontal ticks labelled in the appropriate units for the given visualization. The bound-ary of the region in which the data was located is displayed on the bottom of the framebox. Three dimensional scalar fields were displayed using transparent and opaque con-tour surfaces on such scalar values as water vapour and potential temperature. Currentflows were visualized using the trajectories of particles. The thickness and length of theline implies the speed and duration of the flow. Slices of three dimensional data, texturemapping of generated surfaces and other visual techniques were combined to providethe ability to readily construct animations of the given data. This was combined in theanimation with a small rocking of the scene, to provide additional depth cues.Hibbard and Santek conclude by emphasizing the need for interactive systems. Theusers of the scientific community need both the versatility and speed afforded by ananimation based system, but also need to be able to get up and running with the systemwith a minimum amount of training.Chapter 4. Visualization^ 344.3 Interactive Simulation and VisualizationVisual interactive simulation, a term used by Bell and O'Keefe, describes a designmethodology in which the prime motivation is the user's interaction with the model,as opposed to just the portrayal of the outputs [Be187]. This was discussed in the con-text of discreet systems simulations, but the basic philosophy was just as valid for acontinuous simulation.The notion was to afford the user the ability to interrupt and directly manipulatethe model. This was done graphically though interactions with buttons, icons and othervisual input mechanisms. Controlling parameters can be altered, objects removed fromthe simulation, and the system restarted from the point of interruption. The effect ofthese changes on the system will be apparent from the visual display of the model. Thisapproach allows the user to avoid the problem of restarting the simulation each time aparameter was adjusted.These systems use animation extensively to show the state of the model. Icons whichare readily identifiable as real world components of the model were used within theseanimations. Various statistics were presented in different windows, which changed withthe model state. As the problem was shown in a familiar setting, the user was foundto understood the changes occurring in the simulation, and could react by altering theparameters.Bell and O'keefe reported that users found this type of system very appealing. Itprovided a large amount of freedom for the user to shift attention between differentcomponents of the simulation. The end result is that users feel they are participants,rather than passive observers. It was found to be extremely valuable in the verificationand validation stage of a given model. In the long term, the authors suggest a numberof points that can be followed in defining a methodology for such as system. TheyChapter 4. Visualization^ 35include getting the user involved early in the simulation, displaying the visualizations ofthe system as soon as possible, making interactions general, and transferring control ofsimulation directly to the end user.Peskins and associates discuss a system which was developed to allow the user to steerthe development of a model through visual means in [Pes91]. This system is applicationindependent, and was developed to allow a user to quantitatively measure propertiesof the output data, and present those results visually. The ability to adjust modelparameters incrementally, judge the effect through visualizations, and thus steer thedevelopment of the model, is emphasized.4.4 Visualization of Model ProblemThe work of previous researchers gives a gauge for the type of visualization that will beappropriate for the oil spill simulation system. The intended purpose of the system is todetermine the spatial and temporal distribution of an oil slick over a given domain. Torealize this, the models of the system must be tuned for these tasks, and the composingof oil spills by the user easily achieved.The precursor to the visualization must then be to provide an interactive system,which permits the adjustment of the parameters which dictate how the model will performand react.In terms of visualization, the spatial distribution of the oil over the region will bepresented in the form of an image showing the slick's position and extent over time withinthe tested domain. This will provide the user with a condensed visual representation ofthe state of the system for the given setup. In terms of oil spill pollution, this meansthat a researcher can vary position, time, volume, and other variables crucial to theprocess and guide the system towards various goals by interpretation of the system'sChapter 4. Visualization^ 36visual output. If the need is to determine the areas affected by oil released from a givenposition, then altering the release time and volume will produce a series of images fromwhich the effect can be estimated. These images will consist of the position of a thicknessof oil over the test region's sea surface and shoreline for a given set of parameters andtime.As the state of the environment is so critical to the resulting dispersion of the slick,facilities for visualizing the various environmental components are provided.In this chapter the general field of visualization or graphic data display has beendiscussed. The current trends in this field with an emphasis on the physical sciences andhuman-machine interaction are covered. A number of developed systems and visualiza-tion methodologies have also been reviewed.Chapter 5Previous WorkConsiderable research has been directed towards understanding and modelling the activ-ities that cause oil spills, and the events themselves. In reviewing this area of research,a distinction must be made between fate-transports models, and statistical and impactanalysis studies. While this thesis is primarily interested in the former, the latter helpsto form the justification of applying technology to this problem, and gives a sense of theextent of the problem.5.1 Statistical ModellingStatistical models are concerned with predicting the long and short term likelihood ofoil spills occurring in a given region. Using historical data and statistical methodologies,an attempt is made to arrive at distributions which will allow for the estimation ofthe number of accidents, their magnitude, and their potential affect on the surroundingenvironment.Mackay and Wilks provide a comprehensive analysis of the oil spill problem acrossCanada [Mac78]. The problem is approached in a sector by sector analysis of all petroleum-based activities and the resulting incidents of all types of spill pollution over an eightyear period. The intent their model is to determine long term trends, and from this toestimate the distribution of oil spills over long periods of time. The sectors of interestto coastal marine regions are refineries, storage and ocean tanker traffic. Each sector islA fate-transport model tracks the spatial position of the oil slick over time.37Chapter 5. Previous Work^ 38assigned an exposure variable, which are intuitively based, such as port-calls for tankers,or barrel-stored for refineries. Using historic data on actual oil spills, and performingstatistical analysis on these data, spill frequencies and distributions were determined.Conclusions on the spatial distributions of oil spills across Canada cannot be made asthe data is aggregated for all regions. In addition the environmental conditions were notfactored in, as the analysis was based on occurrence, rather than fate and affect. One ofthe conclusions of this analysis was that marine tanker transport is approximately twiceas accident prone as any other sector of the industry.Pelletier developed a ecological risk estimation model [Pe1861. It was based in the St.Lawrence Estuary, but the model is general enough that it could be extended into anydomain. The analysis considered local oceanographic, topographic, and meteorologicaleffects, but only in a very generalized sense. Given the nature and volume of the oil spill,and the environmental conditions, estimates of the probability of shoreline contaminationand the levels of toxic effect on the local biotas can be estimated. The effect of eachfeature, such as the local currents, proximity to shoreline and amount of spreading, wereassigned a sensitivity value. For example, the probability of shoreline contamination wasmeasured in terms of the spill's proximity and an energy description of the shorelinefeatures. A shoreline feature is classified as being high energy if it has a low probabilityof the oil adhering to its surface. This would include areas such as cliffs or stony beaches.Low energy features are areas such as mudflats or marshes. A combination of thesestatistics was then used for estimating open water, shoreline and direct biological impact.It was concluded that environmental risk was dependent upon the spatial location ofthe spill within the domain. This fact largely determined the spreading, and resultanttransport by local currents.Chapter 5. Previous Work^ 395.2 Transport Fate ModelsTransport fate models are primarily concerned with the spatial distribution of the oilbody over time, its physical and chemical state, and the eventual breakup of the oil spill.[Hua83] and [Spa88] give excellent reviews of the subject area, covering in detail boththe issues involved in oil spill and environmental modelling.Ahlstrom in [A1h75] developed a transport-fate model for industry in the early 70's.It is of primary interest to this study because the model domain was the Birch andSemiamhoo Bay region, which is due south of this system's model domain, and that itdeveloped reliable stochastic methods for estimating turbulent dispersion.Major environmental processes, such as the wind and tidal flows, were empiricallymodelled. Processes such as geostrophic surface currents, fresh water runoff and waveinduced currents are discussed, but are not included in the model formulation. A stochas-tic method for modelling turbulent dispersion is designed. The slick body is viewed asa collection of independent particles, which are displaced by a Lagrangian random pro-cess. Further, the particles are independently advected. This method of modelling theslick body and the dispersion has become more or less a modelling standard in this areaof research, used by researches such as [Ven87, A1R89, Lar88]. Physical spreading wasmodelled using the work of [Fay71]. The physical spreading is only considered in the firstfew minutes of the simulation, and the model depends on turbulent dispersion for theduration. Empirical relationships to determine the maximum size of a slick's area andthe time required to attain this scale are presented and justified. The chemical degrada-tion and dispersal of the slick is not modelled. Rather, the effort is applied to accurateforecast of the advection and surface dispersion.Mackay and Leinonen in [Mac77] produced a set of models to study the behavior ofan oil spill on water. These models were coupled with observational data which was usedChapter 5. Previous Work^ 40to verify the model. Adjustment of the model parameters was also facilitated. Empiricalin nature, it provided a step by step description of the analysis of the chemical processeswhich affect a slick upon entry into the water, and parameterizations of these processes.The system was not designed to be run within a domain, under the influences of spatialand temporal variations of the environment. Rather, the system was generalized to modelenvironmental effects from scalar measures of these processes, such as wind speed andsea state. Other model parameters (such as volume and composition) are varied on acase by case basis, in constructing sets of model outputs to validate the assumptions thatthe models were constructed upon.The system modelled the spreading, dissolution, emulsification, evaporation and ver-tical diffusion processes over differing compositions and volumes of oil, wind speed, seastate, and the variation of certain model parameters. The separation of an oil slick intoa thin and thick slick was formalized. This phenomenon has been observed in many oilspill cases. His models for determining the evaporation and dissolution flux from an oilslick appear to have become the accepted standards, as used by [Ven90, Ven87, A1R89].The works of Venkatesh merged both the transport and fate requirements of thisflavour of oil spill modelling. In [Ven87] the model used and developed by the CanadianAtmospheric Environment Service was discussed. It was developed primarily to aid inshort and long term cleanup of spill events. Designed to work on a general domain, thetest region presented was the Bay of Fundy, and the model results were compared againstobservational data from the Argo Merchant spill of 1977. The CEAS model accuratelymodelled the transport of this oil slick. The Exxon Valdez spill's transport was simulatedwith this system in [Ven90]. The modelled drift was faithful to the actual transport of thismajor spill. The system models the advection, turbulent and physical spreading of theslick, and weathering by evaporation, and emulsification. The environmental processeswere derived by accepting inputs from existing hydrodynamic and meteorological systems.Chapter 5. Previous Work^ 41Al-Rabeh and associates in [A1R89] developed a stochastic simulation of oil spillfate and transport in the Arabian Gulf. A complete analysis and derivation of theproblems of advection, dispersion and chemical degradation were presented. Advectionand turbulent were diffusion modelled as [A1h75]. Surface spreading was applied to boththe individual parcels, and to the oil body as a whole. Both the parcel and spill spreadingwas modelled around the elliptic paradigm. A statistically based method of determiningthe wind flow was used [Cek88]. Tidal and geostrophic currents in the test domain, theArabian gulf, were derived from existing hydrodynamic models. The system modelledthe physical spreading, turbulent diffusion, vertical mechanical dispersion, evaporation,and emulsification of test spills. A case study of an number of oil spills in the Arabiangulf was executed to verify the model.Over the past two decades the quality and breadth of the modelling of oil slick pol-lution has certainly improved. However, as most of the chemical degradation, dispersionand spreading algorithms are still highly empirical in nature and the improvements appearto have come from increased resources, and improved hydrodynamic and meteorologicalmodels.Works by [E1186] and [Nih84] both approach the problems of determining spread anddegradation from a different tack. [E1186] models the slick as a distribution of oil dropletswhich are churned in the ocean by the wind and wave state. Droplets are constantlyagitated into the water column, and are advected by the surface currents. The majorityof the droplets rise to the surface due to buoyancy. As a result the spreading of the slickis described by the distribution of the oil parcels and shear diffusion. [Nih84] developeda non-linear model of the transport and spreading of oil slicks. His work simultaneouslyaccounts for the effects of gravity, surface tension, friction, viscosity and weatheringprocesses. The parameterizations were adapted to field conditions and are reported tohave produced results which accurately modelled actual rates of spreading and advection.Chapter 5. Previous Work^ 42Other researchers who have investigated the broad area of modelling oil spill pollutioninclude [Tuc80, Lar88, Hun90, Yap89].5.3 Visualization of Previous SystemsNone of the literature reviewed discussed output methods for interpreting the variousmodels. However, the nature of the problem suggests that visual methods were mostcertainly applied when evaluating the models.Ahlstrom includes annotated line plots of the distribution of the oil spill through time[A1h75]. [Mac77] provides extensive graphs of the resulting process rates for evaporation,emulsification, drift and spreading. Additionally, the effect that volume and compositionhad on these rates is shown. [Ven87, Ven90] shows the distributions of the particles ratherthan the slick body over time. The effect of varying the rate of the turbulent dispersionis readily apparent through these visualizations. In [Ven90] the grounding of parcels isshown along the Alaskan coastline. [A1R89] provides the most detailed visualization oftheir models output. Contour plots of the oil distribution over the domain region overtime are presented. These displays gives a strong impression of the effect of the driftand spread on the distribution of the oil spill over time. Using similar methods thedistribution of oil in the water column and shoreline are shown.This thesis will take the modelling knowledge provided by these systems, and extendit by adding a robust interface, direct model interactions, and presentation of the modelstates as visualizations. The interface will provide a simple means for the user to employthe system. Further it will shield the user from the internals of the models, and allow fordata transparency and enforce organization on the system. Interactions will make thesystem both participatory and afford latitude in verification of models and of simulations.The visualizations provide a rapid means for interpreting the results, and are directedChapter 5. Previous Work^ 43towards a predictive goal of estimating the gross effects of an oil spill on the domain.This integrated system provides a tool for a broad range of disciplines, from ecology andregional planning, to computer science and geophysics, and yields a system which buildson the paradigm of modelling and simulating a marine oil spill.In this chapter previous work in the areas of modelling and simulating oil spills havebeen reviewed. The flavours of these systems range from the statistical and observational,to the strictly analytical. These techniques and methodologies were used in defining andimplementing the oil spill visualization system.Chapter 6Modelling6.1 ModellingThis chapter covers the derivation of the models which comprise the system. The spe-cific techniques applied to representing and modelling the environment and the oil spillprocess are described in detail. Further, the orderly flow of data through the system, andgraphical issues such as the visual properties of petroleum, and interpolation schemesemployed will be presented.6.1.1 Restriction to the Ocean's SurfaceThe system, and it component models, are designed to simulate only activities that occuron the air/water boundary, and at the land/water boundary at sea level. While many ofthe environmental processes extend high into the atmosphere, or deep into the ocean, forthe purposes of this research we will only consider the resultant net effect at the surface.This boundary does not form a plane, as the surface of the ocean is lifted and depressedby wave actions, the rise and fall of the tide, and the restoring force of gravity. Sincethe size of a domain can extend for many kilometers, and the vertical changes due tothe wave state are small, it will be assumed that we are operating on a plane. Theserestrictions reduce the complexity of the models in both time and space, but concurrentlyreduces the flexibilty and accuracy of the system.44Chapter 6. Modelling^ 456.2 Modelling the EnvironmentThe environmental models will consist of the terrain, tide, wind, and sea state within theconfines of the region of interest.6.2.1 The TerrainA representation of the terrain of the study domain was constructed. The topography ofthe Vancouver region is complex, formed by the Coast mountain range, the Fraser Rivervalley, and a number of fjords, channels and island groups. An accurate model of thesurface of this region was critical to many of the stages of the simulation and visualizationdevelopment.Two representations of the local terrain were developed. A three-dimensional rep-resentation was constructed for use in visualizing the region, and for derivation of theprimary representation which is used exclusively in the simulation.The three dimensional data for the region is collected from three sources, and ulti-mately merged into a irregular triangulation of (x, y, z) data. Most of the data for whichz > 0 1 came from a digital map-sheet of the 92G UTM block, which originated from thegovernment body Energy, Mines and Resources [EMR88]. The data set is comprised ofthematic boundaries and features for the region 2 , and contours at 1000 foot intervals.Data for z < 0 was provided by the Department of Oceanography. It is in the formof a two kilometer grid of depths represented in fathoms. With the inclusion of lowlevel spot heights 3 , all three sets were merged and triangulated using the functionalityof the ARC/INFO geographic information system. This system was used throughoutthe development of this system for all large scale tasks that involved the registration1 Sea level.2 11uman and natural features, such as political boundaries, major roads and highways, forest cover.3To detail portions of the region that are below the 1000 foot contour.Chapter 6. Modelling^ 46and transformation of large blocks of data. The triangulation was output into a generalformat for later use.6.2.2 The Land QueryThe second representation of the local terrain was formed by resampling the three-dimensional model. This resampling is for use in the simulation stage of the model.As we are operating only on the surface of the ocean, there we need to know if a givenpoint is on land, or on this surface. This question or query must be posed frequently inthe operation of the simulation and when sampling from the spatially dependant envi-ronmental models. As a result a criteria must be to answer this query in a fast and spaceefficient manner.The first method considered was a point-in-polygon test. The sea level contour (z = 0)was recovered from the three dimensional model, and constructed into polygons whichwere either land or sea respectively. For an arbitrary point within the domain, to de-termine if it is on land or sea by performing a point in polygon test. This test requiressome search be performed on the vertices, or the edges of the polygons. Preprocessingthe polygons into an efficient structure allows for the query to be answered in time whichis in the average case O(logn) [Pre85]. This cost amortized over a large number of itera-tions of the simulations is far to great. As the land/water query will be posed frequently,and as the contours that the polygons are constructed from may contain thousands ofpoints, this method was not considered a good choice for use in the simulation. Further,the need to extract subsets from the chosen representation exists. While intersecting,extracting and reassembling polygons from two sets of polygons is certainly solvable, thisapproach was not considered.The next method involves a further approximation of the available data in the region.The land-water polygons are labelled and constructed, and sampled into a regular gridChapter 6. Modelling^ 47with a resolution of 100 meters in both x and y. This step has a one time constructioncost. Further, it allows us to extract subsets without the need to reconstruct polygons.In an attempt to compress the size of the required grid, a quadtree representation wasconsidered. This structure could in some cases reduce the storage requirement of thegrid. However, to answer the general query would require traversal of the tree. Theaverage case depth of the resulting tree would be dependent on the complexity of thecoastline at the sampled resolution. In the best case, single queries would suffice whenan entire quadrant was contained in either land or sea. In the worst case, a branch ofdepth log N would be traversed, were N is the the maximum number of grid cells in xor y. Due to the inherent complexity of the shoreline in the region, this representationwould afford no improvement. In an attempt to address the excessive use of memory,a bit representation of the grid was implemented. As each grid cell is either on the seasurface, or not, it can be encoded by a single bit. For each row the cells in groups of eightare packed into a single byte. For an N x M grid, only Iii M + 1 bytes are required. Amapping of each point (x, y) into its grid coordinates (i, j) was implemented. This querywill be refered to as1 1 if z > 0L(x,y) =0 if z < 0This yielded an almost constant time cost for each operation, independent of the sizeof the grid. The major drawback of this representation is the loss of spatial resolution.However, the tradeoff is more than compensated for in performance.The grid method also gives us a compact and fast method of determining the inter-tidal zone. If the maximum amplitude of the tide in the region is determined, then theintertidal zone can be viewed as all surface features for which the position is constrainedby —C > z < ( where ( represents the tidal amplitude or height. If we extract theChapter 6. Modelling^ 48contours from the three dimensional model for z = —C and z = (, assemble the landwater polygons, and construct the grids L_c and Lc, then a point will be contained inthe intertidal zone if the following is true:L_c(x,y) and not Lc(x,y).6.2.3 Grid Data DescriptionA mechanism was associated with this grid which allow registration to the model coordi-nate system. This allows for a correspondence between the grid and the earth's surface. Itconsists of a description of the dimensions of the grid, and translation, scale and rotationparameters. This structure is described byG = {{n,m,h, 0,T, S}, G,}where n and m are the number of rows and columns in the grid, h is the length of thesides of the grid cells, T is the translation vector of the lower left corner of the gridcell (0,0) into the model coordinate system, 0 is the angle with which to rotate the gridcounterclockwise about the lower left corner of the grid cell (0,0) to bring the positivey-axis into alignment with true north, S is a scale factor, where it is assumed the gridhas the same spacing in x and y. G is an n x m grid, where each cell resolves to somestructure. In the case of the terrain grid, defined as G, the structure of the grid is a bit.However, as will be explained, this structure could contain a scalar or vector quantity.From this description any grid based data can now be registered to the model coordi-nate system. This concepts will be extended to all grids throughout the various models.These grids can have different spatial extents and grid cell lengths, and will not upsetthe operations of the system. The wind and tidal data structures, also based on gridstructures, can now be described in a similar vein, which allows for determining the cor-respondence between a wind or tide flow and the underlying surface. In addition, as theChapter 6. Modelling^ 49model coordinate system is based on the Universal Transverse Mercator projection, directregistration of this data to the earth's surface is possible. The contents of the elementsof a grid, can be any hypothetical structure, but are bits, scalars, and vectors for thisimplementation. A number of these grids will be defined. This grid description can bemodified slightly to account for data which is changing temporally as well as spatially.This synchronization data is required for the wind field data. The modified structure isdefined as follows:G = {{n, m, h, 0,T, S},G,, {T0,TN-1, N}}where To and TN-1 correspond to the inital and final sampling time (in hours), and Nrepresents the number of samples. The samples are evenly spaced over time. If thesimulation runs for longer than 24 hours, multiple copies of the sample set are used, andthe first and last samples are assumed adjacent.The grid description mechanism allows for the orderly flow of data into and out ofthe system. Any data set for which this format can be determined can then be employedby the system. Similarly, data used by the system can be transformed into other systemsby analysis of the grid description. The problem of data conversion into and out of thesystem is now external to the model.6.2.4 The TideDue the complex nature of the tidal generating force, and the interactions of forcedwaves with topography, no analytical solution to this system exists. However, there areempirical and numeric models which approximate the actions of the tide.Harmonic analysis involves the decomposition of observed tidal elevation data intoits frequency constituents. This method was the first quantitative method for predictingthe long term rise and fall of the tide. It is in general quite accurate, and is still theChapter 6. Modelling^ 50source of data for tide charts and tables. This method does not adequately describe thesurface currents.Modelling the tidal current from observed data is difficult, especially in regions ofcomplex topography. This is primarily a result of the difficulties involved in the longterm collection of surface current data, and in the effect of residual currents on the theobservation. The background noise introduced by wind and wave currents can obscurethe value of the tidal current. However, if enough horizontal current data exist, and noisein the measurements accounted for, then with an appropriate interpolation scheme anestimate can be made of the current distribution over a region. Many complex interpo-lation schemes have been applied to sparsely located spatial data. Ahlstrom used datacollected from the domain of interest using drogues 4 and fixed current meters to estimatethe tidal current in his domain [A1h75]. Using the observation that the ebb and floodtidal streamlines generally paralleled the shoreline, with small changes of direction fromflood to ebb, a base directional field for the tidal current was constructed. These wererepresented with line segments whose end points represent the direction at flood andebb, and the length of the line representing the relative amplitude of the current betweenthese given points. The amplitude of the tidal current was fixed to oscillate between thetwo end points of the line segment.With the advent of the computer, hydrodynamic models of the tidal forces and theresulting currents have been developed. The majority of the systems surveyed on oilspill simulations determined tidal currents either from a directly coupled hydrodynamicmodel, or interpolated the fields with data generated by such a model [Lar88, A1R89,Ven87, Ven90].The harmonic constituents of the major tidal forces in the Strait of Georgia basin andthe surrounding marine region were provided by the UBC Department of Oceanography.4 Free floating current gauges.Chapter 6. Modelling^ 51These constituents were extracted from a hydrodynamic model of estuarine circulationin the enclosed waterways of the Pacific Northwest developed by [Cre88]. The system inits entirety covers Puget Sound, the Strait of Juan de Fuca, and the Strait of Georgia(and most of the incident channels, sounds and island groups surrounding and containedby these water bodies).The existing three dimensional model represents the culmination of a series of tidalmodels developed over a number of years by Crean and associates. The initial modelwas one dimensional, and was used for determining the M2 tidal constituents throughthe domain. A grid was constructed over the region, with each grid cell representing across section of a bay, straight or channel. The connectivity of the grids was based onthe topology of the conveying channels. The equations of motion with friction presentwere solved using a finite difference scheme. Using many years of observational data, theparameters of the model were adjusted such that the energy losses of the waves in themodel and the real world were in approximate agreement. Through the years the modelwas extended into a two-dimensional and finally a three-dimensional scheme based on aregular two kilometer grid. For a detailed account of the methodologies and mathematicsapplied, see [Cre88].This provides us with a time-based analytical approximation of the varying tidalelevation and surface currents over the test region. To determine the velocity in the udirections at cell (i, j) at time t we need only evaluate the functioncE I3ijk cos(`2irfo + oiik)J...i.is the frequency of constituent k at cell (i,j)is the amplitude of constituent k6 See Chapter on Environment6The positive u direction is east and the positive v direction is north, j represents the evelation aboveor below mean tide.Chapter 6. Modelling^ 52Oink is the phase of constituent kC is the number of harmonic constituentst is the time in hours measured with respect to an absolute clockSimilar evaluations are required for determining v and C. The beauty of this model isthat it is stable, has a compact representation, and is quite accurate. [Cre88] reportedthat discrepancies in the height term c are less than five centimeters.The largest drawback of the model as a whole is that it registers poorly to the realworld. For a given grid cell, there is no direct mapping of the surface of the grid cell toits corresponding region on the earth surface. This had to be expected, as the grid cellsare aligned with the channels and other coastal features. A channel may be representedby a single grid cell in the model, while the actual channel may only be one kilometer inwidth. In addition to coarse spatial resolution, the model was developed at a time whencore memory was at a premium. It requires 7 that care must be taken when transformingdata from the grid cells to the earths surface. A subset of the tidal model was extractedthat covered the study area, and only grid cells that could be aligned with the terrainmodel by rotation, translation and scaling were chosen. This coarse registration wasperformed by hand using the resources of the Geography department.6.2.5 The WindDetermining the distribution of the local winds is the most important component of a oilspill prediction model. As was previously discussed, the wind largely dictates the overallfate of the spill. Its effect is felt in the drift, spreading, and weathering rates, and it isthe primary force driving the sea state.There are two primary methods of determining a wind flow over a region. These'The indexing scheme into the grid is not straightforwardChapter 6. Modelling^ 53methods consist of interpolation of observed data, and the solving of numerical approx-imations to the meteorology of the region. Both of these methods have advantages anddisadvantages. To correctly navigate through this area, both methods must be described.6.2.6 Empirical Wind ModellingWind velocity data tends to be sparse and irregular in space. Further, if the recordingdevice is not continuous, then the sampling will be sparse over time. Each set of obser-vational wind data represents the state of the winds for the particular conditions thatexisted at the time.A diurnal cycle over the domain of study was developed by interpolating a set of winddata from measurement stations over the domain. The data's source was EnvironmentCanada, and the date of the sampling was the 17th of July 1992. On this day a stablehigh pressure ridge had been in place for a number of weeks, resulting in typical daytimeand nightime temperatures for the region. The data was collected from 10 recordingstations spread out over the domain. The positions of these recording stations wereused to generate an irregular network. This network was interpolated using bilinearinterpolation at each of the sampling times, and sampled into a regular grid. With theinclusion of the grid and time description, it was included into the system.Given data sampled over a number of diurnal cycles, a more realistic approximation ofthe wind could be achieved. This problem is addressed in a modelling system developedby [Cek88]. Data from a three dimensional wind model was generated for prevailingwinds and indexed by eight directions and seven wind speeds. This modelled data couldbe replaced by observational data, although a larger bulk of data must be collected. Thismodelled data was interpolated into a regular grid and placed into the database. Usingobservations of the direction and magnitude of the prevailing wind, an estimate of thewind patterns at a given time could be interpolated from the database.Chapter 6. Modelling^ 546.2.7 Numerical Wind ModellingNumerical modelling of wind fields has been an active area of research for a number ofyears. It is an intricate science, from which promising results have been derived. Whilethese models have been developed to determine a wide range of atmospheric motions, thefocus in this thesis will be placed on mesoscale wind modelling. This area is concernedwith determining the wind patterns which are in the scale of one to 100 kilometers, andof which have a large influence on the spread and drift of oil spills'. Wind patterns atthis scale near the earth's surface are largely influenced by channeling and deflectioncaused by topography, friction, and exchanges of heat and momentum between the loweratmosphere and the surface of the earth.Mesoscale numeric wind modelling systems have been developed by many researchers,including [Pie74, Mas85, Ayo87]. Our discussion will focus on generalized models whichdo not consider all of the physics involved. The models rely on solving tendency equa-tions for meteorological flow. These tendency equations consist of formulations of thehorizontal transfer of momentum under the influence of friction, heating, the Coriolisforce and the pressure gradient. There are three types of generalized models currentlydeveloped for diagnosing wind fields over terrain.Mass conservative models rely on observational data for initial conditions, and theninterpolate these values into a regular grid, and apply least squares methods to theconservation of mass throughout the model.Multi-layer vertically integrated models divide the atmosphere into four layers. Thetendency equations of heat and momentum transfer are integrated through the layers.This method is sensitive to boundary conditions and assumes a well mixed atmosphere.The sigma-coordinate method uses one level in pressure coordinates to integrate the8See Chapter 3.Chapter 6. Modelling^ 55tendency equations. These models are not mass conservative. Changes in pressure aredetermined hydrostatically9 by parameterization of the atmospheres structure, along withboundary layer friction and diabaticA sigma-coordinate based wind modeller was was provided by Keith Ayotte of YorkUniversity and is based on [Ayo87]. The system was developed in the Department ofGeography at UBC. This model is of the third type, relying on a vertical parameterizationof the atmosphere. The initial conditions which the model requires includes a grid of theregional elevation, and the heights and temperature of a reference pressure leveln. Thefriction and heating terms are presented to the model as spatial grids.The necessary inputs for the Vancouver and region study domain were generated.To model Vancouver and its surrounds, an elevation grid was sampled from the three-dimensional surface model. Estimates of the Bowen ratio and the roughness length weretaken from [Ayo87]. Heating and cooling levels were adjusted in an attempt to estimatethe wind patterns that would result over a diurnal cycle. The model was found to bevery sensitive to grid size, and to the heating-cooling parameters. The minimum grid celllength for which the model would converge was found to be 18 km, which only coveredthe domain with an eight by six grid.6.2.8 Interpolation of Wind Field DataGiven that a surface wind field is available, the surface current due to the wind must bedetermined. The physics of the shear induced by the wind and the resulting current isdiscussed in the Appendix. Given a wind velocity -17 at a given point, the direction and9The hydrostatic assumption states that the vertical component of the pressure gradient force isbalanced by the force of gravity ie a = —pg where p is the density of the air and g is the gravitationalforce.'°Diabatic forcing results from an exchange of energy'For this thesis it was set at the 850 mb pressure level, which corresponds in elevation to approximately1500 metersChapter 6. Modelling^ 56magnitude of the induced surface current was given by,V* = crMe -1 .7.a = wind coupling coefficientMB = rotation matrixThe deflection angle Odef is set to 0 for the simulations presented. Spaulding reports in[Spa88] that many researchers have had success in modelling the transport of oil slickswhen applying a deflection angle of Interpolation of Wind and TideThe tidal and wind models provide reasonable approximations of the flows over regulargrids. The simulation will need to know the velocity at off grid points.Without loss of generality, interpolation of the flow in the x direction, u, will beconsidered. v and C are determined in the same fashion. A u surface will be constructedfrom a sampling of the flow in the x and y direction by using (x, y, u(x , y)) as the surfacepoints at time t. The grid points determine the corners of each rectangular patch. Eachpatch will be represented parametrically in x and y by s and t, which will range over[0, 1]. For a position (x, y), the patch it is contained in is determined, and the values for sand t found. u(s, t) can then be determined. Due to the mixed grids used in this system,some boundary effects result. The resolution of the land/water grid is much finer thanthe tidal grid. If a grid cell in the tidal model has no constituents, then the flow is set asa default to (0, 0). However, in reality, the flow due to the tidal current will be non-zeroup to the land/sea interface. As a result, there is disproportional weighting towards azero flow near this boundary. This problem is described in Section 6.1.Chapter 6. Modelling^ 57Ui- 1 j- 1WaterFigure 6.1: Bilinear Interpolating of Tide FlowThis interpolation scheme restricts the points of the grid to be on the surface, but it isnot C 1 or C2 continuous. Catmull-Rom interpolation, developed in [Bre77] was appliedto enforce smoothness on the surface. This bi-cubic parametric interpolation schemeenforces C2 continuity, and places the surface through the sample points. To determineu at a point (x, y) the patch it is contained in is determined. From (x, y) the values of .9and t can be determined. The value of u is found by evaluatingu = :§cmcTtTwhere (s) and^are the vectors [s3s 2s1] and [t3t2t1]. C is the Catmull-Rom kernelmatrix, and M is the geometry matrix formed of 16 control points from the surroundingpatches [Fo182, Bre77].Temporal interpolation of the tide field is not necessary as the model is explicitlybased on time. Interpolation over time is required for the wind field. This is achieved bylinearly interpolating between two sampled fields. The wind fields are sampled N timesover a diurnal cycle. For an arbitrary time t at (x, y) the flow is determined by usingbilinear interpolation to find the flow in the previous wind field, Vprey and the flow in theVcurrentv= VVprev + (1 — V)Iinextt — prey=Chapter 6. Modelling^ 58next wind field Vnext • The flow at time t is found by calculating Veurrent wheretnext — tprevwhere tprev is the time the previous wind field was sampled at, and t- next is the samplingtime of the next wind field.6.2.10 Other Advecting FlowsThe domain is sheltered for the most part from large scale geostrophic oceanic flowsoccurring in the Pacific [Cre88, A1h75]. The flow from fresh water runoff is significantin this predominantly estuarine region, largely as a result of the contributions from theFraser River. This fresh water runoff has an effect on the temperature and compositionof the water at various levels in the water column, which is primarily a vertical action.As a result it will not be factored into the system.6.3 Other Environmental FactorsThe state of the sea surface is important in determining the various process rates effectingan oil slick. It determines the turbulence and the amount of available mechanical energyon the ocean's surface. It is largely determined by the prolonged force of wind appliedto the ocean's surface.Qualitative description of the state of the sea have been made by people since theyfirst began travelling by water. Over the years these descriptions have been somewhatquantified in the Beaufort scale. This scale is calibrated against both the velocity ofthe predominant wind, and a description of the apparent wave state. For example, aBeaufort value of 7 defines a range of wind speeds 12 from 22 to 25 ms -1 , and is described12Around 55kh-1Chapter 6. Modelling^ 59as a moderate gale, with large waves, whitecaps and the heaping of the sea onto itself.The Beaufort scale provides a description of the state of the sea at a given moment,but does not give any insight into the affect of a change in the wind patterns. Waveprediction is based on the assumption that surface waves are caused by a wind forcewhich blows at a given speed for a given duration over a length or fetch of the ocean.A fully developed sea is the state that results when an infinite fetch wind blows for along period of time. This fully developed sea will require a stronger wind to attain largerwaves, and a weaker wind to diminish the surface.As an approximation to the state of the sea surface, a parameterization of the BeaufortScale was developed. It was formulated to respond to changes in the average wind speedover the domain. This is achieved over time by the formulationSt = V3(lUt l) + (1 — q)St-Atwhere i is a weight in [0, 1], MI is the average wind speed at time t, and St is the seastate at time t, and B is a function which returns the sea state number based on the windspeed. The sea state is dependent on the decay of the previous state, plus the energycontribution from the wind. Ranging the parameter q changes the effective contributionof the previous sea state and the present.6.4 Modelling of the Oil Spill ProcessThe models assembled will provide an approximation to the general transport equationap = _rjvp+v(ichvp)+ 0at (6. 1)where p is the density of the pollutant, 0 is the net advecting current, Kh is thecombined effects of dispersion and physical spreading, and 0 represents sources and sinksChapter 6. Modelling^ 60[Ven87]. This system can be approximated by solutions to the differential equations thatit is described by. These methods tend to be slow and potentially unstable. As theintent of this system is to provide an interactive mechanism for abstracting oil spills, aLagrangian particle tracking method is employed. This method is more efficient becauseit does not depend on the system converging on a solution, and is stable [A1h75]. It isalso more efficient in cases where a large volume of petroleum is involved, and wherethe particles tend to occupy only a portion of the domain [A1R89]. Each of the termsof Equation 6.1 will be determined and the state of the slick at each point in timedetermined. The first term —0V p represents the advection of the slick from the winddriven and tidal currents. The second term V(KhV p) indicates how much the slickspreads and disperses over time, and the final term, 0, represents the amount of oilwhich is added to the slick from a venting petroleum source, and the amount taken awayby weathering and grounding processes. These terms are determined at each iteration ofthe simulation.The initial conditions are contained in a representation of the oil spill source S.S . {Po,V0,C, f, to}where Po = (x, y, 0) is position, V is the volume, C = (p, M0 ...MN) is the compositionof the oil, and f is the flow rate, and t o is the time when the source will begin to ventinto the marine environment.The position Po can be located anywhere on the ocean surface. The volume Vo isexpressed in liters, and the flow rate f in liters per second. The composition C ofthe petroleum is reduced to its density p, and the respective mass of each hydrocarbonfractions Mo through MN. These parameters are used primarily when determining theweathering processes.The oil slick will not be modelled as a continuum. It will be represented as a setChapter 6. Modelling^ 61of parcels, each of which is independent. It is the nature of the Lagrangian particlesystem to consider the slick as a collection of distinct entities. Each parcel will carrywith it its volume, composition, time stamp, and a distinct label, and is processed by thesimulation as a micro-slick with these properties. The collection of particles is placed ina list structure, ordered by label, which will be refered to as LS.6.4.1 Oil ReleaseThe oil is released at Po at time to . For each iteration a new volume is released, until thesource is depleted. The volume that must be released is determined by V, = fAt whereAt is the time step. This volume is then allocated to particles based on a parameterVmax which is the maximum volume which can be allocated to a particle. The numberof particles released at a given time step by,Not = lv÷..xj +1 if Vs > Vmas1^if Vs < VmaxAs Not will likely not be an integer, at a release where Vmax << Vs the volume ofthe first Not — 1 parcels will have volume V. with a single particle of volume V. —Vmax (NA t — 1). The new particles will each be assigned a unique label, and placed ontothe slick list denoted Ls, and the volume of Vs decremented.6.4.2 TransportationThe transport of the slick at a given step is found by determining the tidal and winddriven current at each position currently occupied by a parcel. Duplicates of Ls arecreated, and passed to the wind and tide interpolation routines. The lists LT and Lwrepresent the tidal and wind flow at each point occupied by particles of Ls. The newChapter 6. Modelling^ 62positions for each particle are determined byPi+ ii = (T(Pii ) + aModef [W(Pii )])Atwhere Pi; is the jth particle position on list Ls at step i, T(P) is the tidal velocity atposition P, W(P) is the wind velocity at P, Rodef is a counter-clockwise rotation of thewind direction at the surface, and o is the wind coupling coefficient (See appendix A).The advection of the particle P is the superposition of the local tide and current due towind stress at that point.6.4.3 Turbulent Diffusion and SpreadingThere have been many empirical models put forward over the years for modelling thephysical spreading and turbulent diffusion that affect a petroleum slick.Early methods by [Fay71] concentrated on determining the physical spread from thespreading force of gravity and the retarding forces of viscosity and inertia. The modelsplits spreading into three regimes. The first regime involves spread under the force ofgravity, with inertia as the retarding force. In the second phase, the inertia is over-whelmed by the viscous force as the primary retarding force. In the final regime, theoil has become so thin that the effect of gravity is no longer felt. Instead the net forceresulting from the difference in tensions between the water, air and oil dominates thespreading. Mackay and Leinonen derive an empirical relationship based on radial spreadusing Fay's formulations for the radius of a slick over time in [Mac77].c/ 3 Ci tVo (pin — po )port = ro + 7rptuwhere ro is the initial radius, C 1 is a constant, t is the time in seconds, Vo the originalvolume and pw and po the densities of the water and oil respectively.Chapter 6. Modelling^ 63Methods which assume the physical spreading to occur radially must also presumethat it is occurring on a calm sea. These models do not account for weathering processesor turbulent shear.To model the spreading of oil in a more realistic way, both physical and turbulentspreading must be modelled. Methods for determining the turbulent diffusion of a slickbody either involve solutions of the differential equations describing this process or usestatistical methods, which assume the diffusion can be described by a random process.The turbulent diffusion model developed by [A1h75] is of the later type. It was selectedfor use because it acts upon discreet quantities rather than a continuum, is rapid andstable, and is easily parameterized. This model assumes that the particles diffuse underthe influences of near Brownian motion resulting from turbulence in the transportingmedium. Diffusion occurs in both the vertical and the horizontal plane, although we willonly consider the horizontal term. A parcel will travel on average a distance describedbyDh = ^47hAtwhere At is the time step, and -yh is the horizontal dispersion coefficient which is acombination of the eddy diffusivity and the spreading diffusivity For an arbitrary parcelthe diffusion step size is determined by a random process. A uniform distribution on[0, N/127d is used to generate these steps. A direction for the step is determined from79 = 271-[4. A parcel will be dispersed by the amount(Dh cos 19, Dh sin 19).For each iteration of the simulation, the diffusion of the slick is determined by assigning arandom diffusion step to each particle. These steps are generated for each particle of Lsand placed onto a new list labelled Ldif , which will eventually be added to the positionsof the particles on Ls. The value of the dispersion coefficient -yh has been reported inChapter 6. Modelling^ 64the range of 1 to 600 m 2 .9 -1 and is mostly dependent on the sea state.This models the diffusion, but physical spreading must still be dealt with. Manyresearchers have observed that oil spreads in an elliptic fashion. Al-Rabeh and associatesuse the following formulation to determine the elliptic spread of individual particles andthe slick body as a whole in [A1R89]. This formulation is shown in figure 6.2A = 7-1-QRPw — Po^1 1Q = C(^)3i VIt4PoR = Q - F aWtttwhere Q and R are the minor and major axis of the slick ellipse, V is the initial volume,and W is the wind speed, and a is the wind coupling coefficient. R is oriented in thedirection of the wind. This spreading is applied to both the individual parcels and theslick as a whole. The slick is considered as a nest of similar ellipses. If we consider thewind to be blowing in the positive x direction, with the center of mass of the spill locatedat the origin (as in figure 6.2) then the outward spread of the individual particles isdetermined bydx = Sixdtdy = Syydt,Si = xdRRSy _ ydQQ •The model employed in the simulation system uses a combination of radial and ellipticspreading. However, elliptic spreading is only applied to the entire collection of particlesto determine the overall spread. The physical spreading at a particle is assumed to beChapter 6. Modelling^ 65Individual Parcels,P/S. ,.^ .. Boundary at t.N......^ .Direction of Wind. Displacement of ParcelBoundary at VI-oft^ By Elliptic SpreadFigure 6.2: Spreading of Model Slick as Concentric Ellipsesradial. This was done to avoid calculation of the long and short axis of the hypotheticalellipse on a parcel by parcel basis. This can be justified based on the observation that theelliptic spread of the individual particles will not greatly affect the overall distribution ofthe oil. Rather, it is the elliptic spread of the entire oil body that will determine the shapeof the slick. One question that remains to be answered is the effect on this algorithm ofthe point selected as the center from which to spread from. At present this is determinedby calculating the center of mass of the particles. However, in the presence of strongcurrents and winds, the center of mass may not even be near any parcels. The spreadingwill often draw the slick bodies out into frayed and linear shapes. Spreading about thispoint seems a bit simplistic. The spread about a medial line could help to relieve thisproblem in certain cases. However, the real problem comes from the representation of theChapter 6. Modelling^ 66slick body as a collection of particles. There is no fast and concise method of determiningthe shape of the slick from the particle representation. The slick may not even be a singleentity, but could be formed of multiple disjoint slicks which separate and rejoin over time.This elliptic spread model can only be applied consistently if the spill vented more orless instantaneously. Further, its effect should be modified when the slick body beginsto elongate in the presence of the local currents.The only other consideration that must be made is the duration of the spreadingapplied to a slick. This was done by halting the spread at the particles once the thicknesshad reached limin, = ,7-.4 = 10' meters, or a hundredth of a millimeter. While many spillsachieve slick thicknesses which are monomolecular, most large spills are reported to havestopped spreading far before before achieving this thickness.To summarize, the spreading process is modelled by applying the diffusion algo-rithm developed by [A1h75], in conjunction with a modified elliptic spread function from[A1R89].6.4.4 Evaporation and DissolutionIn order to determine the evaporative and dissolution flux from a slick body, either apsuedo-component or analytical approach must be used. [Mac77] proposed a psuedo-component approach, which separated the oil composition into its main fractions. Eachcomponent is described by its respective density, boiling point and other parameters.To determine the evaporation flux for the entire oil mass, the flux of each component isdetermined, and summed over the entire massN NN = E IV: = EICe x iP:RTi=1^i=1Chapter 6. Modelling^ 67where Arf is the evaporative flux of component i, ife is the evaporative mass transfercoefficient 13 , xi is the mole fraction of component i, Pia is the vapour pressure of thecomponent i, R is the gas constant and T is the air temperature above the slick. Themass lost to evaporation proceeds into the atmosphere, which is assumed to be an infinitesink. Other methods avoid decomposing the oil body, relying instead on averages to theevaporative flux. These methods have the advantage of requiring less information. Theevaporation rate has been described as a function of area, wind speed, slick thicknessand the roughness of the sea surface.For the purposes of this system, the composition of the oil will be generalized, andMackay's algorithm applied. This algorithm is selected because it considers the majorfactors involved in evaporation, and is compact and stable. The ambient air tempera-ture will be considered constant for the simulation. The components of the oil will begeneralized into a volatile and non-volatile fraction, by boiling point. If a component hasa boiling point lower than the ambient temperature, then it is considered to be volatile.The vapour pressures are averaged by mole fraction.6.4.5 GroundingWhen a parcel drifts into contact with the shoreline, a decision must be made as toits fate. In general the entire mass of the slick will not ground. Only some fraction willground, with the majority of the mass remaining in the ocean surface to continue to drift.The mass that sticks to the foreshore is usually a function of the slope and roughness, andthe composition of the oil. The properties of the shoreline are assumed to be uniform.An arbitrary splitting rate w is set. If parcel Pi with volume V comes in contact withthe shoreline, (1 — w)Vi remains in the parcel, while wVi is assumed to have grounded.The mass of each water-borne slick is monitored, and if it goes below a threshold, then'The mass transfer coefficient is proportional to lUv where v < 1 and IUI is the wind speed in ms-1OriginalDisplacementChapter 6. Modelling^ 68the parcel is removed. The volume that has grounded is kept on a separate list LG. Todetermine if a slick grounds, the superposition of the effects of advection, diffusion andspreading are determined in the vector U. If the addition of this vector to the position ofthe parcel results in a grounding, then a grounding parcel is split off and placed at P + U.The new position of the main parcel is determined by a random process. The vector Uis alternately rotated counterclockwise by .6.0 and constricted by AIM. AO is a randomirvariable sampled from a uniform distribution in the range of [41 . The constriction ofthe vectors length is generated from a [R]lou l uniform distribution. These parameters wereset to these ranges to minimize the number of steps required for a parcel to find a pointcontained on the ocean's surface. This process is shown in figure 6.3.Figure 6.3: Particle-Shoreline Avoidance6.4.6 Weathering Processes not ModelledThere are many processes which are not modelled in this system. Most of the effort wasconcentrated on determining consistent advection, diffusion, spreading and grounding.The primary reason for not developing a more-in-depth weathering model is the levelof detail required to attain a truly consistent system. Complex chemical descriptionsof the petroleum source are required. Methods for determining emulsion formation,Chapter 6. Modelling^ 69dissolution, and vertical mechanical dispersion are both complex and time consuming,and are discussed in [Hua83, Spa88, Mac77].Processes such as biodegradation and photo-oxidation are not modelled because theyhave a small effect on the mass of the slick, and occur at slow rates late into the spill's lifecycle. Further, the processes that drive them are poorly understood, which is reflectedin the lack of available models.6.5 Determining the Time StepAn appropriate time step must be determined for the simulation. Too large a step willdestabilize the oil spill processes, and in particular the diffusion and spreading. Further,the temporal variations of the environmental conditions must also be considered. Toosmall a time step will cause poor performance time. [A1h75] states that if smooth andcontinuous performance is required, then the time step should be limited such that themaximum displacement of a parcel over time be less than the minimum distance betweenthe data points in the velocity field. In our case, these separations are not fixed due tomixed grid sizes. As an alternative the maximum distance that a particle may travel isconstrained by the parameter Dmax. The time step is found by determiningAt =  DmaxiVi.where 11/1„„„ is the maximum speed of the total advecting current. The time stepsgenerated by this constraint are determined on an iteration by iteration basis with goodresults. It was found that the system behaved erratically with large, fixed time steps.As a result of this adaptive time step, the absolute clock is determined byktabs = to + E tii=iwhere k is the iteration number, and ti is the time step determined at the ith iteration.Chapter 6. Modelling^ 70In this chapter models for representing the petroleum, environment and spill havebeen presented. Further, representation and interpolation of data over the spatial andtemporal domain have been elaborated. A subset of these models were selected, developedand integrated into a the oil spill visualization system.Chapter 7The SystemGiven the collection of approximations to the various physical and chemical processes in-volved, systems must be developed for combining and managing these models to producesimulations and visualization of the model problem. These systems consist of a database,a simulation/visualizer and a batch or extended visualizer.7.1 InterfacesGraphical user interfaces were developed for both the database and the visualization sys-tem. It was assumed that a potential user of the system would be aware of the problemsassociated with oil spills and the associated processes in the marine environment, andhas some knowledge of interactions with a computer system. The criteria for develop-ment was to provide a front end that would allow for not only the rapid developmentof model simulations of oil spills, visualization of this event, and variation of key modelparameters, but of more mundane tasks, such as extracting subsets, and saving modelstates. This was facilitated by developing a system using the graphics language GL 1[GL90] which runs under the graphic user interface software of FORMS [Ove91]. Imple-mentation details can be found in the appendix. This application software allows for theinteractive construction of menus, with the inclusion of features such as buttons, sliders,and other iconic input mechanisms familiar to users in general.A hierarchy of menus is developed, each with a strict purpose, such as to visualize tidal'Developed by Silicon Graphics71Chapter 7. The System^ 72and wind flows, or to set up reference grids or other display oriented features. Throughthese menus, many of the tasks involved in setting up and synchronizing the systemremain transparent to the user. The major drawback of the FORMS software is that itextensively uses the Silicon Graphics font manager. This effectively binds FORMS, andthe database and visualizer to operate only on systems which support this software.The extended visualizer is batch oriented. It was developed under the assumptionthat a user would eventually develop animations and detailed renderings of the variousprocesses and oil spill simulations. By using a simple script language, the user can specifythe location and type of data that is to be processed, the output required, how manyframes to generate, and the relative time from which to base these frames. These scriptfiles can be run as background processes.7.2 The DatabaseThe database contains the state data for the tidal model, the grid representation of theterrain, and wind flow data. The front end constructed over the database shows the userwhere the wind and tide models are defined spatially, and allows for the extraction ofsubsets of this data. It was deemed reasonable to provide the functionality for subsetextraction for three reasons. The first is that the overall wind, tide and terrain modelsare quite large. Execution of the entire system is a strain on the fastest of the availableworkstations, at least in the spirit of trying to provide a somewhat interactive system.The second reason is based in the scale of the problem. If a user wishes to spill tenmillion liters of oil into the Gulf of Georgia, then a large portion of the entire modelmay be extracted. However, if they wish to study the effect of spilling ten thousandliters into Vancouver harbour, and to determine possible grounding of oil in this region,then they need only extract that portion of the model which pertains to that study. TheChapter 7. The System^ 73option, however, is available to work on the entire domain or a subset. The final reasonis security. It is a dangerous precedent to allow users to be directly operating on a rootset of data. The possibility for corrupting portions of the database could occur, whichwould then require the administrator of the system to reconstruct the data sets. Instead,the user will select a section of the database, which is then extracted and placed into alocal directory in the user's workspace.7.2.1 Data Storage and ExtractionThe tide is an analytical model, and is spatially tied to grid cells covering the region.This model is time based, and can be sampled over all time to determine an exact valueof the tidal flows and the elevation at the given location. The numerical and empiricalwind models produced data sets which did not have such compact representations. Thenumeric model requires large amounts of processor time before it converges on a solutionfor a given set of initial conditions. As a result this model cannot be directly coupled tothe interactive simulation and visualization system. If it were, it would not only grosslytax the workstation which was running the system, but also would require the user'sattention, in terms of providing initial conditions, and sequencing fields. This wouldmake the inclusion of wind field data far from transparent. As an alternative, the windfields are precomputed, with synchronization embedded within the data representation.The conditions under which these fields are generated is certainly of concern to the user,but it will be a concern independent of the oil spill system. Empirical wind data wasused for the simulation as the numeric model had too coarse a resolution. The wind fielddata is generated at some temporal rate, and presented to the database in the form of agrid description, synchronization data on the various fields, and a set of wind fields. Thedatabase merely extracts a subset of these windfields, and recreates the local grid andsynchronization data in the users copy of the data. The components of the tidal modelChapter 7. The System^ 74and terrain grid are extracted in a similar fashion.This methodology provides users with a consistent means for the inclusion of datafrom other systems. The actual current field data must be placed in the system's regulargrid format, and a spatial and temporal description provided. In the case of presentingnew wind or terrain data, this would be achieved without any changes to the software.Only minor changes would be required to include a new current, such as a river flow or ageostrophic current. Changes in software would be required to replace the current tidalmodel, which has a functional representation, rather than one based on sampling.The data for the system is not directly loaded by the database software. Rather, thegrid description of the data is used. When a user selects a region, all sample grid pointsthat lie within this region (with data padding necessary in some cases due to boundaryeffects) are extracted. The subset is rectangular, and oriented with the horizontal edgesof the rectangle parallel to lines of latitude, and the vertical edges aligned parallel tolines of longitude. The orientation and shape restriction allows for the rapid extractionof subset data.The front end of the database provides little additional functionality, other than theability to zoom and pan around the region, to change the path with which to store thedata, and the naming of a given data set.7.3 The Interactive VisualizerThe interactive simulation and visualization system allows the user to develop, execute,and visualize oil spill scenarios, save the state of these simulations, and to visualize thesurface flows due to the tide and wind interactively. Descriptions of the models, samplingschemas and the flow of the simulation are more extensively covered in chapter 6. Thesystem is presented to the user in the form of an input menu, and a graphic display ofChapter 7. The System^ 75the current domain. The interface is a hierarchy of menus, each performing a specifictask. The route taken by the user and the order of operations is relatively free flowing.Domain selection is performed by selecting one of the many domains previously ex-tracted from the database. The user can rapidly peruse the available regions from abrowser, with the graphic display updated to reflect the current region. The structuresand data of the domain models are not loaded into the system until the user explicitlyperforms this operation. From this junction the user can either execute a simulationwithin the confines of the loaded domain, or visualize the major flows.7.3.1 Running a SimulationThe simulation panel allows the user to set the initial conditions, sample from, andexecute a oil spill simulation. The initial conditions consist ofVo =Vp =f =Po =to =inital volume of the sourcemaximum volume per parcelflow rate in liters per secondposition of the source within the domainstart time of the simulationThe volume, parcel volume and flow rate are set by manipulation of a slider. Themaximum volume that a parcel may contain dictates the number of parcels which will bereleased. If this volume is set to a large value with respect to the volume of the source,then few parcels will be released. The performance of the simulation is proportional to thenumber of parcels released. Sensible bounds are enforced on these sliders. To determinethe initial position of the source, the mouse icon must be placed over the correspondinglocation in the graphic display and the appropriate key depressed. This will place aChapter 7. The System^ 76symbol permanently at that location within the graphic display. This operation is notpermitted either outside the domain, or on land. The starting time of the simulation isprovided through the time menu. It consists of parameters for year, month, day, hourand minute, within appropriate bounds.There is also a set of model parameters external to this menu. There adjustment ispossible though a parameter file that is placed in every users workspace, or changes tointernal variables in the program. The setting of their values is not strictly speaking arequirement for running a spill simulation. Further, they are used in other stages of theoverall model. A description now follows of these parameters, based on their purpose.Te = tidal component switchesa = wind coupling coefficientQT = tidal flow constraint coefficient17 = sea state weighting parameter7', determines which of the harmonic constituents of the tidal model are to be used inthe current simulation. The model includes 14 major constituents, which can be turnedoff and on through T. a is the wind coupling coefficient. It determines the magnitudeof the drift current due to the wind, and is usually set in the range of 0.01 to 0.08. ayis analagous to the wind coupling coefficient. It has no physical basis, and is merely inplace to allow for adjusting the velocity of the tidal current. The parameters a and QTcan be employed to adjust the effect of the wind and tide, if necessary. The weightingterm 9 is used in determining the sea state. It is constrained to range between 0.0 and1.0. A small value for 9 places more weight on the previous state of the sea.Chapter 7. The System^ 777h =Po =Pw =w =Dmax =ATmax =horizontal dispersion coefficientpetroleum densitywater densitysplitting rate for grounding parcelsThe maximum distance a parcel may travel in the advecting flowThe maximum time step allowed for each iterationThe parameter ryh defined in Chapter 6 determines the magnitude of the step sizethat a parcel may undergo under the influence of horizontal dispersion. As it is largelya function of the upper layer turbulence in the ocean, it will be combined with thesea state in Chapter 8. The terms pa and p„, the densities of the oil and water, willdetermine the spreading rate. The parameter w will determine how much of a parcelwill ground when its advection, spreading and diffusion brings it in close proximity tothe shoreline. The parameter Dmar is used to determine the length of a time step on aniteration by iteration basis. The term AT,,,„„ is set to avoid large time steps. Large timesteps cause erratic and visually unappealing events to occur in the simulation. Otherparameters and combinations of parameters will be discuss in Chapter 8 and in contextof the visualization process.7.3.2 Displaying the SimulationOnce the setting of parameters is complete, the user may execute the simulation. The slickis represented as a collection of parcels venting from the source symbol (Figure reffig:eo).For each iteration of the simulation, the parcels are vented, advected, diffused, spreadChapter 7. The System^ 78and potentially grounded. The rate at which the iterations are displayed is dependent ona frame rate. It determines how many iterations to perform before updating the display.Setting this rate sufficiently high speeds up the execution rate, but will appear veryjumpy and discontinuous. Grounded parcels are differentiated from sea borne parcels bycolour. The grounding position is recorded, and redisplayed for each preceding frame.The discussion of this visualization is expanded in Chapter 8.Figure 7.1: A Model Slick as Viewed in the VisualizerAn environmental monitor can be switched on during the simulation. This monitorreports the general state of the environment in a noninteractive panel separate from theactual graphic display. It provides a synoptic view of the state of the environment on anChapter 7. The System^ 79iteration by iteration basis. The current time, and the direction of the prevailing windare displayed as text. The direction of the wind is in reported as a cardinal direction 2 .The magnitude of the wind is reported in meters per second, in the form of a vertical bargraph. The elevation of the tide and the sea state are reported in a similar manner. Theupdate rate of this monitor can be set in a fashion analogous to the frame rate. Thismonitor is a major tax on the available resources of a work-station. To achieve maximumperformance it must be turned off, or its sampling rate set very low.7.3.3 Sampling from a SimulationA key-frame can be defined as a sampling of the system taken at a particular time. Thissampling is sufficient to generate a graphic image of the simulation. A key-frame optionis available on the simulation panel. Prosaic tasks like specifying the path and names forthe resulting slick scripts must be specified. The parameters for sampling areNk = number of key-frames to sampleAt = time step between framestk = time offsetThe system attempts to perform the sampling in steps of size At from initial timeto + tk to to + tk (Nk — 1)At. However, because of the adaptive time step in thesimulation, this is not always possible. It will sample when the clock exceeds or is equalto the next sampling time. Each sample is placed in a script file with an extensiondenoting its position in the series 0 ... Nk - 1. Each script file contains the time of2A cardinal direction is one of the eight points on a compass, ie N,NW,NE etc.Chapter 7. The System^ 80sampling, a default eye position', and the identifier, position, volume and radius of eachparcel, separated as either being seaborne or grounded.7.3.4 Flow VisualizationThe flow visualizer allows the user to view the flows generated by the tide and the windmodels over the currently selected domain in an interactive manner. The input panelallows the variation of the sampling, scaling of the velocity field, and selection of a stylein which to view the flow.The sampling of the flow is performed by laying a regular mesh of N x M samplepoints over the domain. The dimensions of this grid are provided by the user, and Nand M need not be equal, but necessarily N, M > 1. The temporal sampling rate canalso be set.The tidal flow in the domain rarely exceeds 3ms -1 . The velocity of the wind canreach 30 — 40ms -1 if gale force winds were modelled. However, in general velocities donot exceed 10ms -1 . When flows of this magnitude are viewed within the framework ofa domain which may be 100 to 200 square kilometers, the dynamic patterns of thesesystems will not be apparent due to poor resolution. If the domain is shrunk, however,the gross spatial patterns of the flow can be lost. To address this problem, a scaling sliderfor both the wind and tidal flows is provided. The scale of the two flows can be adjustedindependently. This scaling allows the viewer to observe the more subtle interactionsoccurring in and around the channels and other oceanographic features due to the tide.The type of symbol attributed to a point flow can be changed. The most commonmethod of viewing a surface flow is by representing each point flow by a vector. For eachgrid point in the mesh, a vector which represents the magnitude and direction of the'The eye position is based upon the viewing parameters that were set in the interactive visualizerwhen the frame was extractedtk• ..3'^3 --....■_■- -,"4.-■____,,...,6• e-■^ri■-.././1■—.4„,„,,,i.,■■ 1.- ■'s.,-" ,,ete. /..-4- 4 ,..e.1•■;--...re//i."-, 1 1 1/./....,i, ,^1 1 it I 1 1! 1^I .,^ wor , ..4,..,. „.1 1 .1^■I i i 1^% ‘4.• , ..e.:. .". 4, si■NiNAAN„^  / ....wAg..._•-....,....j...■^\ \ NV. \l`•---'...■....;\^\-.4's\\*\i^'4"...:%"^NI"\7"..'__+-■-.■.■•^'''''''''\%\"\l',^, •-1,.. . .. k• ^ .., ,. '1, ..... . 0, -N-4-k A. •„ ift • , 1rrChapter 7. The System^ 81Figure 7.2:^Tidal Flow, Time0900:17-Jul-92Figure 7.3:^Drift Method, Time0900:17-Jul-92flows is rendered. The base of the vector is the sample point on the mesh, with the headdetermined by the addition of the local flow to the sample point. To view some of themore intricate flows in action, the sampling mesh must be increased in resolution. Thiscan however lead to relatively confusing images, especially in areas of slowly changingsmooth flow, caused by the overlapping of vectors. To address this problem, the flowcan be displayed with only the head of the vector displayed. This displays the flow as acollection of offset points.The third method approached both sampling and rendering of the flow in a differentfashion. At each sample point in the mesh, flow samples are released at a fixed timeinterval. It is allowed to advect with the current, and at the next time step, anothersample point is released. Both samples again undergo the process of advection. If theinitial sample point is included, a path of length 2 can now be formed between thesample points. This process is analogous to standing in a vessel at the sample point, and^••Chapter 7. The System^ 82..------,,,-''..-*--1‹..e'...-..-"' ----'^--.r^ ---,- ..,-- 4._ .......:;01-7.P''e""."-rr.--...r^..-4"" —4,4* ..,--4''—?--,— ---- ..--- ■!' '/. ''..,,^r^,,,,,,,,. „.„.,10 ......., .e,,... ,..,,,,er. ", ":orer-je'r.,,,,#.1,-1.er"..... 1-,..4?,...a..4.—••....--4,P#rP,..'<•^4'...^..■^"'"°:,........,••••.F,'4SP.......:*< "^....''''':■riel^.....,'"■•••■0"...,•4,.....r"1.,,rd''''^•■■,•°''..^....:.."..$...,...r:.........":„....,,,#,....r,,■+!•.......r...••••••4.......ia.r.'"'''''rd....r'r■..'''ro•4....r• 'msfi^....... ■.... ....412}".4....0^...... .•Figure 7.4: Wind Flow, Time 0900:17-Jul-92dropping a marker buoy off the side, attached to a string. After a brief period of time,another marker buoy also attached to the string and released. This process is repeated,and over time the pattern traced out by the string of buoys describes the gross patternsof the tidal flow. A similar analogy can be drawn in with describing the wind flow.This method clearly shows how the tide and wind are changing in both magnitude anddirection over time. Examples of these visualizations are shown in figures 7.2, 7.3, and7.4. The time of these snapshots and all other environmental images corresponds to thestandard spill simulation parameters discussed in chapter 8.Chapter 7. The System^ 837.4 The Extended VisualizerThe purpose of the extended visualizer is to produce detailed renderings of componentmodels, and oil spills simulations. These renderings are more computationally intensivethan the simple visualizations presented in the interactive system.To allow the user to communicate with the extended visualizer, a simple script lan-guage was developed. Each script file contains information on the location of the dataset, and the name to give to the resulting images. This name will be suffixed with aframe number. Next, the type of rendering that is requested. At present there are sixtypes supported.1. Render the elevation of the tide.2. Render the magnitude of the horizontal tidal flow.3. Render a thickness map of a set of slick scripts.4. Generate a slice of the slick at a given thickness.5. Generate an Optik format file.6. Interpolate between two slick scripts.These rendering types will be discussed in full detail in the next few subsections. Theresolution of the image to be rendered (where applicable) is specified by (N, M, D) whereN and M are the rows and columns of the result image, and D is the number of bytesused to encode the pixel value. At present D = 1, or 8 bits for the current set, butextensions to include colour could be included easily. The start time of the initial frameand a time step are provided. In the case of tide-specific rendering, the image is rendereddirectly from the tidal model. In slick-specific rendering, the time is encoded in the slickChapter 7. The System^ 84script. Both time and time step can be set to —1 as a signal to the system to used adefault time and time step. The final data that must be provided in the script file isthe start and end frame number. In the case of rendering slick data, a slick file with thesuffix must exist.7.4.1 Determining the Thickness MapIn the interactive stage of the system, the user views the slick as a collection of parcelswhich are acted upon by the local forces. This is a discrete approximation of how an oilslick disperses. The spilling of oil into water is the mixing of two vastly different liquids.As the oil is higher in viscosity and lower in density, it will suspend on the surface, andwill resist mixing with the surrounding ocean until a force is applied. If a slick wasviewed on a calm day, it would remain a single body, spreading radially. In the presenceof stronger local forces, the slick would elongate and break apart into a number of disjointslicks, which may eventually recombine. With very strong local effects, the body maydisperse in particle form throughout the water column, retaining no shape whatsoever.The parcel approach can be viewed as sample points affixed to the slick body, whichspread, advect and disperse in the same manner as the underlying collection of slick bod-ies. The task now is to recover the shape of the oil slicks, and determine an approximationof the distribution of the oil over the sea.As a measure of distribution over the surface of the ocean (and along the shoreline),the thickness of oil at a position is determined from the parcels. Each parcel is describedby its position P, volume V and radius r. A regular mesh of p x q samples is generated overa rectangle which contains some collection of parcels. At each point in the sampling regiona thickness of 0 is assigned. The parcels which represent the slick are then traversed,with the contributions from each parcel being placed into the mesh. The subset of meshpoints that will receive a contribution from a given parcel partially or wholly containedChapter 7. The System^ 85in the rectangular region will satisfy:iPi; — P1 < 7-20 < i < p — 10 < j < q — 1where r is the radius of the parcel P. At each mesh point the contribution of the parcelto its thickness will be accumulated. The thickness at a point is determined by thedistribution of the depth over the radial area. If the parcels have constant thickness,then the thickness at any point within the parcel would be:h— Virr2If the distribution of petroleum over the parcel was conical:h = ho (r — AOVho = 3 irr2Or = iPi.i — PiAfter the parcels have been traversed each location of the mesh has the accumulatedslick thickness recorded. This information is encoded into the p x q image. The thicknessis encoded into a single byte, using a logarithmic encoding scheme. This schema wasdeveloped to preserve the large changes in thickness over the integrated slicks (Figure7.6).The sampling rectangle by default is determined by the intersection of the plane z=0and the four sided pyramid emanating for the viewers eye position (see figure 7.5). Thesize and position of the sampling rectangle can be altered by altering the script file.It is assumed that the thickness is constrained to be in between the range of 10a < 10 bmeters where b > 0 and a < 0. The maximum thickness difference is b — a orders ofChapter 7. The System^ 86ViewpointFigure 7.5: Viewpoint and Sampling^ qRegion Figure 7.6: Sampling of p x q Gridmagnitude. If the thickness t = 0, then 0 is placed in the image. If t > 0 then thepixel is encoded into an 8 bit byte as Wlogio (t10-a) where W = (28 — 1)/(b — a). Thethickness can be recovered by reversing the operation. There will be some loss minor lossin precision, but the range of thicknesses will be preserved. For the system, b = 1 anda = —5, on the assumption that the slick will never be thicker than a meter, nor thinnerthat a hundredth of a millimeter. Both a, b, and W can be changed.This output image can be used both for encoding the thickness of the slick over thesurface and as a visual representation of the slick in the form of a grey scale image whereintensity is proportional to thickness. The boundaries of the radial parcels can in somecases be quite pronounced. To remove this effect a low pass filter can be passed over thethickness map before encoding.7.4.2 Interpolation of the Slick Between FramesGiven two frames of an oil slick sampled at different times, it is necessary to be able tointerpolate an intermediate state. This can be done either by performing a metamorphosisChapter 7. The System^ 87between two images, or by direct interpolation of the script files from which the imageswere generated.To attempt a metamorphosis from one thickness map to another, extraction of theboundary of the slick is necessary. Given the set of pixels of the respective thickness maps,an attempt is be made to find the corresponding positions of these pixels on a frame byframe basis. To accomplish this, methods of surface reconstruction were considered.Contours of equal value (in this case slick thickness) are determined from two thicknessmap separated by some time step, and an attempt is made to match them. From thismatch, an intermediate contour can be determined. Repeating this process over somerange of contours on both thickness maps, an intermediate state can be interpolated andreassembled. These methods are discussed in [Slo87]. The main problem is the matchingof points from different contours. The problem reduces to a search on a class of graphscalled toroidal4 .Contours were extracted from the slick image. The zero contour (the boundary be-tween no pollution and some pollution) was extracted. Each pixel in the image was setto 0 if it was void of pollution, and 1 if any trace of the pollutant existed. This provides anoise free high contrast image of contaminated and uncontaminated areas in the region.From this the boundary of the slick may be determined either by contour tracing oredge detection. Contour tracing, such as discussed in [Pav82], involves determining thechain of pixels which form the boundary of a patch of given intensity in an image. Edgedetection seeks to determine the edges or curves where rapid changes in the brightness(thickness in our case) or in the derivatives of the brightness occur. Given that we aredealing with a binary image, the two operations can be seen as more or less equivalent.4A toroidal graph consists of two collections of vertices Ci and C2 where the edges between verticeswithin a given set C1 or C2 form only cycles. These cycles form the contours. Edges between the setsC1 and C2 form the match between the contoursChapter 7. The System^ 88The Laboratory for Computational Intelligence has a developed set of image process-ing tools called VISTA. This software was used as an alternative to developing softwareto perform image processing tasks. The grey scale image was generalized into a binaryimage, with non-zero regions denoting regions contaminated by oil. The binary imagewas converted into the LCI format, and the components of VISTA used to extract theboundaries. As the image was binary and synthetic, no noise was present. As a resultlow pass filtering was not required. The edges were found by application of a Canny edgedetector [Tor82], followed by a routine which linked up edge pixels. These links wereconverted into an generic ASCII format. The recovery of edges was quite good, withvery few breaks occurring in the recovered boundaries.Unfortunately, attempts to enforce correspondence on a frame by frame basis couldnot be accomplished without human intervention. This was due to the fact that thenumber of components from frame to frame was not fixed. Applying the methods ofsurface reconstruction from collections of contours with this property tends to result inself intersecting surfaces. This results in a many to many relationship from N to Mcontours, where N M. Further, as the slick was the result of the reinitegration of a setof circular parcels, each one could potentially represent an independent slick (in practicevery few parcels were found to be distinct from a main body).Therefore, effort was redirected to a direct interpolation of the slick scripts. Overtime, the slick scripts will consist of seaborne and grounded parcels. No interpolation isrequired with grounded parcels, as their positions do not change with time. The uniqueidentifier assigned to each parcel will now allow for correspondence between frames.Care must be taken as parcels are added over time, the source continues to vent into theenvironment, and parcels are removed as grounding renders them insignificant in volume.When interpolating between any two key-frames, one starts with two lists of parcels,Lk and Lk+17 separated by the time At. Let P be an arbitrary parcel from Lk. If Pu isChapter 7. The System^ 89the identifier of P, then if a parcel 13 ' exists on Lk+1 such that .13.„ is its identifier, thenP and .13' match. The position, volume and radius of the interpolated parcel .P" can beinterpolated at time tkII where tk < tk" < tk+1 from:tk+1 —4+1 - tksPk + (1 — s)Pk+1V" = sVk + (1 — s )Vk+ ir = srk + (1 — s)rk+1If no point with the identifier Pu exists, then the parcel vanished in going from frame kto frame k + 1. This parcel will be interpolated at tk" from:tk-} i — tk' S =tk+1^tk13" = PkV = SVkr = srkThis process traverses both lists, removing matches from both lists. Parcels in Lk whichdo not match in Lk +i are also removed. When Lk is exhausted, parcels may remain onLk+1 . These parcels are the result of venting from the oil source which occurred betweentk and 4+1, and are not interpolated. This interpolation is demonstrated in figure Other Rendering TypesThe tidal model presented an opportunity to render in detail certain properties of thesurface currents and tidal elevation. These are specified with the particular renderingoptions in the script files.The tidal elevation over a domain can be rendered, using colour as the visual cue. Ina manner similar to the thickness map, samples are laid over the extent of the domainChapter 7. The System^ 90t=ot = t`t = 13Figure 7.7: Interpolation of an Oil Parcel Over Timewith each sample point corresponding to a pixel in the resulting image. A tidal elevationis sampled at each location. This is mapped into a colour by using the HSV colourmodel. The colour subspace for V = 1 corresponds to a circles with the hues along thecircumference. A specific colours is specified by its angle 7, with complementary colourslocated 180 degrees opposite of each other. The value of S, or saturation, is a ratio whichis 1 at the circumference of the circle and 0 at the center. Varying S has the effect ofinterpolating between the colour on the boundary and white, which is the colour for S=0.This colouring could also be achieved using an RGB colour space, but the transformationwould not be as straight-forward as HSV. If the tidal elevation is 77 then the colouring isbased upon:H = flifi7>0H = 180 — ,8 if ri < 0S =^Tiqmax51n the HSV model, V=1 corresponds to maximum intensity6 Alternately defined as hue7With red located at 0 degreesChapter 7. The System^ 9117 = 1where # is the current hue. Using a single hue, or the more familiar grey scale, torepresent a thickness or range is a relatively intuitive approach. Using the HSV systemallows the positive and negative range of the tidal elevation to be concurrently displayed.If alternately a grey scale was used, with full intensity representing the maximum heightof the tide, and black or no intensity to represent the minimum tidal elevation, then mid-tide or zero tidal elevation would be represented by the middle grey value. This cannotbe readily discerned. A function to map HSV to RGB completes the colouring. Theeffect of this system is to set the opposite tide amplitudes to complementary colours.As the cycle from high to low occurs, the colour will be white when the elevation isapproximately mean tide. In an animation, this shows the rise and fall of the tide if anadequate time step and number of frames is used. However, as the change in elevation isvery smooth over time, no local spatial variation in elevation on a frame by frame basisis discernable.While little variation in elevation is evident, this is not true with the horizontalcurrents (figures 7.8 and 7.9). For the given domain a sample of the horizontal flow ateach point corresponding to a pixels is determined. The maximum flow in the region mustbe determined or estimated for the region. At present a grey scale value is determinedat a location based upon:pp = (28 1) Tfmax — Tf Tim=The resulting image displays areas of large current flow as bright white, with areas ofrelative calm as a dark shade of grey, or black in the case of no flow at all. In the contextof an animation, the effects of channeling and shallow water can be seen clearly. Thiswill be expanded further in chapter 8.The ability to output a representation of an oil spill into the Optik ray tracing systemChapter 7. The System^ 92Figure 7.8: Current Flow in Howe^Figure 7.9: Current Flow in VancouverSound Harbour[Ama90] was provided. For every parcel in a slick, a corresponding Optik object iscreated. At present these objects are either particles or spheres. The object is translatedto its model coordinate position, and the radius of the object8 is set to the parcel'sradius (assuming the parcel is a sphere). An oil spill is represented by a collection ofintersecting spheres or disks. A sphere can be scaled into an ellipsoid with its two majoraxis described by the radius, and the minor axis representing the thickness of the parcel.Further, a collection of slick scripts can be composed into a single Optik script file, withz used to differentiate between time slices.7.4.4 Texture Mapping an Oil SpillOne of the intended purposes of generating the thickness map in the form of an image wasto use it as a texture map in a pre-existing rendering system, such as Optik. The slickat a given point will produce a colour dependent on its thickness, state of degradation,8With respect to a particle or sphere.Chapter 7. The System^ 93and the orientation of the surface at that point. Ideally for a given spill scenario, animage of the scene should be rendered which shows both the position and the generaloptical properties of the oil slick. Texture mapping is the geometric mapping of a twodimensional texture onto a three dimensional surface. In the context of this discussion,the texture is a two dimensional representation of the thickness of an oil slick over aregion, and the surface is the ocean's surface and the local foreshore which correspondto the spatial location of the slick. In our case, the surface is a plane. However, given anapproximation of the surface of the ocean resulting from wave action, such as developedin [Fou86], the scene can be made more consistent with reality. The oil slick is not atexture in the strict sense. When the slick is mapped onto the surface, at each pointwhere oil is present we would like the corresponding surface colour in the scene to beaffected. However, at each point where there is no contamination, then the radiance of thecorresponding region should remain unchanged. Most graphic systems support texturemaps in some form. Unfortunately, two problems surface on closer examination. Theseproblems are the registration between the texture map and the underlying surface, andcompositing the thickness map onto the surface. In Vertigo, Wavefront and Alias, it wasnot possible to bring the three dimensional model into the system without performingscaling, rotation and translation operations. As a result, the series of transformationsthat the object underwent would have to be applied exactly to the oil slick texture map,if there is to be correspondence between the position of the texture map and the locationsgenerated in the oil slick modelling system. The Optik rendering system had no suchconstraint. This texture mapping will be discussed further in chapter 8.The various models and custodial functions which comprise the oil spill simulationwere introduced in this chapter. Integrating the components developed during modellingyields a system which provides a means for defining and visualizing a marine oil spill.Chapter 8Experiments and Examples8.1 Standard Slick FormulationFor the testing stage, a fixed set of spatial, temporal, and physical parameters for thesimulation are set as a fixed model, from which the performance of adjusted models willbe measured. The parameters for this fixed model are listed in table 8.1. The resultingoil spill is shown in figure 8.1.Parameter Value Commentst July 17, 1992, 6 am Date of Sampling for Interpolated Wind DataP Fixed at test time Spatial Location of Sourcel'h 5.0 Horizontal Dispersion Coefficienta 0.025 Wind Coupling CoefficientcrT 1.00 Tide CoefficientP. 0.88 Specific Gravity of Test PetroleumPw 1.0275 Specific Gravity of Ocean WaterAt Time Step Adaptive or FixedVo 106 Litersfr 10 Liters per secondTable 8.1: Default Slick Model8.2 Testing for Sensitivity to Space and TimeThe first collection of tests will determine the sensitivity of the simulation to initialconditions in space and time. This will be a good measure of the quality of the overallintegrated simulation. For a fixed set of slick parameters, executed at time t and position94Chapter 8. Experiments and Examples^ 95Figure 8.1: The Standard Oil SpillP, the simulation after some period of time TD will produce a result. Visually, it will bemeasured by the final distribution achieved by the parcels within the spatial domain. Thespread, size, orientation, number of components, and the amount of shoreline affectedcan all be measured qualitatively by visual inspection of the rendered scene. If theslick parameters are held fixed, and the position P altered by some small amount AP,or the time of execution adjusted by a time increment At, then the expected outcomeof the simulation for time TD should be similar if the system is not chaotic. If theresult is unacceptably different, then the simulation is exhibiting chaotic behavior, or asensitivity to initial spatial and temporal conditions. If the system does not have thischaotic property, then small changes in position and time will result in a small changesin the distribution of the particles over time and space.Chapter 8. Experiments and Examples^ 96The means of determining this similarity can be qualitative or quantitative. The visualinspection of the two results is the most direct route. An intervening human in most caseswill be able to determine if, in general, the simulation is behaving in a stable manner, andthat patterns occurring in one simulation occur in the other. Quantitatively, there arevarious means by which the stability of the system could be determined. The analysisof the system on a parcel by parcel basis would be difficult. The random processesinvolved in the dispersion and grounding processes remove the possibility of retainingsensible correspondence between parcels from different simulations, despite the fixed spillparameters. Rather, generalized statistics would need to be derived, such as the centerof mass of the respective parcel clusters, the difference in slick area, and the volumeremaining in the slick, and the volume and center of mass of grounded oil. Further,statistical analysis by least squares or principal components of the parcels of the twosimulations could be performed, and the results accepted or rejected based on comparisonof these statistics.Figure 8.2: Position 200m West^Figure 8.3: Position 200m EastChapter 8. Experiments and Examples^ 97A collection of simulations were run on the system and it was visually confirmed thatsmall changes in position and time did not result in vastly different parcel distributions.The position was varied from 1 meter to 1000 meters radially about the initial positionand the simulation executed for 120 minutes. A change in the order of 500 meters usuallyresulted in markedly different advection and spread patterns at the end of the simulation.This must however be balanced by the initial position selected, as the resulting flows varyspatially. The starting time was varied between 1 and 60 minutes, with changes over 25minutes usually producing vastly different results. This would be expected, as the tideis changing with time, as are the winds. Figures 8.2, 8.3, 8.4, and 8.5 show thedistribution patterns achieved by varying space and time.Figure 8.4: Time Set 15m Back Figure 8.5: Time Set 15m Ahead8.3 Testing for Spreading and AdvectionThe spreading model was tested by altering the dispersion coefficient -yh and alternatelyturning the physical spreading off and on. The dispersion coefficient was varied fromChapter 8. Experiments and Examples^ 980.5 to 25 m2s -1 . The larger the dispersion coefficient, the greater the diffusion of theparticles over the surface of the ocean. The results of these tests are shown in figures8.8 and 8.9.Figure 8.6: No Turbulent Diffusion Figure 8.7: No Physical SpreadingThe physical spreading was tested by executing an instantaneous spill. The rate offlow was set to a value sufficiently high as to allow for the entire volume to vent in onetime step, and the effect of drift was canceled by setting a and ay very small.The physical spread was found to be consistent in this case. The results can be seenin figures 8.10, 8.6, and 8.7, which show the spreading resulting from just physical, justturbulent, and finally both forms of spreading.The wind coupling coefficient a and the tide coupling parameter ay, were adjusted.The literature revealed that the value for a is in the range of 0.01 and 0.08. A number oftests were run with values ranging from 0.001 and 0.10. The influence of this parameteron the spread and advection patterns is very apparent, with the distance induced bydrift proportional to this value. The aT parameter, while having no physical basis, wasChapter 8. Experiments and Examples^ 99Figure 8.8: -yh = 1.0m 2s -1 Figure 8.9: yh = 25.0m 2 8 -1intended to allow for varying the magnitude of the tidal velocity field. Spread patternsexclusively due to the wind, or to the tide, could then be observed. The results of thisexperimentation is shown in figures 8.11, and 8.12.Visualizing only the effects of advection of the particles was achieved by switching offthe physical spread, and setting -yh to be very small. This would be analogous to spillinga collection of buoys rather than oil. It shows how the system can be adjusted to serveother purposes, such as tracking a drifting vessel. A result of this type is shown in figure8.13.It was originally stated that the time step used in the simulation was adaptive. How-ever, for the purposes of exact synchronization during the testing and experimentation,the time step was fixed at 60 seconds. The simulations were rerun using the adaptivetime step, and a large fixed time step of 300 seconds. The results of this are shown infigure 8.14 and 8.15. They should be compared with figure 8.1. These images showhow the time step greatly affects the result.Chapter 8. Experiments and Examples^ 100Figure 8.10: No Wind or Tide, Both Spreading Types On8.4 Testing the Oil SourceThe effect of the density of the oil was tested. The difference was measured by visualinspection of the thickness maps generated by the extended visualizer. Key-frames ofthe simulation were captured involving petroleum with different specific gravities. Theeffects of petroleum density are manifested with respect to the amount of radial spreadoccurring at a particle level, and to the overall physical spreading of the oil body. Figures8.16 and 8.17 show the effect of marked differences in the densities introduced to themodel.One feature that is apparent is that the scale of diffusion of the high density slick isfar too large (figure 8.16). The fact that many of the spread discs do not overlap showsChapter 8. Experiments and Examples^ 101Figure 8.11: a = 0.08 Figure 8.12: a j = 0.10that the turbulent diffusion process is overwhelming the spreading process far to earlyin the simulation. This can be corrected by making -yh depend on the density of thepetroleum.8.5 Interpolation Between FramesTwo methods were applied to interpolating between two key-frames of slick simulationdata. The first method involved direct interpolation between the two thickness maps.This involves generating a slice of the thickness map by threshholding on thickness. Thisresults in a high contrast binary image, which has an edge detector applied. All edgeswhich are found are extracted and transformed back into the models coordinate system.The results are shown in figures 8.18 and 8.19.The second method applies linear interpolation directly between the key frames. Theresult of interpolating a slick is shown in figure 8.20.Chapter 8. Experiments and Examples^ 102Figure 8.13: No Effects from Spreading or Diffusion8.6 Altering Parameters Based on TestsThe parameter l'it was determined by a function of both sea state and oil density. Thiswas discussed in chapter 6. The magnitude of this parameter depends on the turbulenceand the properties of the oil. -yh was parameterized over sea state and density by applying:7h = (S + 1.0)L'10where S is the sea state reported as a Beaufort number, Po is the density of the sourcepetroleum, and -yo is the horizontal dispersion coefficient estimated for an oil with thesame density as water on a flat calm sea under no influence from currents. While thisformulation is empirical, it is constructed to be consistent with the processes in action.The term -yo/po grows with decreasing density. The term S +1.0 applies a linear increaseChapter 8. Experiments and Examples^ 103Figure 8.14: Adaptive Time Step Figure 8.15: A t = 300 Secondsin diffusion based on the surface turbulence, as estimated by the sea state model. Thehorizontal diffusion grows most rapidly when a low density petroleum is introduced intoa heavy sea. Using the values for the standard slick set, the result of this experiment areshown in figure Texture Mapping a Slick into the 3 -Dimensional ModelA thickness map was converted into the Wavefront system to render a consistent imageof an oil spill in the test domain. A second texture map to approximate the surfaceproperties of the domain was developed. Gross features of the region, such as forest cover,areas of human activity, and shoreline features were generated by hand in a rectangularcolour map corresponding to the modelled region. Alternately, if satellite imagery wasavailable, it could be used as the texture map. The 3-dimensional model was ported intoWavefront, and then both texture maps were registered to this model. The image of theoil distribution was coloured using a brownish-green shade, modelled after descriptionsChapter 8. Experiments and Examples^ 104Figure 8.16: Density po = 0.75 Figure 8.17: Density po = 0.99of thick oil slick colours described in [Hor72, Kaw70[. For each pixel of the texture map,if there is no oil then the corresponding area of the ocean's surface will have no changein properties. The result of this work is shown in 8.22.Chapter 8. Experiments and Examples^ 105Figure 8.18: Slick Thickness Map^Figure 8.19: Resulting EdgesFigure 8.20: Results of InterpolationChapter 8. Experiments and ExamplesFigure 8.21: Slick Produced Using Adaptive Spread Algorithm106Chapter 8. Experiments and Examples^ 107Figure 8.22: An Oil Slick in Vancouver HarbourChapter 9Conclusion9.1 Future WorkA system of this ilk can rarely be termed complete. There is always one more devel-opment, or level of refinement which could be implemented. This chapter will brieflydiscuss the future of systems of this nature, and the short and long term changes thatcould be made to improve this particular system.This system operates strictly on the ocean's surface. While slick pollution is initiallya surface problem, the chemical degradation extends the pollution into the atmosphereand down into the water column A system of this nature must eventually address thisthird dimension, not only with respect to the oil source, but also in relation to theenvironmental processes. Subsurface currents induced by wind and tides, the surfacewind current, and the topography of the ocean surface, could be modelled, increasing theaccuracy of the system. This extension would not only increase the consistency of themodel, but would change the nature of the visualization that could be achieved. Volumerendering of the resulting subsurface and atmospheric plumes could be considered.The model would be improved by a more complete description and modelling ofthe weather process in action. This would require the collection of extensive data ona number of major petroleum crudes and distillates, and their chemical and physicalproperties. The inclusion of detailed weathering models would degrade the performanceof the simulation, and could potentially render it non-interactive. It should also be108Chapter 9. Conclusion^ 109noted that weather systems and fluid dynamics remain a open area of research. Eventhe leading edge of this science would be insufficient to completely describe all weatherrelated processes.The greatest weakness of this system is its lack of concrete validation. This is primarilya result of the lack of hard data on oil spills in the domain. While spill occurrence datais available, it does not provide any information on the spatial distribution of the slickover time, nor the breakup time of the spill. A possibility for validating the system wouldbe to acquire the required data for modelling in a domain which does have detailed spilldata from various events. Generation of the coastal terrain, tides and winds of the areaat the time of the actual spill could be drawn into the system using the grid mechanism.The parameters of the models could then be adjusted to produce results similar to theactual slick event.In terms of software, the system should be redeveloped to use an system independentinterface. This could be done using the X windows paradigm [Nye90], or by developinga simple and portable interface. The component software could use a re-evaluation interms of both general functionality, and time efficiency. With a careful redesign, the sys-tem could be redeveloped to run on most existing platforms (assuming adequate graphicscapabilities). Further, it could be extended into a general library of particle functionsand visualization. The generalization of the system into a generic particle system wouldrequire the addition of the third dimension in terms of both advective fields and topogra-phy. Vertical processes would then require more attention. With a system of this naturevisualization of the flows under the ocean's surface and into the atmosphere could berealized. This general system could be applied to any user defined environment.As a final extension to the system, the modelling of oil spill fighting methods shouldbe developed. Barriers, such as oil spill booms and skimmers could be introduced toblock the flow of the spill. Chemical dispersants, which increase the diffusion rate ofChapter 9. Conclusion^ 110the spill, could be included by altering the existing chemical degradation models. Theaddition of these features would allow for comparison of the amount of damage causedby an oil spill with or without contemporary oil spill fighting equipment. This could beused as a measure of the relative effect of these techniques. Further, theoretical spillfighting methods could be tested within the framework of this system.9.2 ConclusionThis thesis set out to develop a consistent system for modelling and visualizing thedynamic process of an oil spill in the marine environment. The criteria for developmentwas to produce a system which was interactive, both in terms of human interaction, andin the performance of the system, and possessed a generic mechanism for the input ofsource data, thus achieving domain independence.The system has achieved consistency at the level dictated by the assumptions andgeneralizations under which the models were developed. As has been stressed in previousdiscussions, the issue of modelling both the oil spill and the environment requires somegeneralizations. Attempts have been made to include all features that directly affect theprocess of oil spill pollution. While the inherent complexity of a system is not a validreason for choosing a model, or for not including the process, many of the componentsof the system have been generalized to a point where further condensation would beginto remove validity.The visualization developed provide a base level tool for determining the spatialdistributions of oil spill pollution over time. Additional functionality, in the form offlow visualization and in-depth rendering of particular environmental processes, allowsthe user of the system to analyze visually not only the resulting spill, but the majorinfluences on the spill.Chapter 9. Conclusion^ 111Over the last decade environmental awareness has spread and became a mainstreamtopic from the coffee shop to the board room. Whether or not this apparent concernwill lead to a more sensible, rational and forward looking approach to development andthe usage of ever-diminishing resources in the future is a question for philosophers andhistorians. Be that as it may, this system has attempted in a small way to contribute tothe body of knowledge that can be applied to arriving at this lofty goal.Appendix AAppendixA.0.1 Determining the Wind DriftTo determine the direction and magnitude of the wind-induced drift current, the solutionof the equations of motion with friction present must be solved. The system presentedis due to [Pon83]. Two terms must be considered: the wind as a result of the horizontalpressure gradient, and the vertical friction. Note that this derivation is intended for thenorthern hemisphere. Minor adjustments are required for the southern hemisphere. Thecomponents are solved separately and added together, yielding:f v = f (v. + v E )wheref = 2S2sinO is the Coriolis parameter, or planetary vorticity.ft = The angular speed of earths rotation about its axis.0 = the geographic latitude.v, = the wind velocityvE = the Ekman velocity associated with vertical shear frictionThis separation is possible as the equations are linear. The wind induced drift currentis determined by the superposition of these terms. Assumptions made include a homo-geneous water body, a flat surface, an infinitely deep ocean avoiding bottom friction,and infinite horizontal extent to avoid associated boundary problems. Given this, andassuming the wind to be blowing in the y direction the Ekman equations becomes:u = lio cos(i + r±E z)eff z112Appendix A. Appendix^ 113v = Vo sin(i + if.-z)e* zwhereVo is the total surface current or Ekman flow,DE is the Ekman depth of frictional influence.If we examine the system where z = 0, the solutions become:u = Vo cos(i)v = Vo sin(i)This implies that the wind induced surface current should flow at 45 degrees counter-clockwise to the wind direction. To obtain a numerical relationship between the surfacecurrent Vo , observational experiments determined:Tn = PaCD I W I, the wind stress magnitude,wherepa is the density of air,CD is a dimensionless constant,W I is the wind speed,The surface current and the wind speeds are related by:1411 _  0.0127 W VIC01)yieldingDE =  4.3W N/sin101)From this if the wind in ms -1 is known, Vo and DE can be found at any depth, andin particular on the surface. At the 45th parallel the ratio CF = iv,, sometimes calledthe wind coupling coefficient, has been estimated at 0.015.Appendix A. Appendix^ 114A.0.2 System SpecificsThe system was developed under ANSI C on various Silicon Graphics platforms. Thecode is split into three groups, user interfaces, graphics and simulation-utility. The userinterface was developed using the FORMS interfaces library. The graphics routines usesthe Silicon Graphics GL embedded programming language for graphics initialization andprimitives. The simulation-utility group was developed primarily by the author, withsome matrix and vector routines borrowed from the MxMatrix C library [Buc90].Running the system is a multi step process. The user must first set a search path tofind the oilsim executables, which are located in "/nfs/grads10/grads10/dvanb/binfsgi".Next, the user must extract a set of data upon which a simulation may be run. Thedatabase program does not require any support files, but will require that permission isset for read/write for the graphics group. To extract a dataset type database. This willbring up a graphic display of the domain, and an interactive panel. Specify the path andname of the dataset, and extract. When the extraction is complete the user can inspectthe resultant directory. It will contain files with the names g100b, tide.; tide.u, tide.vwhich are respectively the terrain and tidal representation. The directory wind containsdata files describing the wind.To run the oil simulator, the user must have a copy of the file parameters. dat. Thisfile contains the defaults of the physical and simulation parameters used by the system.The default parameters.dat contains:1. wind coupling coefficient, default 0.032. tide coupling coefficient, default 1.003. horizontal dispersion coefficient, default 1.04. vertical dispersion coefficient, default 0.001Appendix A. Appendix^ 1155. tidal component flag 1, default 16. tidal component flag 2, default 17. tidal component flag 3, default 18. tidal component flag 4, default 19. tidal component flag 5, default 110. tidal component flag 6, default 111. tidal component flag 7, default 112. tidal component flag 8, default 113. tidal component flag 9, default 114. tidal component flag 10, default 115. tidal component flag 11, default 116. tidal component flag 12, default 117. tidal component flag 13, default 118. tidal component flag 14, default 119. hue for HSV, default 75.020. upper range for elevation, default 10.021. lower range for elevation, default 0.022. denominator and numerator of sea state coefficient, default 10,1Appendix A. Appendix^ 11623. fixed drop point flag, default 024. fixed drop point, default (55000,33400)25. iterations before freezing, default 0 (Dont freeze)The user can sensibly alter these parameters to any value. To execute the systemtype oilsim. The input panels are self describing.Bibliography[A1h75] S. W. Ahlstrom, A Mathematical Model for Predicting the Transport of Oil SlicksIn Marine Waters, Battelle Pacific Northwest Laboratories, Richland, Washington,1975.[A1R89] A. H. Al-Rabeh, H. M. Cekirge, and N. Gunay, A Stochastic Simulation of OilSpill Fate and Transportation, Applied Mathematical Modelling, Vol. 13, pp 322-342,June 1989.[Ama90] J. Amanatides, A. Woo, J. Buchanan, P. Poulin, C. Houle, and L. Slipp, OptikUsers Manual, Version 2.1, 1990[ASP83] American Society of Photogrammetry, Manual of Remote Sensing, 2nd ed., FallsChurch, Va., 1983[Ara82] K. P. Aravamudan, P. Raj, J. Ostlund, E. Newman, and W. Tucker, Break upof Oil on the Sea - Simplified Models and Step by Step Calculations, Report No.CG-D-28-82, U. S. Coast Guard, NTIS No. ADA118 372, 1982.[Ayo87] K. Ayotte, The Development and Application of a One-Level Mesoscale Atmo-spheric Model over the Lower Fraser Valley of British Columbia, Masters Thesis,University of British Columbia, 1986[Be187] P. C. Bell, and R. M. O'Keefe, Visual Interactive Simulation - History, recentdevelopments, and major issues, Simulation, Vol. 49, Number 3, pp 109-116, 1987.[Ben77] A. P. Bentz, Effects of Weathering on Oil, Workshop on Pattern RecognitionApplied to Oil Identification, IEEE Cat. 76, 1977.[Bow83] K. F. Bowden, Physical Oceanography of Coastal Waters, Ellis Horwood Lim-ited, 1983.[Bra92] P. Bragatto, N. Mazzino, and P. Palamidese, Animated Visualization of ScalarFields, Second Eurographics Workshop on Animation and Simulation, Vienna, pp115-127, 1992.[Bre77] J. A. Brewer, and D. C. Anderson, Visual Interaction with Overhauser Curvesand Surfaces, Quarterly Report of Siggraph, Vol. 11, Number 2, pp 132-137, 1977.117Bibliography^ 118[Bri91] M. Briscolini, and P. Santangelo, Animation of computer simulations of two-dimensional turbulence and three-dimensional flows, IBM Journal of Research andDevelopment, Vol. 35, Number 1/2, pp 119-137, March 1991.[Bro90] F. P. Brooks Jr., Ming Ouh-Young and James J. Batter, Project GROPE - HapticDisplays for Scientific Visualizations, Proceedings ACM Siggraph, pp 177-185, 1990.[Buc90] J. Buchanan., MxMatrix Math Routines, Imager Laboratory, University ofBritish Columbia, 1990.[Cek88] H. M. Cekirge, A. H. Al-Rabeh, and N. Gunay, Determining the Wind DrivenSurface Currents for Prediction of Movements of Oil Slicks in the Arabian Gulf,Computers Mathematical Applications, Vol. 17, Number 11, pp 1449-1453, 1988.[Cra86] M. A. Crane, Introduction to Ocean Modelling, Advance Physical OceanographicNumerical Modelling, D. Riedel Publishing Company, 1986.[Cre88] P. B. Crean, T. S. Murty, and J. A. Stronach, Mathematical Modelling of Tidesand Estuarine Circulation: The Coastal Seas of Southern British Columbia andWashington State, Springer-Verlag, Berlin, 1988.[EMR88] Digital Topographic Map Data, Topographical Mapping Division, Departmentof Energy, Mines and Resources, Canada Center for Mapping, Ottawa, Canada, 1988[Eli77a] A. Eliassen, and K. Pedersen, Meteorology, an Introductory Course, Volume I,Physical Processes and Motion, Universitetsforlaget, Oslo, Norway, 1977.[Eli77b] A. Eliassen, and K. Pedersen, Meteorology, an Introductory Course, Volume II,Application to Weather and Weather Systems, Universitetsforlaget, Oslo, Norway,1977.[E1186] A. J. Elliott, N. Hurford, and J. C. Penn, Shear Diffusion and the Spreading ofOil, Marine Pollution Bulletin, Vol. 17, Number 7, pp 308-313, 1986.[Enc87] Encyclopedia of Chemical Technology, Third Ed., Vol 17, John Wiley and Sons,New York, 1987.[Far87] E.J. Farrel, Visual Interpretation of Complex Data, I.B.M. Systems Journal, Vol.26, No. 2, pp 174-199, 1987.[Fay71] J. A. Fay, Physical Processes in the spread of oil on a water surface, Proceedingsof the Joint Conferences on the Prevention and Control of Oil Spills, AmericanPetroleum Institute, Washington, D.C., pp 463-467, 1971.Bibliography^ 119[Fo182] J. D. Foley, and A. van Dam, Fundamentals of Interactive Computer Graphics,Addison Wesley, Boston, 1982.[Fou86] Alain Fournier and W.T. Reeves, Aug 86, A Simple Model of Ocean Waves,SIGGRAPH 86, Vol 20, Num 4, pp 75-84, 1986.[Fou89] A. Fournier, The Modelling of Natural Phenomena, Graphics Interface 89, pp191-199, 1989.[GL90] Graphics Library Programming Guide, Silicon Graphics Inc., Document Number007-1210-010, 1989[Gun87] E.R. Gundlach, Oil-Holding Capacities and Removal Coefficients for DifferentShoreline Types, Proceedings, Oil Spill Conference, A.P.I, U.S.C.G and E.P.A, 1987[Hab90] R. B. Haber, and D. A. McNabb, Visualization Idioms: A Conceptual Modelfor Scientific Visualization Systems, IEEE Computer Society Press Tutorial, LosAlamitos, pp 74-92, 1990.[Hec86] P. S. Heckbert, Survey of Texture Mapping, Proceedings Graphics Interface, pp56-67, 1986.[Hib89] W. Hibbard, and D. Santek, Visualizing Large Data Sets in the Earth Sciences,Computer, pp 53-57, August 1989.[Ho179] J. R. Holton, An Introduction to Dynamic Meteorology, Second Edition, Aca-demic Press, San Diego, 1979.[Hor72] B. Horstein, The Appearance and Visibility of Thin Oil Films on Water, EPA-R2-72-039, August 1972.[Hua83] J. C. Huang, Review of the State of the Art of Oil Spill Fate-Behavior Models.,Proceedings: 1983 Oil Spill Conference, American Petroluem Institute, pp 313-322,1983.[Hun90] Hung Tao Shen, Poojitha D. Yapa, De Sheng Wang, and Xiao Qing Yang, AMathematical Model for Transport and Mixing of Oil Slicks in Rivers, Proceedingsof National Conference, American Society of Civil Engineers, Vol. 1, pp 75-81, 1990.[Jor80] R. E. Jordan, and J. R. Payne, Fate and Weathering of Petroleum Spills in theMarine Environment, Ann Arbor Science, 1980.[Kaw70] Kawahara and Ballinger, Characterization of Oil Slicks on Surface Waters, In-dustrial Engineer and Chemical Products Research Development, Vol. 9, Number 4,pp 553-558, 1970.Bibliography^ 120[Lar88] R. W. Lardner, W. L. Lehr, R. J. Fraga, and M. A. Sarhan, A model of residualcurrents and pollutant transport in the Arabian Gulf, Applied Mathematical Mod-elling, Vol. 12, pp 379-390, 1988.[Mac77] D. Mackay, and P. J. Leinonen, Mathematical Models of the Behavior of OilSpills on Water with Natural and Chemical Dispersion, Economic and TechnicalRevue Report, University of Toronto, 1977.[Mac78] Mackay and Wilks, A Statistical Analysis of Oil Spills in Canada, ResearchReport No. 50, Government of Canada. June, 1978.[Mac80] D. Mackay, S. Paterson, and K. Trudel, A Mathematical Model of Oil SpillBehavior, Environmental Protection Service, Fisheries and Environment Canada,1980.[Mas85] C. F. Mass, and D. P. Dempsey, A One-Level, Mesoscale Model for DiagnosingSurface Winds in Mountainous and Coastal Regions, Monthly Weather Review, July1985.[Nae79] A. Naess, The Mixing of Oil Spills Into the Sea by Breaking Waves, Journal ofPetroleum Technology, pp 1113-1122, 1980.[Nie89] G. M. Nielson, Visualization in Scientific Computing, Computer, pp 10-25, Au-gust 1989.[Nih84] C. J. Nihoul, A Non - Linear Mathematical Model for the Transport and Spread-ing of Oil Spills, Ecological Modelling, Number 22, pp 325-339, 1984.[Nye90] A. Nye, Volume 1 : Xlib Programming Manual for Version II, O'Reilly andAssociates Inc, Sebastopol, California, April 1990.[Obr85] J. O'Brien, Advanced Physical Oceanographic Numerical Models, NATO ASISeries, 1985.[Ove91] M.H. Overmars, Forms : A Graphical User Interface Toolkit for Silicon GraphicsWorkstations, Version 1.5, 1991[Pav82] T. Pavlidis, Graphics and Image Processing, Computer Science Press, Maryland,U.S.A, 1982.[Pay85] J. R. Payne, and C. R. Phillips, Petroleum Spills in the Marine Environment,Lewis Publishing Incorporated, 1985.[Pe186] E. Pelletier, Oil Spill in the St. Lawrence Estuary: A Preliminary Approach toa Risk Estimation Model, Proceedings International Symposium Natural and ManMade Hazards, D. Reidel Publishing Company, 1986.Bibliography^ 121[Pes91] R. L. Peskins, S. S. Walther, A. M. Froncioni, and T. I. Boubez, InteractiveQuantitative Visualization, I.B.M. Journal of Research and Development, Vol. 35,Number 1/2, pp 205-226, 1989.[Pie74] R. A. Pielke, A Three-Dimensional Numerical Model of the Sea Breezes OverSouth Florida, Monthly Weather Review, Vol. 102, February, pp 115-134, 1974.[Pon83] S. Pond, and G. L. Pickard, Introductory Dynamical Oceanography, PergamonPress, Toronto, 1983.[Pre85] F.P. Preparata, M.I. Shamos, Computational Geometry:An Introduction,Springer-Verlag New York, 1985.[Pur57] G. A. Purdy, Petroleum, Prehistoric to Petrochemical, Copp Clark PublishingCompany, Toronto, 1957.[Rap89] D. Rapaport, Visualizing Physics, Computers in Physics, Vol. 3, pp 19-29,September 1989.[S1o87] K. R. Sloan, Jr. and J. Painter, From Contours to Surfaces: Testbed and InitialResults, Graphics Interface 87, 1987.[Ste88] D. G. Steyn, and I. G. McKendry, Quantitative and Qualitative Evaluation ofa Three-Dimensional Mesoscale Numerical Model Simulation of a Sea Breeze inComplex Terrain, Monthly Weather Review, Vol. 116, Number 10, pp 1915-1926,October 1988.[Smi90] R. B. Smith, Why Can't Stably Stratified Air Rise over High Ground?, Atmo-spheric Processes over Complex Terrain, Meteorological Monographs AMS, Vol. 23,Number 45, June, 1990.[Spa88] M. L. Spaulding, A State-of-the-Art Review of Oil Spill Trajectory and FateModelling, Oil and Chemical Pollution, Number 4, pp 39-55, 1988.[Tip82] V.T. Tippie, and D. R. Kester, Impact of Marine Pollution on Society, J. F.Bergin Publishers, New York, 1982.[Tor82] V. Torre, and T. A. Poggio, On Edge Detection, IEEE Transactions on PatternMatching and Machine Intelligence, Vol 8, No 2, pp 147-163, 1986[Tuc80] M. J. Tucker, The Importance of Oceanography to the Offshore Oil Industry,PETROMAR 80, pp 497-511, 1980.[Tuf83] E. 11. Tufte, The Visual Display of Quantitative Information, Graphics Press,Cheshire, Connecticut, 1983.Bibliography^ 122[Ven87] S. Venkatesh, The Oil Spill Behavior Model of the Canadian Environment Ser-vice, Part 1: Theory and Model Evaluation, Atmosphere and Ocean, Vol. 22, pp93- 108, 1987.[Ven90] S. Venkatesh, Model Simulations of the Drift and Spread of the Exxon ValdezOil Spill, Atmosphere and Ocean, Vol. 28, pp 90 - 105, 1990.[Viz73] K. N. Vizy, Detecting and Monitoring Oil Slicks with Aerial Photographs, East-man Kodak Company, 1973.[Yap89] P. D. Yapa, R. J. Thomas Jr., R. S. Rutherford, and H. T. Shen, MicrocomputerModel for Oil Spill Simulation, Journal of Computing in Civil Engineering, Number3, pp 33-46, 1989.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items