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 VISUALIZING MARINE OIL SPILLS By David van Blankenstein B. Sc. (Computing Science) Simon Fraser University  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF SCIENCE  in THE FACULTY OF GRADUATE STUDIES COMPUTER SCIENCE  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA  December 1992 © David van Blankenstein, 1992  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  Department of  Co tr F si+ t se ,  ,  6  GI  The University of British Columbia Vancouver, Canada  Date^7;')  DE-6 (2/88)  --  I.^Q2.-  e  Abstract  The main goals of this thesis are to develop a system to simulate and visualize the dynamics of a marine oil spill. This necessitates the development of computational models of the properties of petroleum, the containing environment, and a melding of them into a consistent simulation of the fate and transport of oil spills. Techniques of computer graphics will be applied to present the state of this system as a visual output. Through a large set of accessible parameters and data, the variation of the system in terms of spill simulation, and model tuning, will be visualized over both time and space. The resulting system is a combination of research from a spectrum of scientific and technological areas. It provides a means of abstracting the problem of pollution through an accessible and consistent collection of mathematical models, and presents its results for interpretation through visualization. As a practical illustration of the system, it will be tested and verified in the region of southwestern British Columbia.  ii  Table of Contents  ii  Abstract^  iii  Table of Contents^  viii  List of Tables^ List of Figures^  ix  Acknowledgement^  xi  1  2  Introduction  1  1.1  Oil Spill Pollution ^  1  1.2  Visualization and Simulation ^  2  1.3  This Thesis ^  3  1.3.1^The Database ^  4  1.3.2^The Interactive Visualizer ^  4  1.3.3^The Extended Visualizer ^  4  1.3.4^Testing the System ^  5  Oil Spills  6  2.1  Oil Spills ^  6  2.2  The Properties of Petroleum ^  7  2.3  The Physics of an Oil Spill ^  9  2.3.1^Spreading ^  9  iii  3  2.3.2^Drift or Transport ^  12  2.3.3^Evaporation ^  13  2.3.4^Vertical Dispersion, Dissolution and Emulsions ^  13  2.3.5^Grounding and Sinking ^  14  2.3.6^Other Processes ^  15  2.4^Spatial and Time Scales of an Oil Spill ^  16  2.5^Visual Characteristics of Petroleum on Water ^  16  The Environment  19  3.1^The Environment ^  19  3.1.1^Study Area ^ 3.2^Current Inducing Processes ^  4  5  19 20  3.2.1^The Tides ^  21  3.2.2^The Winds ^  23  3.2.3^Other Oceanographic Flows ^  25  3.2.4^Waves and the Sea State ^  25  3.3^Processes Not Considered ^  26  Visualization  29  4.1  The Role of Visualization ^  29  4.2  Current Trends in Visualization ^  31  4.3  Interactive Simulation and Visualization ^  34  4.4 Visualization of Model Problem ^  35  Previous Work  37  5.1  Statistical Modelling ^  5.2 Transport Fate Models ^ iv  37 39  5.3 Visualization of Previous Systems ^ 6 Modelling^ 6.1 Modelling ^ 6.1.1 Restriction to the Ocean's Surface ^ 6.2 Modelling the Environment ^  42 44 44 44 45  6.2.1 The Terrain ^  45  6.2.2 The Land Query ^  46  6.2.3 Grid Data Description ^  48  6.2.4 The Tide ^  49  6.2.5 The Wind ^  52  6.2.6 Empirical Wind Modelling ^  53  6.2.7 Numerical Wind Modelling ^  54  6.2.8 Interpolation of Wind Field Data ^  55  6.2.9 Interpolation of Wind and Tide ^  56  6.2.10 Other Advecting Flows ^  58  6.3 Other Environmental Factors ^  58  6.4 Modelling of the Oil Spill Process ^  59  6.4.1 Oil Release ^  61  6.4.2 Transportation ^  61  6.4.3 Turbulent Diffusion and Spreading ^  62  6.4.4 Evaporation and Dissolution ^  66  6.4.5 Grounding ^  67  6.4.6 Weathering Processes not Modelled ^  68  6.5 Determining the Time Step ^  69  7  The System  7.1  Interfaces ^  71  7.2  The Database ^  72  7.2.1^Data Storage and Extraction ^  73  The Interactive Visualizer ^  74  7.3.1^Running a Simulation  75  7.3  7.4  8  9  71  ^  7.3.2^Displaying the Simulation ^  77  7.3.3^Sampling from a Simulation ^  79  7.3.4^Flow Visualization ^  80  The Extended Visualizer ^  83  7.4.1^Determining the Thickness Map ^  84  7.4.2^Interpolation of the Slick Between Frames  86  ^  7.4.3^Other Rendering Types ^  89  7.4.4^Texture Mapping an Oil Spill ^  92  Experiments and Examples  94  8.1  Standard Slick Formulation  8.2  Testing for Sensitivity to Space and Time ^  94  8.3  Testing for Spreading and Advection ^  97  8.4  Testing the Oil Source ^  100  8.5  Interpolation Between Frames ^  101  8.6  Altering Parameters Based on Tests ^  102  8.7  Texture Mapping a Slick into the 3-Dimensional Model ^  103  ^  Conclusion  94  108  9.1  Future Work ^  108  9.2  Conclusion ^  110 vi  Appendices^  111  A Appendix^  112  A.0.1 Determining the Wind Drift ^  112  A.0.2 System Specifics ^  114 117  Bibliography^  vii  List of Tables  2.1 Various Oils and Their Properties ^  9  3.1 Semi-Diurnal Components of the Tide ^  22  8.1 Default Slick Model ^  94  viii  List of Figures  2.1  Schematic of Processes at Work on Oil Slick ^  10  2.2  Time Scale of the Various Processes at Work on a Large Oil Spill ^. . .  17  3.1  Tidal Forces Resulting from the Moon  20  3.2  Vancouver and Surrounds ^  28  6.1  Bilinear Interpolating of Tide Flow ^  57  6.2  Spreading of Model Slick as Concentric Ellipses ^  65  6.3  Particle-Shoreline Avoidance  68  7.1  A Model Slick as Viewed in the Visualizer ^  78  7.2  Tidal Flow, Time 0900:17-Jul-92 ^  81  7.3  Drift Method, Time 0900:17-Jul-92 ^  81  7.4  Wind Flow, Time 0900:17-Jul-92  ^  82  7.5  Viewpoint and Sampling Region  ^  86  7.6  Sampling of p x q Grid  7.7  Interpolation of an Oil Parcel Over Time  7.8  Current Flow in Howe Sound  7.9  Current Flow in Vancouver Harbour  8.1  The Standard Oil Spill ^  95  8.2  Position 200m West  ^  96  8.3  Position 200m East  ^  96  8.4  Time Set 15m Back  ^  97  ^  ^  ^ ^  ^  ix  ^  86 90 92 92  8.5 Time Set 15m Ahead ^  97  8.6 No Turbulent Diffusion ^  98  8.7 No Physical Spreading ^  98  = 1.0m 2 s -1  ^  99  8.9 -yh = 25.0m 2 .9 -1  ^  99  8.8  7h  8.10 No Wind or Tide, Both Spreading Types On ^  100  8.11 a =0.08 ^  101  8.12 at = 0.10 ^  101  8.13 No Effects from Spreading or Diffusion ^  102  8.14 Adaptive Time Step ^  103  8.15 A i = 300 Seconds ^  103  8.16 Density p o = 0.75 ^  104  8.17 Density p o = 0.99 ^  104  8.18 Slick Thickness Map ^  105  8.19 Resulting Edges ^  105  8.20 Results of Interpolation ^  105  8.21 Slick Produced Using Adaptive Spread Algorithm ^ 106 8.22 An Oil Slick in Vancouver Harbour ^  107  Acknowledgement  The work presented here would not have been possible without the help from a number of people scattered across four Departments and two Universities. I would like to thank my supervisor Alain Fournier, for his frank and helpful discussions and guidance. Thanks too, to my second reader Brian Klinkenberg of the Department of Geography. Appreciation is extended to Guido Marionne of the Department of Oceanography for his help in acquiring the tidal model and for assisting me in many matters revolving around oceanography. Thanks also to Douw G. Steyn in the Department of Atmospheric Sciences, and to Keith Ayotte of York University for providing meteorological advice and the wind model. In the Department of Computer Science I would like to thank the many characters who lurk about the graphics labs, in particular John Buchanan, Peter Cahoon, Chris Healey, Pierre Poulin, Chris Romanzin, Rob Scharein, and Marcus Tessmann. This goes double 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 her support and encouragement over my long years of schooling.  xi  Chapter 1 Introduction  1.1 Oil Spill Pollution All industrialized societies use vast quantities of petroleum for fuel and materials. As the point of origin and the point of consumption are rarely the same, petroleum must be transported both for refining, and for distribution to consumers. This thesis will be concerned with accidents that occur during the transport and use of petroleum in a marine environment, and in particular with the pollution which results, termed a surface oil spill. The frequency of movement, and the great distances covered by shipping, results in the inevitable occurrence of marine accidents. Further, land-based human activities along coastal waterways also contribute to the chronic problem of oil spill pollution. This thesis will model these seaborne events. Massive spills such as with the Exxon Valdez, where fifty million liters of Prudhoe Bay crude was spilled into Prince William Sound, are a major stress on the regional ecosystems, over both time and space. Seemingly innocuous events, such as a few liters of diesel fuel spilling over a vessel's side are far less devastating, of short duration, and yet are far more common. In general, the frequency of oil spill events are inversely related to the volume of the spill. In the end both contribute to the long and short term degradation of the local marine environment. An oil spill is a complex event over spatial and temporal ranges. The circumstance of the spill (in terms of its proximity to the local coastline, size and composition), and the  1  Chapter 1. Introduction^  2  state of the marine environment (described by the winds, wave, tide and other natural processes), will dictate where the slick will transport, spread and eventually disperse into the water column The consistent modelling of this transportation and fate was a major goal of this thesis. 1.2  Visualization and Simulation  Computer graphics has fostered the growth of visualization as a mainstream method of interpreting the results of experimentation. With the advent of the digital computer a tool was provided which could both analyze, transform and display data in large quantities. Over the years, as processors and display devices have increased in power and numbers, and the techniques and algorithms of computer graphics have improved and multiplied, analysis of data by visual means has become more and more attractive. This thesis draws from this body of knowledge, and models, simulates, and visually presents the problem of oil spill pollution. The system is designed for consistency where possible, but due to the intrinsic complexity of the subject area many broad assumptions are made. The domain is restricted to be a small subset of the earth's surface, where approximations are made of the major oceanographic, meteorological and chemical processes which directly effect the petroleum slick's transport and degradation. Within this modelled region, and in context to the assumptions made, estimates of the fate of spills in this region, and their potential effect on the marine biota, can be realized. The methods involve a hybrid combination of both numerical, analytical and empirical models, and empirical and estimated data. From the resulting scheme, an inference towards a course of action can be made from interpretation of the potential damage caused in a simulated environment by a given case of an oil spill. This is judged primarily  Chapter 1. Introduction^  3  by visual means, through representations of the oil spill and its interactions with its environment. The primary flavours of visualization are in the animation of the transport, degradation and breakup of the spill, representations of the flows and processes that form the marine environment, and the rendering of affected areas. These representations range from simplistic images of component models, to realistic animations of specific spill scenario. This allows not only for the inference of the potential damages as the result of a given oil spill, but can also present the problem in a setting which is familiar to people in general. To view an animation of an oil spill spoiling a local beach can bring home the possible result of petroleum use and misuse, without actually causing damage to the environment. Variation of the many parameters of the component models, and visualization of the results have been facilitated. Viewing the model and the spill event not only over time but over the variation of certain parameters provides information on both the oil spill process, and on the merit of the component models. This allows the user of the system to adjust certain parameters of the component models for their specific needs. 1.3 This Thesis This thesis investigates the problem of oil spill pollution, and how it interacts and is affected by the environment it is contained in, by providing a visual representation for the 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 trends in visualization of scientific and industrial phenomenon will be presented. Computational models of a subset of the driving meteorological, and oceanographic forces 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^  4  1.3.1 The Database A database of state data for the various models is developed. A graphical front end will allow the selection of subsets of the spatial domain, and the portions of the various models which 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 to provide consistent and reliable data sets for testing and composing oil spill scenarios. 1.3.2 The Interactive Visualizer Interactive simulation and visualization environment are constructed. Their task is to allow the user to initialize, execute and visualize oil spill scenarios. It also allows the user to save the state of these simulations over time, visualize the major components of the environmental models, and allow the straightforward variation of particular model parameters. Flow visualization, and the slick distribution over time, will be the primary focus of the visual portion of this system. Sampling parameters, selection of various spatial subsets, and altering of graphics settings is also supported. 1.3.3 The Extended Visualizer The final environment allows for extended sampling of the system and the post processing of simulation data. Some forms of visualization are too computationally intensive to include in an interactive system. The visualization tasks allocated to the interactive system involve the display of model output readily extracted or generated from the simulation in action. Transformation of simulation output into more complex visual representations is the task of the extended visualizer. Other duties of this subsystem include the conversion of simulation data to other (graphic, visualization and animation) systems. The synchronization of time sampled output sets is formalized in a simple script language,  Chapter 1. Introduction^  5  which also provides the mechanism for generating output from this subsystem. 1.3.4 Testing the System A series of scenarios are generated from this system. These are intended to verify the models, and to allow for a broad range of views of the outputs. These visualization and graphic images range from the crude to the complex, and will cover not only oil spill simulations, but also the major environmental models. The thesis concludes with a discussion of the oil spill problem and the knowledge gained from modelling and visualization, and directions for this type of multidisciplinary application of computer graphics and computational modelling.  Chapter 2 Oil Spills  2.1 Oil Spills Marine oil spills are the intentional or accidental venting of some volume of petroleum into a coastal or oceanic system. It is a chronic problem which plagues all nations which rely 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 is a result of slight and perennial seepage of oil into the environment. Offshore drilling, onshore storage, pipelines, petroleum refineries and areas of intense marine activity usually cause this problem to some degree. This form of petroleum pollution usually manifests itself 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 contamination [Pe186, Tip82]. The more abrupt and certainly more recognizable form of oil pollution is the spilling of a raw or refined petroleum in its liquid state from a vessel or processing station. Spills of large 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 platform or grounded vessel, where the oil escapes from a rupture over time. In other cases, the vessel 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 of Prince William Sound in Alaska. Approximately 50,00,000 liters of Prudhoe Bay crude oil  6  Chapter 2. Oil Spills^  7  spilled into the sound, spoiling extensive stretches of the Alaskan coastline. Ecosystems, primarily the intertidal zone and various estuarine regions, were heavily contaminated by the slick, which travelled over 500 kilometers through the western Alaskan islands before finally 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 Amoco Cadiz spill (which occurred of the coast of northern France in March 16, 1978 spilling 180,000,000 liters of crude oil), resulted in similar short and long term environmental devastation of a large coastal region. The environmental effect of these spills is the disruption of the inter-tidal and estuarine biota. Toxic level contamination of the water column, and the total smothering of the intertidal zone are the most drastic environmental shock realized. Further, many human activities such as commercial fisheries, and tourism suffer as a result of these spills [Tip82]. These large spills tend to receive a large amount of attention, but are relatively infrequent. This statement is meant in no way to belittle the magnitude of the damage caused by large spills, but only that they are not the entire problem. A far greater volume of oil finds its way into the marine water system from smaller, more frequent spills [Tip82, Mac78]. Vessels flushing their bilges, careless refueling, and other small scale accidents introduce levels of petroleum contamination on a regular basis in heavily traveled coastal areas. While seemingly insignificant relative to 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 Petroleum The term oil refers to a variety of organic and inorganic liquids. Petroleum, which is latin for rock oil, is distinct from oil in general in that it is formed in rock formations over millions of years. For the remainder of this thesis the terms oil and petroleum will  Chapter 2. Oil Spills^  8  be 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 and shales). It is extracted and refined into various fuels, lubricants and petrochemicals. Petroleum is a colloidal dispersion of gas, liquid and solids, formed primarily of short and long chain hydrocarbon molecules, with impurities in the form of sulphur, oxygen, nitrogen and miscellaneous metals [Jor80, Pur57]. The fraction of the various long and short chain hydrocarbons largely determines the physical properties of a given petroleum compound. Short chain molecules tend to be light, with low boiling points, while the long chain molecules tend to be denser, with higher boiling points. Some important measures assigned to petroleum are the fractions of the various hydrocarbons and impurities, its density, viscosity and specific gravity. A given petroleum type or product can be described by providing the percentage or fraction that each distinct compound forms. This is usually stated as a mole fraction' The density and viscosity  2  describe the general property of the petroleum. The density is usually reported as a specific gravity 3 . Low viscosity oils, such as gasoline and kerosene have a viscosity less than 30 centistokes, with crudes and light bunker oils in the range of 30-2000, and high viscosity petroleum, such as heavy bunkers and weathered crudes, having values greater than 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 the component per molecules of the substance. 2 Viscosity of a solid or liquid is internal resistance to spread. 3 Specific gravity is the mass per unit volume of a compound over the mass per unit volume of water.  Chapter 2. Oil Spills^  Petroleum Type Arabian Heavy Crude Arabian Light Crude Iranian Heavy Crude Kerosene Light Residual Heavy Residual  9  Specific Gravity 0.886 0.831 0.966 0.787 0.924 0.984  Viscosity (Centistokes) 35.8 5.65 840 4.74 146 5400  Table 2.1: Various Oils and Their Properties 2.3 The Physics of an Oil Spill When petroleum is spilled into a large body of water, a complex series of competing process begin to act immediately. Occurring at both large and small scales, these processes lead to the eventual dispersal of the oil into the water body, atmosphere and the incident shorelines. They are driven by the composition of the oil, and the state of the marine environment, as illustrated in figure 2.1. 2.3.1 Spreading The most important process which an oil spill undergoes in its early stages is dispersion under the influence of turbulent and physical spreading. This is the lateral dilation of the spill under the influence of number of competing forces. This is one of the most important process which an oil spill undergoes, as it will largely determine the area of contamination, and the rate of the various weathering processes. Initially oil spreads under the force of gravity. This force is the dominant spreading force for up to the first week of an oil spill's life [Pay85, Jor80, Hua83]. This force is related to the thickness of the oil film, and the difference in density of the oil and the water. Gravitational spreading is counterbalanced in the first hours by inertia, which is rapidly overwhelmed by friction. The frictional term becomes the dominant force in  Chapter 2. Oil Spills  10  Figure 2.1: Schematic of Processes at Work on Oil Slick restricting 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 diminishing film thickness. The rate of this spreading is affected by the thickness gradient across the slick body. Differential thickness results in different spreading rates across the oil. Spreading under the influence of surface tension dominates during the later stages of an oil spill. This tends to occur only with slicks that last more than seven days. The differences between air/water and air/oil surface tensions are the prime cause of this force, which is independent of slick thickness. It tends to affect only those oils which  Chapter 2. Oil Spills^  11  have 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 spreading is uniform and radial about the source. Dependent on the density and viscosity of the oil, a minimum thickness will be reached, and spreading will stop. However, this is an extremely ideal case which would never occur outside of a laboratory testing tank. In the presence of wind and surface currents the body of the slick will distort. The velocity gradients 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 to grow at a far greater rate than the short axis, and is oriented in the direction of the surface winds [A1R89]. Further, most of the volume tends to accumulate near the leading edge (with respect to the local wind currents) of the slick. This process is called shear diffusion [E1186, Spa88], and is strongly influenced at this scale by high frequency, short term wind patterns and turbulent surface currents. The oil's viscosity and density, and the ambient and water temperatures also are involved in determining the rate of the spreading. A decrease in viscosity and density results in a decrease in spreading. An increase in temperature, either from the atmosphere 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 heavier petroleum compounds, such as unrefined crudes and residual oils, as their constituents are less susceptible due to their higher boiling points. However, the effects are readily observed in lighter oils and refined fuels, as their constituent are affected by temperature changes present in the environment, and undergo greater changes in composition over time [Mac77, Hua83, Spa88].  Chapter 2. Oil Spills^  12  2.3.2 Drift or Transport The drift is the large scale movement of the center of mass of the slick body. It is primarily driven by large scale forcing from wind, wave trains and surface currents. This process 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, long period wind currents. Wind shear stresses at the air/water interface induce a surface current which will cause the slick body to drift. The magnitude of this drift has been estimated in relation to the magnitude of the wind, and has been found to be in the range of one to five percent of the speed of the wind. Further, a drift angle of up to 45 degrees counterclockwise of the wind direction' is associated with this current (This is expanded in the Appendix). These drift currents extend into the water column, resulting in subsurface currents. The transport of the slick body is influenced by local tidal, geostrophic and wave induced 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 transport. However in regions where shallows, islands, reefs and other forms of topography cause blocking and deflection of these currents, the effect on the transport is far more pronounced [Bow83, Cre88]. The net advecting current is the superposition of these flow fields, although some research suggests that this simple linear sum may not be correct, due to non-linear effects in the system [A1h75, E1186]. The transport of the slick by this current is the primary cause of shoreline contamination, one of the most serious consequences of the slick's entry into the marine environment. 4  1n the northern hemisphere.  Chapter 2. Oil Spills^  13  2.3.3 Evaporation Oil is composed of various hydrocarbon compounds, which when considered as separate components, have different properties such as boiling points, specific gravity, vapour pressure and densities. The differences in composition of these components results in differential evaporation rates throughout the petroleum body. This thesis is primarily concerned with the lighter fractions of an oil body. Evaporation from light crudes and refined fuels' can result in the removal of as much as 75 percent of the slick mass. With heavier 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 viscosity and 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 Emulsions Dispersion and dissolution occur when oil from the spill enters the water column in a dissolved or dispersed form. In the dissolution process the lighter components are lost to the water column in particle form, which is a result of surface agitation and wind speed. The mass lost to this process tends to be small, as it is competing with an higher evaporation rate for the same fraction of the oil body. The amount of loss due to dissolution has been reported to 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^  14  oil body under the surface in the form of oil droplets. The majority of these droplets rise and rejoin the main body, while the smaller particles diffuse downwards into the water column. This process is extremely important, in that it will determine the long term lifetime of the oil spill. On a calm sea, this process is very weak, and the slick body, after the 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 mechanical dispersion, it can persist, increasing the possibility of the slick grounding. In a rough sea, a large amount of the slick volume can be dispersed into the water column. This process is a function of oil slick thickness, the interfacial tension between the oil and the 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 then disperse under the influence of evaporation and mechanical dispersion. In other cases oil-in-water and water-in-oil emulsions will form as a result of mechanical actions from surface turbulence and wave actions. The resulting compounds do not appear to disperse readily. 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 spreading is retarded, and a stable emulsion will form [Nae79]. The emulsion will take much longer to dissipate, and will require an applied force in the form of high energy wave actions to disperse the slick body. Further, the formation of these emulsions can increase the effective volume of the spill up to five fold [Hua83].  2.3.5 Grounding and Sinking Sinking occurs when the oil slick has achieved a specific gravity greater than that of the surrounding water body. This results partially from weathering processes, but also from the binding of sediments to the oil body.  Chapter 2. Oil Spills^  15  Grounding is caused by the deposition of the oil spill onto the inter-tidal zone. The rate and amount of oil which will collect on the shoreline is dependent on the ebbing and flow 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 important features associated with a shoreline are its slope, and its surface characteristics. Flat or low sloping areas tend tend to allow the oil to collect upon them. Tidal flats, marsh areas and gently sloping beaches have this property. If the grade is is more severe, then the oil will not as readily adhere to it. The surface properties of the shoreline can be described in terms of its roughness and the porous nature of the shoreline. Rocky, gravelly, and sandy coastlines allow oil to adhere in differing amounts. Further, the mixing of oil and shoreline sediments will also determine the amount of oil which is deposited on a shoreline. 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 this topic can be found in [Gun87]. 2.3.6 Other Processes  Photoxidation of hydrocarbons, caused largely by solar radiation, results in the production of water soluble compounds from a slick body. The physical processes at work are poorly understood [Hua83, Spa88]. It has been observed to be a process which takes weeks to months to occur, and affects only a small portion of the slick volume. Biodegradation is the slow process which results from marine organisms decomposing the oil slick. It is also a poorly understood process, with most knowledge derived from laboratory 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^  16  2.4 Spatial and Time Scales of an Oil Spill The spatial extent of an oil spill and the range of its transport, and the time scale under which an oil spill will exist, is clearly determined by the size and composition of the oil source, and the state of the environment. See figure 2.2. This must be tempered in the case of a large-scale accidental spill where the oil source is introducing petroleum into the environment continually. Hence, the time that will elapse from the entry of the oil into the marine environment, up until the spill body has completely dispersed, is highly variable. 2.5 Visual Characteristics of Petroleum on Water The 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 properties of the underlying water. If we assume a calm sea, a high solar altitude, and no bottom reflectance, then in the visible spectrum the radiance of an oil body is slightly higher than the surrounding ocean. The radiant properties of clean sea water and oil contrast highly in ultraviolet and the infrared bands of the spectrum. The reflectance properties of oil and water are almost identical in the green and red bands. This of course must be tempered by the composition of the water body. Turbid waters which are carrying large sediment loads reflect more light than clean water, and this will affect the contrast between the oil and water. The solar altitude largely determines the contrast between oil and water. With a low solar altitude oil is all but indistinguishable from the surrounding water [Viz73] To be able to reliably distinguish between oil and water readily other properties must be observed. Oil has a dampening effect on the ocean's surface. Small scale surface waves, or  Chapter 2. Oil Spills^  17  Figure 2.2: Time Scale of the Various Processes at Work on a Large Oil Spill capillary waves are removed by the affect of the riding oil body. The range of capillary waves which will be dampened depends on the thickness of the oil slick. Large oil spills have been observed to dampen out even wind driven waves. It is this property which makes 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 oilin-water emulsions yields a compound which ranges from a muddy brown to an almost white appearance. In the event of emulsion formation, the contrasts between the oil and  Chapter 2. Oil Spills^  18  water 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. Any two rays of light with the same wavelength emanating from a distant light source will be in phase. The difference in the refractive indices of the oil, air and water results in a phase shift between rays that were parallel on entry. Destructive and constructive interference occurs. Destructive interference tends to remove certain wavelengths, while constructive 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 colourless film, 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 in oils of thickness greater than 2000 nm. In this chapter the general properties of a marine oil spill were discussed. This included the properties, sources, and effects of oil spill pollution, and the affect of the regional environment on the oil spill.  Chapter 3  The Environment  3.1 The Environment The rates governing the spread, degradation and eventual breakup of the oil spill are a result of the state of the marine environment, and the resulting proximities of the slick to land and other features. The mechanics and dynamics of the atmosphere and ocean are extremely complex and variable. Those processes which are of concern to the model will be identified. Further refinement will specify those activities that occur in the domain within which the model will be developed. The scale of the intended research will not include 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 Area The study area or model domain will consist of North America from the 49th to the 50th latitude, and from the 122nd to 124th longitude (Figure 3.2). This includes the western portion of Greater Vancouver, the mouth of the Fraser River, Howe Sound, the south eastern section of the Strait of Georgia, and the Northern Gulf islands. This is a region of complex topography, formed primarily by the uplift of the Coast mountain range, and the 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 mouth of a river; rather, it is a semi-enclosed coastal region with a free connection to an ocean within which the sea water is considerably diluted with fresh water  19  Chapter 3. The Environment^  20  d  Moon Earth  Figure 3.1: Tidal Forces Resulting from the Moon Human activity in the region, centered around Vancouver, a rapidly growing city of approximately 1.5 million people, results in a high marine traffic density in both commercial 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 water column in the Georgia basin. Many other smaller rivers and streams contribute to the freshwater flow into the system. The region has a semi-diurnal unequal  2  tide. Winds in  the region are dominated by two large atmospheric pressure systems, the North Pacific high and the Aleutian low [Cre88]. Winds from the south east occur in the eastern portion of the domain when the low dominates, with winds from the north west occurring in the western half. The high pressure system induces westward blowing winds along the western portion. In general, the low dominates in the winter, and the high in the summer. 3.2 Current Inducing Processes An oil body is transported and dispersed in an estuarine region by the combined effect of surface-driven wind currents, tidal currents, geostrophic and residual currents, and wind 2 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 minima  Chapter 3. The Environment^  21  induced waves. The processes at work in the atmosphere and ocean are extremely complex, occurring at planetary and molecular scales. The laws of thermodynamics, fluid dynamics and planetary geophysics describe these systems. This thesis cannot hope to completely describe these systems, and so the interested reader is directed to the works of [Pon83] for a general discussion of dynamic oceanography, and [Ho179] for similar work on atmospheric processes.  3.2.1 The Tides The 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 to forces imparted on the earth by the moon and sun. The rise and fall is merely the consequence of the inward and outward flow of water. This tide generating force was shown by Newton to arise from the law of gravitation. If a system comprised of the earth and moon is considered, neglecting the rotation of both bodies around the sun, and assuming they are fixed spheres, then the two bodies would move as a result of mutual gravitational attraction. The force imparted by the moon on the earth is F1 at point P on the surface, and is F2 at the center of the earth C. This is shown in diagram 3.1. The term F is the tidal generating force, which is the horizontal component of F3 = F1 - F2. It can determined by the system 3 ME a  F = –2 mm — ( ci ) 3 g sin 20 where ME and Mn, are the masses of the moon and the earth, a is the radius of the earth, d is the distance between their centers, and 0 is the latitude of the point P. The relative positions of the moon and sun with respect to the earth vary with the orbits and rotations of these three bodies, yielding a number of constituents, each described  Chapter 3. The Environment^  22  by 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 than 24 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 individual harmonic constituents. Some of these constituents are listed in table 3.1. Constituent Principal Lunar Principal Solar Large Lunar Elliptic Luni-Solar  Symbol M2 S2 N2 K2  Period (hours) 12.42 12.00 12.66 11.97  Relative Magnitude 100 47 19 13  Table 3.1: Semi-Diurnal Components of the Tide There are two major theories on the ocean's response to tidal forces, the equilibrium theory and the dynamic theory. In the equilibrium theory, proposed by Newton, it is assumed that the earth is covered by a sea of fixed depth and density. The ocean responds to the applied tidal force by converging on the applied force, lifting the surface of the ocean in one location, and depressing it in others. A new oceanic equilibrium is attained when the hydrostatic pressure forces from the sloping ocean surface balances the geopotential forces. This theory underestimates the magnitudes of the oceanic tides, and cannot account for phase lags in constituents. The dynamic theory, first derived by Laplace, also assumes a world girdled by a uniform ocean. However, the tidal forces are seen to generate waves with the same period as the constituent. The movement of these waves yields surface currents, while the intersection of the waves with coastal features results in the rising and falling of the sea level we perceive. The M2 tidal force generates a forced wave with a period of 12.42 hours, over the ocean surface of the model planet. In principal, this theory should provide a means for calculating the tides at any given point on the ocean's surface. In reality oceanic  Chapter 3. The Environment^  23  boundaries and subsurface topography complicate the constituent equations. These features affect both the phase and magnitude of the many forced waves. In constricted channels a tidal wave will respond by accelerating, a result of the changing width and depth it encounters. If the channel narrows and decreases its depth, then the wave will build up on itself, increasing its amplitude. This phenomenon is called funneling. Further, boundary effects can be found in the vicinity of subsurface features, and points of land, which cause the deflection of gross tidal currents. In bays, tidal resonance can cause abnormal fluctuations in the tidal elevation, which can be markedly different than the nearby open ocean. Further, tidal waves are strongly influenced by friction in shallow and constricted waterways. These topographic affects on the tides abound in the Vancouver region. Clearly the horizontal flow of the tide has a large effect on the transport and spreading of a slick in a marine coastal region. Further, the rise and fall of the sea level will alternately expose and protect the shoreline, allowing a slick to affect the entire intertidal zone. The effect of the tide on petroleum slicks is small on the open ocean. The oscillations of the tide over the diurnal cycle tend to result in a small net displacement, and the magnitude 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 the tidal current is more pronounced, resulting in complex surface currents. For further reading on tidal phenomenon see [Bow83, Cra86, 0br85]  3.2.2 The Winds The planet earth is surrounded by a layer of air called the atmosphere. Kept in place by gravity, it is a low density mixture of gases such as nitrogen, oxygen, argon, water vapour and particulate matter. Although the earth and atmosphere receive as much energy as  Chapter 3. The Environment^  24  they radiate, this heat balance is only true on a global scale. Locally there are variations in the exchange of heat, and these differences place the atmosphere into motion. This thesis will be concerned with the net processes that occur in the circulation of the lower atmosphere or troposphere, and ultimately in the horizontal surface winds. The primary interest is in diagnosing the winds in the study domain. Wind on a planetary scale is generated primarily by the pressure force and the Coriolis force. The balancing of these forces is termed the geostrophic flow. The pressure force 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 in heat 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. The geostrophic flow determines the global circulation patterns, and flows from the west in the 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 the order 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 unstable air, and wind patterns induced by friction, topography, and heat exchanges between the land and atmosphere, at what is called the planetary boundary layer. These mesoscale winds are characteristic of the region they influence. The primary wind flow that must be investigated is induced by boundary effects caused when the geostrophic wind comes into contact with the terrain of a region. In the Pacific Northwest, the mountainous terrain greatly modifies the weather patterns in the lower troposphere [Mas85]. Channeling, blocking and deflection induce patterns which in no way resemble the synoptic patterns which drive them. These patterns are also induced 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 generated  Chapter 3. The Environment ^  25  primarily by vertical shear at the ground. The turbulent mixing which transfers heat and energy away, and momentum towards the surface during the diurnal cycle of heating and cooling, is mechanically, not thermally, driven. Further discussion of this topic will be deferred to chapter 6 where modelling of the mesoscale terrain induced wind currents will be covered. See [Ste88, Smi90, Eli77a] for further readings on the nature and properties of weather systems and in particular wind systems.  3.2.3 Other Oceanographic Flows Geostrophic or oceanic currents are primarily the result of the rotation of the earth about its axis, which is manifested in the form of the Coriolis force. This force results in the counterclockwise 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, as the semi-enclosed nature of this marine body provides a topographic buffer, and so these currents will not be considered.  3.2.4 Waves and the Sea State Waves in the marine environment range over a wide spectrum. The concern in this thesis is primarily with surface waves, which are the result of a sustained wind stress over a patch of the ocean's surface. Larger waves, such as internal and planetary waves will not be considered. Tidal waves are discussed in a previous section of this chapter. Surface waves have relatively short periods, ranging from 0.1 to 30 seconds, with wavelengths from one centimeter to hundreds of meters. The vertical accelerations of these waves can approach that of gravity with both horizontal and vertical accelerations overcoming geostrophic forces. These waves largely determine the sea conditions of a coastal region.  Chapter 3. The Environment ^  26  Surface waves can be found far from the area of influence which initially generated them in the form of swell. A swell will decay over time as it travels beyond its origins. If such a wave enters shallow water, then the bottom will effectively rejuvenate it, causing the swell to peak. The swell will eventually disperse, either within the oceanic body, or on the shoreline as surf. Wave analysis is a complex and fascinating field, but for the purposes of this research we are interested only in the general state of the sea surface, and not with individual waves. The sea state is a rather qualitative term used to describe the wave state of a given sea surface. Surfaces waves have both horizontal and vertical kinetic, and potential energy as a result of vertical displacement. A measure of the density of energy on the ocean's surface is given by E = pgH 2 /8 , where p is the density of the water, and H is the wave height 3 and g is the acceleration due to gravity at sea level. This energy will actively 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 in chapter 6.  3.3 Processes Not Considered There are many components of the environment not considered. The effect of fresh and salt water mixing is an important factor in an estuarine setting. At the mouth of a river there will be a distinct layering as the two liquids mix. The rise and fall of the tide helps in determining if the fresh water forms the top or bottom layer. The mixing is primarily vertical. Its effect on the net horizontal surface currents away from the mouth of the river is slight, and as such it will not be considered. The effect of spatial and temporal 3 Wave height, which is measured from peak to trough, is meant in the statistical sense. One methods of estimating wave height is to determine the significant wave height, or the mean height of the highest one third of the waves. In many cases this statistic is based purely on visual estimates  Chapter 3. The Environment^  27  variation of the ambient and ocean temperatures is important to processes such as large scale weather patterns, and the exchange of heat in the atmosphere. At the scale of the model problem they are not as important. The major processes which effect the movement and degradation of oil spills in a marine environment have been discussed. These processes include the terrain, wind and tidal currents, and are specific to the test region.  Chapter 3. The Environment^  Figure 3.2: Vancouver and Surrounds  28  Chapter 4 Visualization  Visualization in the context of computer graphics and science can be loosely defined as the transformation of data from some initial state into a visual representation via the processing of a digital system. A more frank definition is that visualization is merely the display of data [Tuf83]. From this representation, interpretation of the data and determination of the underlying processes at work is possible. The transformation can represent both the mathematical conversion of the data from one state into another through a model, or to the intermediate filtering and composing of the data into a visual form. This chapter will focus on the role of visualization or data display in scientific and industrial applications. A review of previous visualization research in related areas will be covered. The function of interaction in these systems will be emphasized. A discussion of the application of visualization to the model problem of oil spill pollution will complete the chapter.  4.1 The Role of Visualization Visualization 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 often can do more to describe the actual workings of the root process to a human than either the raw 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 becomes woefully inadequate. 29  Chapter 4. Visualization^  30  Digital 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 from both observation and mathematical modelling. On the other hand it is producing devices and systems which create or collect ever increasing volumes of data. A method must be applied to clearing the backlog of data produced by these systems. Thus scientific visualization is presented with a number of major tasks that it must fulfil. The first task that visualization encounters is the interpreting of large volumes of data. High bandwidth data sources such as satellites and radio telescopes, and medical imaging systems, to name but a few, produce voluminous quantities of data which must be summarized rapidly. Further, computational models of physical and chemical processes of ever increasing complexity are being developed. These models produce a boundless supply of data. The output from these sources must in some way be refined to allow for interpretation without removing the underlying patterns of the data. Complex and multidimensional data needs to be presented in a form which allows for the determination of trends and order within these sets. We exist in a multidimensional 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 many dimensions, 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 on a given data set in the form of graphic presentations. It must provide a base of functionality, which can be either general or application specific, and that can be applied to transforming a set of potentially multidimensional inputs into a visual output. This provides a window into local and global properties of the data set.  Chapter 4. Visualization^  31  4.2 Current Trends in Visualization Scientific visualization has branched out over the years to encompass many fields in the course of transforming data into visual image. The methods of computer graphics, vision, computer aided design, signal processing, and computer human interactions to name but a few, are applied to producing visualizations of data from a wide variety of applications. Researchers reiterate the problem associated with attempting to present too much data to a user [Hib89, Nie89, Be187, Hab90]. When this occurs, the user's ability to interpret the resulting image decreases and the effort required to generate the image increases. Each image should be derived to show no more than the effect of two or three parameters. If more information is shown, the usefulness of the visualization begins to degrade. While in many cases exotic graphics can be visually appealing, they quite often do not provide any assistance in understanding the processes that they were generated from. The application of animations to the visualization process has largely been realized through the rapid advancements of processing and display hardware. The frames of these animations are usually but not necessarily dependent on time. In the case of the model problem, an animation is formed of the spatial movements of the slick body over time. However, what variable the frames of the animation are based upon will vary from application to application. Researchers who emphasize the use of animation as a tool for visualization include [Rap89, Bri91, Be187, Bra92, Hib89, Nie89]. Many researchers discuss the linking of visualization to what can be termed as a computational experiment. Rather than deal with static data generated from an external model or observation, these systems have models coupled directly to the system. The state of the system is initialized, and the model executed, with the progress and results presented visually.  Chapter 4. Visualization^  32  A computational experiment begins with a loosely stated problem. For this thesis the statement is "What would the effect of a oil spill be on the local environment." The next step is to generalize and transform this statement into a well-posed mathematical problem, from which a reliable and adequate approximation can be developed. The simulation is executed and the output of the model is interpreted. The user can then adjust various quantities, both for the purposes of verifying models and for considering different setups in these models. The evaluation of these models by visual means greatly accelerates the process. An example of a visualization system developed to be run as a virtual lab was the system developed by Briscolini and Sanangelo in [Bri91]. A visualization system was developed to study the properties of two dimensional turbulence and three dimensional flows. The main form of visualization was the animation of flows over time, with the various physical parameters being altered over each set. The primary value measured by visual means was the vorticity. Rather than render lines of isovorticity, the effort was placed in determining the isovorticity of the surface. The blue and red components denoted a clockwise and counterclockwise flow respectively, with areas of low vorticity enhanced by the inclusion of white and green respectively. When extending their work to three dimensional flows, they reported no conclusive method could be found for visualizing these flows. Three dimensional vector representations of flows are discussed, and the 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 are used to represent low and high velocity flow. Intermediate flows grade from the blue into green and then red. They concluded that the two dimensional animations allowed for the observation and understanding of the mechanisms that control turbulent flows. Hibbard and Santek developed a visualization system specific to the data and problems associated with the earth sciences [Hib89]. The range of data sources included  Chapter 4. Visualization^  33  meteorological data, image data, and the results of numerical simulations. The emphasis of the system was on the processing of extremely large sets of data, some containing in the order of one billion observations. The numerical models, which were a prominent source 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 reliable mechanisms for bringing data into the system. The data, once in the system, could be converted into various data structures could be used in the visualizations. These structures included grids, images, paths, and structures of irregular data. Tools to interpolate these structures over time and space were included. The resulting system provided a general 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 and horizontal ticks labelled in the appropriate units for the given visualization. The boundary of the region in which the data was located is displayed on the bottom of the frame box. Three dimensional scalar fields were displayed using transparent and opaque contour surfaces on such scalar values as water vapour and potential temperature. Current flows were visualized using the trajectories of particles. The thickness and length of the line implies the speed and duration of the flow. Slices of three dimensional data, texture mapping of generated surfaces and other visual techniques were combined to provide the ability to readily construct animations of the given data. This was combined in the animation with a small rocking of the scene, to provide additional depth cues. Hibbard and Santek conclude by emphasizing the need for interactive systems. The users of the scientific community need both the versatility and speed afforded by an animation based system, but also need to be able to get up and running with the system with a minimum amount of training.  Chapter 4. Visualization^  34  4.3 Interactive Simulation and Visualization Visual interactive simulation, a term used by Bell and O'Keefe, describes a design methodology 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 context of discreet systems simulations, but the basic philosophy was just as valid for a continuous simulation. The notion was to afford the user the ability to interrupt and directly manipulate the model. This was done graphically though interactions with buttons, icons and other visual input mechanisms. Controlling parameters can be altered, objects removed from the simulation, and the system restarted from the point of interruption. The effect of these changes on the system will be apparent from the visual display of the model. This approach allows the user to avoid the problem of restarting the simulation each time a parameter was adjusted. These systems use animation extensively to show the state of the model. Icons which are readily identifiable as real world components of the model were used within these animations. Various statistics were presented in different windows, which changed with the model state. As the problem was shown in a familiar setting, the user was found to understood the changes occurring in the simulation, and could react by altering the parameters. Bell and O'keefe reported that users found this type of system very appealing. It provided a large amount of freedom for the user to shift attention between different components 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 verification and validation stage of a given model. In the long term, the authors suggest a number of points that can be followed in defining a methodology for such as system. They  Chapter 4. Visualization^  35  include getting the user involved early in the simulation, displaying the visualizations of the system as soon as possible, making interactions general, and transferring control of simulation directly to the end user. Peskins and associates discuss a system which was developed to allow the user to steer the development of a model through visual means in [Pes91]. This system is application independent, and was developed to allow a user to quantitatively measure properties of the output data, and present those results visually. The ability to adjust model parameters incrementally, judge the effect through visualizations, and thus steer the development of the model, is emphasized.  4.4 Visualization of Model Problem The work of previous researchers gives a gauge for the type of visualization that will be appropriate for the oil spill simulation system. The intended purpose of the system is to determine the spatial and temporal distribution of an oil slick over a given domain. To realize this, the models of the system must be tuned for these tasks, and the composing of 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 perform and react. In terms of visualization, the spatial distribution of the oil over the region will be presented in the form of an image showing the slick's position and extent over time within the tested domain. This will provide the user with a condensed visual representation of the state of the system for the given setup. In terms of oil spill pollution, this means that a researcher can vary position, time, volume, and other variables crucial to the process and guide the system towards various goals by interpretation of the system's  Chapter 4. Visualization^  36  visual output. If the need is to determine the areas affected by oil released from a given position, then altering the release time and volume will produce a series of images from which the effect can be estimated. These images will consist of the position of a thickness of oil over the test region's sea surface and shoreline for a given set of parameters and time. 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 been discussed. The current trends in this field with an emphasis on the physical sciences and human-machine interaction are covered. A number of developed systems and visualization methodologies have also been reviewed.  Chapter 5  Previous Work  Considerable research has been directed towards understanding and modelling the activities 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 impact analysis studies. While this thesis is primarily interested in the former, the latter helps to form the justification of applying technology to this problem, and gives a sense of the extent of the problem. 5.1 Statistical Modelling  Statistical models are concerned with predicting the long and short term likelihood of oil 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 of the number of accidents, their magnitude, and their potential affect on the surrounding environment. Mackay and Wilks provide a comprehensive analysis of the oil spill problem across Canada [Mac78]. The problem is approached in a sector by sector analysis of all petroleumbased activities and the resulting incidents of all types of spill pollution over an eight year period. The intent their model is to determine long term trends, and from this to estimate the distribution of oil spills over long periods of time. The sectors of interest to coastal marine regions are refineries, storage and ocean tanker traffic. Each sector is lA fate-transport model tracks the spatial position of the oil slick over time.  37  Chapter 5. Previous Work^  38  assigned 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 performing statistical analysis on these data, spill frequencies and distributions were determined. Conclusions on the spatial distributions of oil spills across Canada cannot be made as the data is aggregated for all regions. In addition the environmental conditions were not factored in, as the analysis was based on occurrence, rather than fate and affect. One of the conclusions of this analysis was that marine tanker transport is approximately twice as 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 any domain. The analysis considered local oceanographic, topographic, and meteorological effects, 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 contamination and the levels of toxic effect on the local biotas can be estimated. The effect of each feature, such as the local currents, proximity to shoreline and amount of spreading, were assigned a sensitivity value. For example, the probability of shoreline contamination was measured in terms of the spill's proximity and an energy description of the shoreline features. A shoreline feature is classified as being high energy if it has a low probability of 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 these statistics was then used for estimating open water, shoreline and direct biological impact. It was concluded that environmental risk was dependent upon the spatial location of the spill within the domain. This fact largely determined the spreading, and resultant transport by local currents.  Chapter 5. Previous Work^  39  5.2 Transport Fate Models Transport fate models are primarily concerned with the spatial distribution of the oil body 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 both the 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 and Semiamhoo Bay region, which is due south of this system's model domain, and that it developed reliable stochastic methods for estimating turbulent dispersion. Major environmental processes, such as the wind and tidal flows, were empirically modelled. Processes such as geostrophic surface currents, fresh water runoff and wave induced currents are discussed, but are not included in the model formulation. A stochastic method for modelling turbulent dispersion is designed. The slick body is viewed as a collection of independent particles, which are displaced by a Lagrangian random process. Further, the particles are independently advected. This method of modelling the slick body and the dispersion has become more or less a modelling standard in this area of research, used by researches such as [Ven87, A1R89, Lar88]. Physical spreading was modelled using the work of [Fay71]. The physical spreading is only considered in the first few minutes of the simulation, and the model depends on turbulent dispersion for the duration. Empirical relationships to determine the maximum size of a slick's area and the time required to attain this scale are presented and justified. The chemical degradation and dispersal of the slick is not modelled. Rather, the effort is applied to accurate forecast of the advection and surface dispersion. Mackay and Leinonen in [Mac77] produced a set of models to study the behavior of an oil spill on water. These models were coupled with observational data which was used  Chapter 5. Previous Work ^  40  to verify the model. Adjustment of the model parameters was also facilitated. Empirical in nature, it provided a step by step description of the analysis of the chemical processes which 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 spatial and temporal variations of the environment. Rather, the system was generalized to model environmental effects from scalar measures of these processes, such as wind speed and sea state. Other model parameters (such as volume and composition) are varied on a case by case basis, in constructing sets of model outputs to validate the assumptions that the models were constructed upon. The system modelled the spreading, dissolution, emulsification, evaporation and vertical diffusion processes over differing compositions and volumes of oil, wind speed, sea state, and the variation of certain model parameters. The separation of an oil slick into a thin and thick slick was formalized. This phenomenon has been observed in many oil spill cases. His models for determining the evaporation and dissolution flux from an oil slick 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 this flavour of oil spill modelling. In [Ven87] the model used and developed by the Canadian Atmospheric Environment Service was discussed. It was developed primarily to aid in short and long term cleanup of spill events. Designed to work on a general domain, the test region presented was the Bay of Fundy, and the model results were compared against observational data from the Argo Merchant spill of 1977. The CEAS model accurately modelled the transport of this oil slick. The Exxon Valdez spill's transport was simulated with this system in [Ven90]. The modelled drift was faithful to the actual transport of this major spill. The system models the advection, turbulent and physical spreading of the slick, and weathering by evaporation, and emulsification. The environmental processes were derived by accepting inputs from existing hydrodynamic and meteorological systems.  Chapter 5. Previous Work^  41  Al-Rabeh and associates in [A1R89] developed a stochastic simulation of oil spill fate and transport in the Arabian Gulf. A complete analysis and derivation of the problems of advection, dispersion and chemical degradation were presented. Advection and turbulent were diffusion modelled as [A1h75]. Surface spreading was applied to both the individual parcels, and to the oil body as a whole. Both the parcel and spill spreading was modelled around the elliptic paradigm. A statistically based method of determining the wind flow was used [Cek88]. Tidal and geostrophic currents in the test domain, the Arabian gulf, were derived from existing hydrodynamic models. The system modelled the physical spreading, turbulent diffusion, vertical mechanical dispersion, evaporation, and emulsification of test spills. A case study of an number of oil spills in the Arabian gulf was executed to verify the model. Over the past two decades the quality and breadth of the modelling of oil slick pollution has certainly improved. However, as most of the chemical degradation, dispersion and spreading algorithms are still highly empirical in nature and the improvements appear to have come from increased resources, and improved hydrodynamic and meteorological models. Works by [E1186] and [Nih84] both approach the problems of determining spread and degradation from a different tack. [E1186] models the slick as a distribution of oil droplets which are churned in the ocean by the wind and wave state. Droplets are constantly agitated into the water column, and are advected by the surface currents. The majority of the droplets rise to the surface due to buoyancy. As a result the spreading of the slick is described by the distribution of the oil parcels and shear diffusion. [Nih84] developed a non-linear model of the transport and spreading of oil slicks. His work simultaneously accounts for the effects of gravity, surface tension, friction, viscosity and weathering processes. The parameterizations were adapted to field conditions and are reported to have produced results which accurately modelled actual rates of spreading and advection.  Chapter 5. Previous Work^  42  Other researchers who have investigated the broad area of modelling oil spill pollution include [Tuc80, Lar88, Hun90, Yap89]. 5.3 Visualization of Previous Systems None of the literature reviewed discussed output methods for interpreting the various models. However, the nature of the problem suggests that visual methods were most certainly 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 composition had on these rates is shown. [Ven87, Ven90] shows the distributions of the particles rather than the slick body over time. The effect of varying the rate of the turbulent dispersion is readily apparent through these visualizations. In [Ven90] the grounding of parcels is shown along the Alaskan coastline. [A1R89] provides the most detailed visualization of their models output. Contour plots of the oil distribution over the domain region over time are presented. These displays gives a strong impression of the effect of the drift and spread on the distribution of the oil spill over time. Using similar methods the distribution of oil in the water column and shoreline are shown. This thesis will take the modelling knowledge provided by these systems, and extend it by adding a robust interface, direct model interactions, and presentation of the model states as visualizations. The interface will provide a simple means for the user to employ the system. Further it will shield the user from the internals of the models, and allow for data transparency and enforce organization on the system. Interactions will make the system both participatory and afford latitude in verification of models and of simulations. The visualizations provide a rapid means for interpreting the results, and are directed  Chapter 5. Previous Work^  43  towards 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 and regional planning, to computer science and geophysics, and yields a system which builds on the paradigm of modelling and simulating a marine oil spill. In this chapter previous work in the areas of modelling and simulating oil spills have been reviewed. The flavours of these systems range from the statistical and observational, to the strictly analytical. These techniques and methodologies were used in defining and implementing the oil spill visualization system.  Chapter 6  Modelling  6.1 Modelling This chapter covers the derivation of the models which comprise the system. The specific techniques applied to representing and modelling the environment and the oil spill process are described in detail. Further, the orderly flow of data through the system, and graphical issues such as the visual properties of petroleum, and interpolation schemes employed will be presented. 6.1.1 Restriction to the Ocean's Surface The system, and it component models, are designed to simulate only activities that occur on the air/water boundary, and at the land/water boundary at sea level. While many of the environmental processes extend high into the atmosphere, or deep into the ocean, for the 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 depressed by wave actions, the rise and fall of the tide, and the restoring force of gravity. Since the size of a domain can extend for many kilometers, and the vertical changes due to the wave state are small, it will be assumed that we are operating on a plane. These restrictions reduce the complexity of the models in both time and space, but concurrently reduces the flexibilty and accuracy of the system.  44  Chapter 6. Modelling^  45  6.2 Modelling the Environment The environmental models will consist of the terrain, tide, wind, and sea state within the confines of the region of interest. 6.2.1 The Terrain A representation of the terrain of the study domain was constructed. The topography of the Vancouver region is complex, formed by the Coast mountain range, the Fraser River valley, and a number of fjords, channels and island groups. An accurate model of the surface of this region was critical to many of the stages of the simulation and visualization development. Two representations of the local terrain were developed. A three-dimensional representation was constructed for use in visualizing the region, and for derivation of the primary representation which is used exclusively in the simulation. The three dimensional data for the region is collected from three sources, and ultimately merged into a irregular triangulation of (x, y, z) data. Most of the data for which  z > 0 1 came from a digital map-sheet of the 92G UTM block, which originated from the government body Energy, Mines and Resources [EMR88]. The data set is comprised of thematic 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 form of a two kilometer grid of depths represented in fathoms. With the inclusion of low level spot heights 3 , all three sets were merged and triangulated using the functionality of the ARC/INFO geographic information system. This system was used throughout the development of this system for all large scale tasks that involved the registration 1 Sea  level. and natural features, such as political boundaries, major roads and highways, forest cover. 3 To detail portions of the region that are below the 1000 foot contour. 2 11uman  Chapter 6. Modelling^  46  and transformation of large blocks of data. The triangulation was output into a general format for later use.  6.2.2 The Land Query The second representation of the local terrain was formed by resampling the threedimensional 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 given point is on land, or on this surface. This question or query must be posed frequently in the operation of the simulation and when sampling from the spatially dependant environmental models. As a result a criteria must be to answer this query in a fast and space efficient 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 which were either land or sea respectively. For an arbitrary point within the domain, to determine if it is on land or sea by performing a point in polygon test. This test requires some search be performed on the vertices, or the edges of the polygons. Preprocessing the polygons into an efficient structure allows for the query to be answered in time which is in the average case O(logn) [Pre85]. This cost amortized over a large number of iterations 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 of points, 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, this approach 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 grid  Chapter 6. Modelling^  47  with a resolution of 100 meters in both x and y. This step has a one time construction cost. 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 was considered. This structure could in some cases reduce the storage requirement of the grid. However, to answer the general query would require traversal of the tree. The average case depth of the resulting tree would be dependent on the complexity of the coastline at the sampled resolution. In the best case, single queries would suffice when an entire quadrant was contained in either land or sea. In the worst case, a branch of depth log N would be traversed, were N is the the maximum number of grid cells in x or y. Due to the inherent complexity of the shoreline in the region, this representation would 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 sea surface, or not, it can be encoded by a single bit. For each row the cells in groups of eight are packed into a single byte. For an N x M grid, only Iii i M + 1 bytes are required. A mapping of each point (x, y) into its grid coordinates (i, j) was implemented. This query will be refered to as  L(x,y) =  1 1 if z > 0 0 if z < 0  This yielded an almost constant time cost for each operation, independent of the size of 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 intertidal zone. If the maximum amplitude of the tide in the region is determined, then the intertidal zone can be viewed as all surface features for which the position is constrained by —C > z < ( where ( represents the tidal amplitude or height. If we extract the  Chapter 6. Modelling^  48  contours from the three dimensional model for z = —C and z = (, assemble the land water polygons, and construct the grids L_c and Lc, then a point will be contained in the intertidal zone if the following is true:  L_c(x,y) and not Lc(x,y). 6.2.3 Grid Data Description A mechanism was associated with this grid which allow registration to the model coordinate system. This allows for a correspondence between the grid and the earth's surface. It consists of a description of the dimensions of the grid, and translation, scale and rotation parameters. This structure is described by  G = {{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 the sides of the grid cells, T is the translation vector of the lower left corner of the grid cell (0,0) into the model coordinate system, 0 is the angle with which to rotate the grid counterclockwise about the lower left corner of the grid cell (0,0) to bring the positive y-axis into alignment with true north, S is a scale factor, where it is assumed the grid has the same spacing in x and y. G is an n x m grid, where each cell resolves to some structure. 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 coordinate 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 upset the operations of the system. The wind and tidal data structures, also based on grid structures, can now be described in a similar vein, which allows for determining the correspondence between a wind or tide flow and the underlying surface. In addition, as the  Chapter 6. Modelling^  49  model coordinate system is based on the Universal Transverse Mercator projection, direct registration of this data to the earth's surface is possible. The contents of the elements of a grid, can be any hypothetical structure, but are bits, scalars, and vectors for this implementation. A number of these grids will be defined. This grid description can be modified 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 is defined 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 N represents the number of samples. The samples are evenly spaced over time. If the simulation runs for longer than 24 hours, multiple copies of the sample set are used, and the first and last samples are assumed adjacent. The grid description mechanism allows for the orderly flow of data into and out of the system. Any data set for which this format can be determined can then be employed by the system. Similarly, data used by the system can be transformed into other systems by analysis of the grid description. The problem of data conversion into and out of the system is now external to the model.  6.2.4 The Tide Due the complex nature of the tidal generating force, and the interactions of forced waves with topography, no analytical solution to this system exists. However, there are empirical and numeric models which approximate the actions of the tide. Harmonic analysis involves the decomposition of observed tidal elevation data into its frequency constituents. This method was the first quantitative method for predicting the long term rise and fall of the tide. It is in general quite accurate, and is still the  Chapter 6. Modelling^  50  source of data for tide charts and tables. This method does not adequately describe the surface currents. Modelling the tidal current from observed data is difficult, especially in regions of complex topography. This is primarily a result of the difficulties involved in the long term collection of surface current data, and in the effect of residual currents on the the observation. The background noise introduced by wind and wave currents can obscure the value of the tidal current. However, if enough horizontal current data exist, and noise in the measurements accounted for, then with an appropriate interpolation scheme an estimate can be made of the current distribution over a region. Many complex interpolation schemes have been applied to sparsely located spatial data. Ahlstrom used data collected from the domain of interest using drogues 4 and fixed current meters to estimate the tidal current in his domain [A1h75]. Using the observation that the ebb and flood tidal streamlines generally paralleled the shoreline, with small changes of direction from flood to ebb, a base directional field for the tidal current was constructed. These were represented with line segments whose end points represent the direction at flood and ebb, and the length of the line representing the relative amplitude of the current between these given points. The amplitude of the tidal current was fixed to oscillate between the two end points of the line segment. With the advent of the computer, hydrodynamic models of the tidal forces and the resulting currents have been developed. The majority of the systems surveyed on oil spill simulations determined tidal currents either from a directly coupled hydrodynamic model, 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 and the surrounding marine region were provided by the UBC Department of Oceanography. 4 Free  floating current gauges.  Chapter 6. Modelling^  51  These constituents were extracted from a hydrodynamic model of estuarine circulation in the enclosed waterways of the Pacific Northwest developed by [Cre88]. The system in its 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 contained by these water bodies). The existing three dimensional model represents the culmination of a series of tidal models developed over a number of years by Crean and associates. The initial model was one dimensional, and was used for determining the M2 tidal constituents through the domain. A grid was constructed over the region, with each grid cell representing a cross section of a bay, straight or channel. The connectivity of the grids was based on the topology of the conveying channels. The equations of motion with friction present were solved using a finite difference scheme. Using many years of observational data, the parameters of the model were adjusted such that the energy losses of the waves in the model and the real world were in approximate agreement. Through the years the model was extended into a two-dimensional and finally a three-dimensional scheme based on a regular two kilometer grid. For a detailed account of the methodologies and mathematics applied, see [Cre88]. This provides us with a time-based analytical approximation of the varying tidal elevation and surface currents over the test region. To determine the velocity in the u directions at cell (i, j) at time t we need only evaluate the function c  E I3ijk cos(`2irfo + oiik) J...i. is the frequency of constituent k at cell (i,j) is the amplitude of constituent k 6 See  Chapter on Environment positive u direction is east and the positive v direction is north, or below mean tide. 6 The  j  represents the evelation above  Chapter 6. Modelling^  52  Oink is the phase of constituent k C is the number of harmonic constituents t is the time in hours measured with respect to an absolute clock Similar evaluations are required for determining v and C. The beauty of this model is that it is stable, has a compact representation, and is quite accurate. [Cre88] reported that 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 real world. For a given grid cell, there is no direct mapping of the surface of the grid cell to its corresponding region on the earth surface. This had to be expected, as the grid cells are aligned with the channels and other coastal features. A channel may be represented by a single grid cell in the model, while the actual channel may only be one kilometer in width. In addition to coarse spatial resolution, the model was developed at a time when core memory was at a premium. It requires 7 that care must be taken when transforming data from the grid cells to the earths surface. A subset of the tidal model was extracted that covered the study area, and only grid cells that could be aligned with the terrain model by rotation, translation and scaling were chosen. This coarse registration was performed by hand using the resources of the Geography department. 6.2.5 The Wind  Determining the distribution of the local winds is the most important component of a oil spill prediction model. As was previously discussed, the wind largely dictates the overall fate of the spill. Its effect is felt in the drift, spreading, and weathering rates, and it is the 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 straightforward  Chapter 6. Modelling^  53  methods consist of interpolation of observed data, and the solving of numerical approximations to the meteorology of the region. Both of these methods have advantages and disadvantages. To correctly navigate through this area, both methods must be described. 6.2.6 Empirical Wind Modelling  Wind velocity data tends to be sparse and irregular in space. Further, if the recording device is not continuous, then the sampling will be sparse over time. Each set of observational wind data represents the state of the winds for the particular conditions that existed at the time. A diurnal cycle over the domain of study was developed by interpolating a set of wind data from measurement stations over the domain. The data's source was Environment Canada, and the date of the sampling was the 17th of July 1992. On this day a stable high pressure ridge had been in place for a number of weeks, resulting in typical daytime and nightime temperatures for the region. The data was collected from 10 recording stations spread out over the domain. The positions of these recording stations were used to generate an irregular network. This network was interpolated using bilinear interpolation at each of the sampling times, and sampled into a regular grid. With the inclusion 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 of the wind could be achieved. This problem is addressed in a modelling system developed by [Cek88]. Data from a three dimensional wind model was generated for prevailing winds and indexed by eight directions and seven wind speeds. This modelled data could be replaced by observational data, although a larger bulk of data must be collected. This modelled data was interpolated into a regular grid and placed into the database. Using observations of the direction and magnitude of the prevailing wind, an estimate of the wind patterns at a given time could be interpolated from the database.  Chapter 6. Modelling^  54  6.2.7 Numerical Wind Modelling Numerical modelling of wind fields has been an active area of research for a number of years. It is an intricate science, from which promising results have been derived. While these models have been developed to determine a wide range of atmospheric motions, the focus in this thesis will be placed on mesoscale wind modelling. This area is concerned with determining the wind patterns which are in the scale of one to 100 kilometers, and of which have a large influence on the spread and drift of oil spills'. Wind patterns at this scale near the earth's surface are largely influenced by channeling and deflection caused by topography, friction, and exchanges of heat and momentum between the lower atmosphere 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 which do not consider all of the physics involved. The models rely on solving tendency equations for meteorological flow. These tendency equations consist of formulations of the horizontal transfer of momentum under the influence of friction, heating, the Coriolis force and the pressure gradient. There are three types of generalized models currently developed for diagnosing wind fields over terrain. Mass conservative models rely on observational data for initial conditions, and then interpolate these values into a regular grid, and apply least squares methods to the conservation of mass throughout the model. Multi-layer vertically integrated models divide the atmosphere into four layers. The tendency 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 the 8 See  Chapter 3.  Chapter 6. Modelling^  55  tendency equations. These models are not mass conservative. Changes in pressure are determined hydrostatically 9 by parameterization of the atmospheres structure, along with boundary layer friction and diabatic A sigma-coordinate based wind modeller was was provided by Keith Ayotte of York University and is based on [Ayo87]. The system was developed in the Department of Geography at UBC. This model is of the third type, relying on a vertical parameterization of the atmosphere. The initial conditions which the model requires includes a grid of the regional elevation, and the heights and temperature of a reference pressure leveln. The friction 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 threedimensional surface model. Estimates of the Bowen ratio and the roughness length were taken from [Ayo87]. Heating and cooling levels were adjusted in an attempt to estimate the wind patterns that would result over a diurnal cycle. The model was found to be very sensitive to grid size, and to the heating-cooling parameters. The minimum grid cell length for which the model would converge was found to be 18 km, which only covered the domain with an eight by six grid. 6.2.8 Interpolation of Wind Field Data Given that a surface wind field is available, the surface current due to the wind must be determined. The physics of the shear induced by the wind and the resulting current is discussed in the Appendix. Given a wind velocity 17 - at a given point, the direction and 9 The hydrostatic assumption states that the vertical component of the pressure gradient force is balanced by the force of gravity ie a = —pg where p is the density of the air and g is the gravitational force. '°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 approximately 1500 meters  Chapter 6. Modelling^  56  magnitude of the induced surface current was given by,  V* = crMe 17 -..  a = wind coupling coefficient MB = rotation matrix  The deflection angle  Ode f  is set to 0 for the simulations presented. Spaulding reports in  [Spa88] that many researchers have had success in modelling the transport of oil slicks when applying a deflection angle of 0. 6.2.9 Interpolation of Wind and Tide The tidal and wind models provide reasonable approximations of the flows over regular grids. 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 be considered. v and C are determined in the same fashion. A u surface will be constructed from a sampling of the flow in the x and y direction by using (x, y, u(x , y)) as the surface points at time t. The grid points determine the corners of each rectangular patch. Each patch 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 s and 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 than the tidal grid. If a grid cell in the tidal model has no constituents, then the flow is set as a default to (0, 0). However, in reality, the flow due to the tidal current will be non-zero up to the land/sea interface. As a result, there is disproportional weighting towards a zero flow near this boundary. This problem is described in Section 6.1.  57  Chapter 6. Modelling^  Water  Ui 1 j 1 -  -  Figure 6.1: Bilinear Interpolating of Tide Flow This interpolation scheme restricts the points of the grid to be on the surface, but it is not C 1 or C 2 continuous. Catmull-Rom interpolation, developed in [Bre77] was applied to enforce smoothness on the surface. This bi-cubic parametric interpolation scheme enforces C 2 continuity, and places the surface through the sample points. To determine u at a point (x, y) the patch it is contained in is determined. From (x, y) the values of .9 and t can be determined. The value of u is found by evaluating  u = :§cmcTtT where (s) and^are the vectors [s 3 s 2 s1] and [t 3 t 2 t1]. C is the Catmull-Rom kernel matrix, and M is the geometry matrix formed of 16 control points from the surrounding patches [Fo182, Bre77]. Temporal interpolation of the tide field is not necessary as the model is explicitly based on time. Interpolation over time is required for the wind field. This is achieved by linearly interpolating between two sampled fields. The wind fields are sampled N times over a diurnal cycle. For an arbitrary time t at (x, y) the flow is determined by using bilinear interpolation to find the flow in the previous wind field, Vprey and the flow in the  Chapter 6. Modelling^  58  next wind field Vnext • The flow at time t is found by calculating Veurrent where Vcurrent  v  = VVprev + (1 — V)Iinext  =  t—  prey  tnext — tprev  where t prev is the time the previous wind field was sampled at, and  t- next  is the sampling  time of the next wind field. 6.2.10 Other Advecting Flows The domain is sheltered for the most part from large scale geostrophic oceanic flows occurring in the Pacific [Cre88, A1h75]. The flow from fresh water runoff is significant in this predominantly estuarine region, largely as a result of the contributions from the Fraser River. This fresh water runoff has an effect on the temperature and composition of 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 Factors The state of the sea surface is important in determining the various process rates effecting an oil slick. It determines the turbulence and the amount of available mechanical energy on the ocean's surface. It is largely determined by the prolonged force of wind applied to the ocean's surface. Qualitative description of the state of the sea have been made by people since they first began travelling by water. Over the years these descriptions have been somewhat quantified in the Beaufort scale. This scale is calibrated against both the velocity of the predominant wind, and a description of the apparent wave state. For example, a Beaufort value of 7 defines a range of wind speeds 12 from 22 to 25 ms -1 , and is described 12 Around  55kh-1  Chapter 6. Modelling^  59  as 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. Wave prediction is based on the assumption that surface waves are caused by a wind force which 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 a long period of time. This fully developed sea will require a stronger wind to attain larger waves, and a weaker wind to diminish the surface. As an approximation to the state of the sea surface, a parameterization of the Beaufort Scale was developed. It was formulated to respond to changes in the average wind speed over the domain. This is achieved over time by the formulation St = V3(lUt l) + (1 — q)St-At where  i  is a weight in [0, 1], MI is the average wind speed at time t, and St is the sea  state at time t, and B is a function which returns the sea state number based on the wind speed. The sea state is dependent on the decay of the previous state, plus the energy contribution from the wind. Ranging the parameter q changes the effective contribution of the previous sea state and the present. 6.4 Modelling of the Oil Spill Process The models assembled will provide an approximation to the general transport equation  ap = _rjvp+v(ichvp)+ 0 at  (6. 1)  where p is the density of the pollutant, 0 is the net advecting current, Kh is the combined effects of dispersion and physical spreading, and 0 represents sources and sinks  Chapter 6. Modelling^  60  [Ven87]. This system can be approximated by solutions to the differential equations that it is described by. These methods tend to be slow and potentially unstable. As the intent of this system is to provide an interactive mechanism for abstracting oil spills, a Lagrangian particle tracking method is employed. This method is more efficient because it does not depend on the system converging on a solution, and is stable [A1h75]. It is also more efficient in cases where a large volume of petroleum is involved, and where the particles tend to occupy only a portion of the domain [A1R89]. Each of the terms of Equation 6.1 will be determined and the state of the slick at each point in time determined. The first term —0V p represents the advection of the slick from the wind driven and tidal currents. The second term V(KhV p) indicates how much the slick spreads and disperses over time, and the final term, 0, represents the amount of oil which is added to the slick from a venting petroleum source, and the amount taken away by weathering and grounding processes. These terms are determined at each iteration of the simulation. The initial conditions are contained in a representation of the oil spill source S.  S . {Po,V0,C, f, t o } where Po = (x, y, 0) is position, V is the volume, C = (p, M0 ...MN) is the composition of the oil, and f is the flow rate, and t o is the time when the source will begin to vent into the marine environment. The position Po can be located anywhere on the ocean surface. The volume Vo is expressed in liters, and the flow rate f in liters per second. The composition C of the petroleum is reduced to its density p, and the respective mass of each hydrocarbon fractions Mo through MN. These parameters are used primarily when determining the weathering processes. The oil slick will not be modelled as a continuum. It will be represented as a set  Chapter 6. Modelling^  61  of parcels, each of which is independent. It is the nature of the Lagrangian particle system to consider the slick as a collection of distinct entities. Each parcel will carry with it its volume, composition, time stamp, and a distinct label, and is processed by the simulation as a micro-slick with these properties. The collection of particles is placed in a list structure, ordered by label, which will be refered to as LS. 6.4.1 Oil Release The oil is released at Po at time t o . For each iteration a new volume is released, until the source is depleted. The volume that must be released is determined by V, = fAt where At is the time step. This volume is then allocated to particles based on a parameter Vmax which is the maximum volume which can be allocated to a particle. The number of particles released at a given time step by,  Not =  lv÷..xj +1 if Vs > Vmas 1^if Vs < Vmax  As No t will likely not be an integer, at a release where Vmax << Vs the volume of the 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 onto the slick list denoted Ls, and the volume of Vs decremented. 6.4.2 Transportation The transport of the slick at a given step is found by determining the tidal and wind driven current at each position currently occupied by a parcel. Duplicates of Ls are created, and passed to the wind and tide interpolation routines. The lists LT and Lw represent the tidal and wind flow at each point occupied by particles of Ls. The new  Chapter 6. Modelling^  62  positions for each particle are determined by  Pi+ i i  = (T(Pii )  + a Modef [W (Pii )]) A t  where Pi; is the jth particle position on list Ls at step i, T(P) is the tidal velocity at position P, W(P) is the wind velocity at P, Ro def is a counter-clockwise rotation of the wind 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 to wind stress at that point. 6.4.3  Turbulent Diffusion and Spreading  There have been many empirical models put forward over the years for modelling the physical spreading and turbulent diffusion that affect a petroleum slick. Early methods by [Fay71] concentrated on determining the physical spread from the spreading force of gravity and the retarding forces of viscosity and inertia. The model splits spreading into three regimes. The first regime involves spread under the force of gravity, with inertia as the retarding force. In the second phase, the inertia is overwhelmed by the viscous force as the primary retarding force. In the final regime, the oil has become so thin that the effect of gravity is no longer felt. Instead the net force resulting from the difference in tensions between the water, air and oil dominates the spreading. Mackay and Leinonen derive an empirical relationship based on radial spread using Fay's formulations for the radius of a slick over time in [Mac77].  c/ 3 Ci tVo (pin — po )po rt = ro + 7rptu  where ro is the initial radius, C 1 is a constant, t is the time in seconds, Vo the original volume and p w and p o the densities of the water and oil respectively.  Chapter 6. Modelling^  63  Methods which assume the physical spreading to occur radially must also presume that it is occurring on a calm sea. These models do not account for weathering processes or turbulent shear. To model the spreading of oil in a more realistic way, both physical and turbulent spreading must be modelled. Methods for determining the turbulent diffusion of a slick body either involve solutions of the differential equations describing this process or use statistical 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 selected for use because it acts upon discreet quantities rather than a continuum, is rapid and stable, and is easily parameterized. This model assumes that the particles diffuse under the influences of near Brownian motion resulting from turbulence in the transporting medium. Diffusion occurs in both the vertical and the horizontal plane, although we will only consider the horizontal term. A parcel will travel on average a distance described by Dh = ^47hAt  where At is the time step, and -yh is the horizontal dispersion coefficient which is a combination of the eddy diffusivity and the spreading diffusivity For an arbitrary parcel the 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 from 79 = 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 a random diffusion step to each particle. These steps are generated for each particle of Ls and placed onto a new list labelled Ldif , which will eventually be added to the positions of the particles on Ls. The value of the dispersion coefficient -yh has been reported in  Chapter 6. Modelling^  64  the 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. Many researchers have observed that oil spreads in an elliptic fashion. Al-Rabeh and associates use the following formulation to determine the elliptic spread of individual particles and the slick body as a whole in [A1R89]. This formulation is shown in figure 6.2 A = 71--QR Pw — Po Q = C(^)3 VIt4 Po R = Q - F aWttt  ^1 i 1  where 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 the direction of the wind. This spreading is applied to both the individual parcels and the slick as a whole. The slick is considered as a nest of similar ellipses. If we consider the wind to be blowing in the positive x direction, with the center of mass of the spill located at the origin (as in figure 6.2) then the outward spread of the individual particles is determined by  dx = Sixdt dy = Sy ydt ,Si = xdR R ydQ  Sy _  Q•  The model employed in the simulation system uses a combination of radial and elliptic spreading. However, elliptic spreading is only applied to the entire collection of particles to determine the overall spread. The physical spreading at a particle is assumed to be  Chapter 6. Modelling^  65  Individual Parcels  ,  P /  S.  .^  .^Boundary N.  ....  .^  Boundary at VI-oft^  at  t.  .  Direction of Wind  ,  .  . Displacement of Parcel By Elliptic Spread  Figure 6.2: Spreading of Model Slick as Concentric Ellipses radial. This was done to avoid calculation of the long and short axis of the hypothetical ellipse on a parcel by parcel basis. This can be justified based on the observation that the elliptic spread of the individual particles will not greatly affect the overall distribution of the oil. Rather, it is the elliptic spread of the entire oil body that will determine the shape of the slick. One question that remains to be answered is the effect on this algorithm of the point selected as the center from which to spread from. At present this is determined by calculating the center of mass of the particles. However, in the presence of strong currents and winds, the center of mass may not even be near any parcels. The spreading will often draw the slick bodies out into frayed and linear shapes. Spreading about this point seems a bit simplistic. The spread about a medial line could help to relieve this problem in certain cases. However, the real problem comes from the representation of the  Chapter 6. Modelling^  66  slick body as a collection of particles. There is no fast and concise method of determining the shape of the slick from the particle representation. The slick may not even be a single entity, 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 or less instantaneously. Further, its effect should be modified when the slick body begins to elongate in the presence of the local currents. The only other consideration that must be made is the duration of the spreading applied to a slick. This was done by halting the spread at the particles once the thickness had reached li m i n, = ,47-. = 10' meters, or a hundredth of a millimeter. While many spills achieve slick thicknesses which are monomolecular, most large spills are reported to have stopped spreading far before before achieving this thickness. To summarize, the spreading process is modelled by applying the diffusion algorithm developed by [A1h75], in conjunction with a modified elliptic spread function from [A1R89].  6.4.4 Evaporation and Dissolution In order to determine the evaporative and dissolution flux from a slick body, either a psuedo-component or analytical approach must be used. [Mac77] proposed a psuedocomponent approach, which separated the oil composition into its main fractions. Each component 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 is determined, and summed over the entire mass  N=  N N  E IV: = EIC x P:RT i=1^i=1  e i  Chapter 6. Modelling^  67  where Arf is the evaporative flux of component i, ife is the evaporative mass transfer coefficient  13  ,  xi is the mole fraction of component  i, Pia is the vapour pressure of the  component i, R is the gas constant and T is the air temperature above the slick. The mass lost to evaporation proceeds into the atmosphere, which is assumed to be an infinite sink. Other methods avoid decomposing the oil body, relying instead on averages to the evaporative flux. These methods have the advantage of requiring less information. The evaporation rate has been described as a function of area, wind speed, slick thickness and the roughness of the sea surface. For the purposes of this system, the composition of the oil will be generalized, and Mackay's algorithm applied. This algorithm is selected because it considers the major factors involved in evaporation, and is compact and stable. The ambient air temperature will be considered constant for the simulation. The components of the oil will be generalized into a volatile and non-volatile fraction, by boiling point. If a component has a 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 Grounding  When a parcel drifts into contact with the shoreline, a decision must be made as to its fate. In general the entire mass of the slick will not ground. Only some fraction will ground, 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, and the 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 with the 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-1  Chapter 6. Modelling^  68  the parcel is removed. The volume that has grounded is kept on a separate list LG. To determine if a slick grounds, the superposition of the effects of advection, diffusion and spreading are determined in the vector U. If the addition of this vector to the position of the 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 U is alternately rotated counterclockwise by .6.0 and constricted by AIM. AO is a random ir variable sampled from a uniform distribution in the range of [41 . The constriction of the vectors length is generated from a [R]lou l uniform distribution. These parameters were set to these ranges to minimize the number of steps required for a parcel to find a point contained on the ocean's surface. This process is shown in figure 6.3. Original Displacement  Figure 6.3: Particle-Shoreline Avoidance  6.4.6 Weathering Processes not Modelled  There are many processes which are not modelled in this system. Most of the effort was concentrated on determining consistent advection, diffusion, spreading and grounding. The primary reason for not developing a more-in-depth weathering model is the level of detail required to attain a truly consistent system. Complex chemical descriptions of the petroleum source are required. Methods for determining emulsion formation,  Chapter 6. Modelling^  69  dissolution, 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 they have a small effect on the mass of the slick, and occur at slow rates late into the spill's life cycle. Further, the processes that drive them are poorly understood, which is reflected in the lack of available models. 6.5 Determining the Time Step An appropriate time step must be determined for the simulation. Too large a step will destabilize the oil spill processes, and in particular the diffusion and spreading. Further, the temporal variations of the environmental conditions must also be considered. Too small a time step will cause poor performance time. [A1h75] states that if smooth and continuous performance is required, then the time step should be limited such that the maximum displacement of a parcel over time be less than the minimum distance between the data points in the velocity field. In our case, these separations are not fixed due to mixed grid sizes. As an alternative the maximum distance that a particle may travel is constrained by the parameter  Dmax.  The time step is found by determining  At =  Dmax  iVi.  where 11/1„„„ is the maximum speed of the total advecting current. The time steps generated by this constraint are determined on an iteration by iteration basis with good results. 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 by k  tabs = to +  E ti i=i  where k is the iteration number, and ti is the time step determined at the ith iteration.  Chapter 6. Modelling^  70  In this chapter models for representing the petroleum, environment and spill have been presented. Further, representation and interpolation of data over the spatial and temporal domain have been elaborated. A subset of these models were selected, developed and integrated into a the oil spill visualization system.  Chapter 7  The System  Given the collection of approximations to the various physical and chemical processes involved, systems must be developed for combining and managing these models to produce simulations and visualization of the model problem. These systems consist of a database, a simulation/visualizer and a batch or extended visualizer. 7.1 Interfaces Graphical user interfaces were developed for both the database and the visualization system. It was assumed that a potential user of the system would be aware of the problems associated with oil spills and the associated processes in the marine environment, and has some knowledge of interactions with a computer system. The criteria for development was to provide a front end that would allow for not only the rapid development of model simulations of oil spills, visualization of this event, and variation of key model parameters, but of more mundane tasks, such as extracting subsets, and saving model states. 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]. Implementation details can be found in the appendix. This application software allows for the interactive 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 Graphics  71  Chapter 7. The System^  72  and wind flows, or to set up reference grids or other display oriented features. Through these menus, many of the tasks involved in setting up and synchronizing the system remain transparent to the user. The major drawback of the FORMS software is that it extensively uses the Silicon Graphics font manager. This effectively binds FORMS, and the database and visualizer to operate only on systems which support this software. The extended visualizer is batch oriented. It was developed under the assumption that a user would eventually develop animations and detailed renderings of the various processes and oil spill simulations. By using a simple script language, the user can specify the location and type of data that is to be processed, the output required, how many frames to generate, and the relative time from which to base these frames. These script files can be run as background processes. 7.2 The Database The database contains the state data for the tidal model, the grid representation of the terrain, and wind flow data. The front end constructed over the database shows the user where the wind and tide models are defined spatially, and allows for the extraction of subsets of this data. It was deemed reasonable to provide the functionality for subset extraction for three reasons. The first is that the overall wind, tide and terrain models are quite large. Execution of the entire system is a strain on the fastest of the available workstations, 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 ten million liters of oil into the Gulf of Georgia, then a large portion of the entire model may be extracted. However, if they wish to study the effect of spilling ten thousand liters 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. The  Chapter 7. The System^  73  option, however, is available to work on the entire domain or a subset. The final reason is security. It is a dangerous precedent to allow users to be directly operating on a root set of data. The possibility for corrupting portions of the database could occur, which would 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 a local directory in the user's workspace. 7.2.1 Data Storage and Extraction The 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 value of the tidal flows and the elevation at the given location. The numerical and empirical wind models produced data sets which did not have such compact representations. The numeric model requires large amounts of processor time before it converges on a solution for a given set of initial conditions. As a result this model cannot be directly coupled to the interactive simulation and visualization system. If it were, it would not only grossly tax the workstation which was running the system, but also would require the user's attention, in terms of providing initial conditions, and sequencing fields. This would make the inclusion of wind field data far from transparent. As an alternative, the wind fields 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 was used for the simulation as the numeric model had too coarse a resolution. The wind field data is generated at some temporal rate, and presented to the database in the form of a grid description, synchronization data on the various fields, and a set of wind fields. The database merely extracts a subset of these windfields, and recreates the local grid and synchronization data in the users copy of the data. The components of the tidal model  Chapter 7. The System^  74  and terrain grid are extracted in a similar fashion. This methodology provides users with a consistent means for the inclusion of data from other systems. The actual current field data must be placed in the system's regular grid format, and a spatial and temporal description provided. In the case of presenting new 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 a geostrophic current. Changes in software would be required to replace the current tidal model, 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, the grid description of the data is used. When a user selects a region, all sample grid points that lie within this region (with data padding necessary in some cases due to boundary effects) are extracted. The subset is rectangular, and oriented with the horizontal edges of the rectangle parallel to lines of latitude, and the vertical edges aligned parallel to lines of longitude. The orientation and shape restriction allows for the rapid extraction of subset data. The front end of the database provides little additional functionality, other than the ability to zoom and pan around the region, to change the path with which to store the data, and the naming of a given data set. 7.3 The Interactive Visualizer The 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 the surface flows due to the tide and wind interactively. Descriptions of the models, sampling schemas and the flow of the simulation are more extensively covered in chapter 6. The system is presented to the user in the form of an input menu, and a graphic display of  Chapter 7. The System^  75  the current domain. The interface is a hierarchy of menus, each performing a specific task. 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 extracted from the database. The user can rapidly peruse the available regions from a browser, with the graphic display updated to reflect the current region. The structures and data of the domain models are not loaded into the system until the user explicitly performs this operation. From this junction the user can either execute a simulation within the confines of the loaded domain, or visualize the major flows.  7.3.1 Running a Simulation The simulation panel allows the user to set the initial conditions, sample from, and execute a oil spill simulation. The initial conditions consist of Vo  = inital volume of the source  Vp  = maximum volume per parcel  f = Po  flow rate in liters per second  = position of the source within the domain  t o = start time of the simulation The volume, parcel volume and flow rate are set by manipulation of a slider. The maximum volume that a parcel may contain dictates the number of parcels which will be released. 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 the number of parcels released. Sensible bounds are enforced on these sliders. To determine the initial position of the source, the mouse icon must be placed over the corresponding location in the graphic display and the appropriate key depressed. This will place a  Chapter 7. The System^  76  symbol permanently at that location within the graphic display. This operation is not permitted either outside the domain, or on land. The starting time of the simulation is provided through the time menu. It consists of parameters for year, month, day, hour and minute, within appropriate bounds. There is also a set of model parameters external to this menu. There adjustment is possible though a parameter file that is placed in every users workspace, or changes to internal variables in the program. The setting of their values is not strictly speaking a requirement for running a spill simulation. Further, they are used in other stages of the overall model. A description now follows of these parameters, based on their purpose.  Te = tidal component switches a = wind coupling coefficient QT = tidal flow constraint coefficient 17  = sea state weighting parameter  7', determines which of the harmonic constituents of the tidal model are to be used in the current simulation. The model includes 14 major constituents, which can be turned off and on through T. a is the wind coupling coefficient. It determines the magnitude of the drift current due to the wind, and is usually set in the range of 0.01 to 0.08. ay is analagous to the wind coupling coefficient. It has no physical basis, and is merely in place to allow for adjusting the velocity of the tidal current. The parameters a and QT can be employed to adjust the effect of the wind and tide, if necessary. The weighting term 9 is used in determining the sea state. It is constrained to range between 0.0 and 1.0. A small value for 9 places more weight on the previous state of the sea.  Chapter 7. The System^  7h  77  = horizontal dispersion coefficient  Po = petroleum density Pw = water density w = splitting rate for grounding parcels Dmax  = The maximum distance a parcel may travel in the advecting flow  ATmax = The maximum time step allowed for each iteration  The parameter ryh defined in Chapter 6 determines the magnitude of the step size that a parcel may undergo under the influence of horizontal dispersion. As it is largely a function of the upper layer turbulence in the ocean, it will be combined with the sea state in Chapter 8. The terms p a and p„, the densities of the oil and water, will determine the spreading rate. The parameter w will determine how much of a parcel will ground when its advection, spreading and diffusion brings it in close proximity to the shoreline. The parameter D mar is used to determine the length of a time step on an iteration by iteration basis. The term AT,,,„„ is set to avoid large time steps. Large time steps cause erratic and visually unappealing events to occur in the simulation. Other parameters and combinations of parameters will be discuss in Chapter 8 and in context of the visualization process.  7.3.2 Displaying the Simulation Once the setting of parameters is complete, the user may execute the simulation. The slick is 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, spread  Chapter 7. The System^  78  and potentially grounded. The rate at which the iterations are displayed is dependent on a 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 very jumpy and discontinuous. Grounded parcels are differentiated from sea borne parcels by colour. 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 Visualizer An environmental monitor can be switched on during the simulation. This monitor reports the general state of the environment in a noninteractive panel separate from the actual graphic display. It provides a synoptic view of the state of the environment on an  Chapter 7. The System^  79  iteration by iteration basis. The current time, and the direction of the prevailing wind are 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 bar graph. The elevation of the tide and the sea state are reported in a similar manner. The update rate of this monitor can be set in a fashion analogous to the frame rate. This monitor is a major tax on the available resources of a work-station. To achieve maximum performance it must be turned off, or its sampling rate set very low. 7.3.3 Sampling from a Simulation A key-frame can be defined as a sampling of the system taken at a particular time. This sampling is sufficient to generate a graphic image of the simulation. A key-frame option is available on the simulation panel. Prosaic tasks like specifying the path and names for the resulting slick scripts must be specified. The parameters for sampling are  Nk = number of key-frames to sample  At = time step between frames tk = time offset  The system attempts to perform the sampling in steps of size At from initial time t o + tk to t o + tk (Nk — 1)At. However, because of the adaptive time step in the simulation, this is not always possible. It will sample when the clock exceeds or is equal to the next sampling time. Each sample is placed in a script file with an extension denoting its position in the series 0 ... Nk - 1. Each script file contains the time of 2A  cardinal direction is one of the eight points on a compass, ie N,NW,NE etc.  80  Chapter 7. The System^  sampling, a default eye position', and the identifier, position, volume and radius of each parcel, separated as either being seaborne or grounded.  7.3.4 Flow Visualization The flow visualizer allows the user to view the flows generated by the tide and the wind models over the currently selected domain in an interactive manner. The input panel allows the variation of the sampling, scaling of the velocity field, and selection of a style in which to view the flow. The sampling of the flow is performed by laying a regular mesh of N x M sample points over the domain. The dimensions of this grid are provided by the user, and N and M need not be equal, but necessarily N, M > 1. The temporal sampling rate can also be set. The tidal flow in the domain rarely exceeds 3ms  -1 .  The velocity of the wind can  reach 30 — 40ms -1 if gale force winds were modelled. However, in general velocities do not exceed 10ms -1 . When flows of this magnitude are viewed within the framework of a domain which may be 100 to 200 square kilometers, the dynamic patterns of these systems 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 slider for both the wind and tidal flows is provided. The scale of the two flows can be adjusted independently. This scaling allows the viewer to observe the more subtle interactions occurring 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 common method of viewing a surface flow is by representing each point flow by a vector. For each grid 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 visualizer when the frame was extracted  Chapter 7. The System  ^  81  k t  • ..3'^3 --....■_■ - "4.-■____,,..., • -■^ri■-.././1■—.4„,„,,,i. e -  ,■■ 1.- ■'s.,-"  ,  ,,  6  ete. /..-4-  i 1•■;--...re//i."-,  it I 1  1 1 .1^■  .  , ,^1 4,  1!  1  4 ,..e. 1/./...., wor , ..4,..,..„.  1  1^I .,^  i i 1^% ‘4.• , ..e.:. ." si■NiNAAN„^ / ....wAg I 1  ..._•-....,....j...■^\ \ NV. \l`• --  rr  '...■....;\^\-.4's \\*\  '4"...:%"^NI"\7"..' __+-■-.■.■• i^ ^ '''''''''\%\"\l' ,^,  -1,. . . .. k• ^ .., . '1,  ..... . „  0, - N-4-k  ift  •  ,  1  •  , A. •  Figure 7.2:^Tidal Flow, Time 0900:17-Jul-92  Figure 7.3:^Drift Method, Time 0900:17-Jul-92  flows is rendered. The base of the vector is the sample point on the mesh, with the head determined by the addition of the local flow to the sample point. To view some of the more intricate flows in action, the sampling mesh must be increased in resolution. This can however lead to relatively confusing images, especially in areas of slowly changing smooth flow, caused by the overlapping of vectors. To address this problem, the flow can be displayed with only the head of the vector displayed. This displays the flow as a collection of offset points. The third method approached both sampling and rendering of the flow in a different fashion. At each sample point in the mesh, flow samples are released at a fixed time interval. It is allowed to advect with the current, and at the next time step, another sample point is released. Both samples again undergo the process of advection. If the initial sample point is included, a path of length 2 can now be formed between the sample points. This process is analogous to standing in a vessel at the sample point, and  Chapter 7. The System^  ..-----,,,-''..-*----,- ..,--  -  1-7 .P''e"" ."-rr .--...r^..-4""  0  ! '/. '  ,,,er. , r  "  4,P#  .,,^r^  4._ . . . . 1  ‹..e'...-.. "' ----'^--. r^ :;  82  —4,4  * ..,--4''—?--,— ---- ..--- ■  „ „.,10 ......., .e .  ,,,,,,,,.  ''.  ,,... ,..,  :orer-je'r.,,, 1,-1.er"..... 1 ^..■^"'"°:,........,••••  "  -,..4?,...a..4.—••....--  ,#.  P ,..'<•^4'...  .F,'4SP.......:*<  ,.....r"1  "  ^....''''':■riel^.....,'"■•••■0"...,•4  .,,rd''''^•■■,•°''..^  ....:..".. $  ...,...r:.........":„...., i a.r.'"'''''rd....r' r■..'''ro•4....r•  ,,  #,....r,,■+!• .......r...•••••• msf^i' ..... ...■....  4.......  ....412}".4....0^......  .•  ^• •  Figure 7.4: Wind Flow, Time 0900:17-Jul-92 dropping 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 patterns of 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 and direction over time. Examples of these visualizations are shown in figures 7.2, 7.3, and 7.4. The time of these snapshots and all other environmental images corresponds to the standard spill simulation parameters discussed in chapter 8.  Chapter 7. The System^  83  7.4 The Extended Visualizer The purpose of the extended visualizer is to produce detailed renderings of component models, and oil spills simulations. These renderings are more computationally intensive than the simple visualizations presented in the interactive system. To allow the user to communicate with the extended visualizer, a simple script language was developed. Each script file contains information on the location of the data set, and the name to give to the resulting images. This name will be suffixed with a frame number. Next, the type of rendering that is requested. At present there are six types 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. The resolution of the image to be rendered (where applicable) is specified by (N, M, D) where  N and M are the rows and columns of the result image, and D is the number of bytes used to encode the pixel value. At present D = 1, or 8 bits for the current set, but extensions to include colour could be included easily. The start time of the initial frame and a time step are provided. In the case of tide-specific rendering, the image is rendered directly from the tidal model. In slick-specific rendering, the time is encoded in the slick  Chapter 7. The System^  84  script. Both time and time step can be set to —1 as a signal to the system to used a default time and time step. The final data that must be provided in the script file is the start and end frame number. In the case of rendering slick data, a slick file with the suffix must exist.  7.4.1 Determining the Thickness Map In the interactive stage of the system, the user views the slick as a collection of parcels which are acted upon by the local forces. This is a discrete approximation of how an oil slick 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, and will resist mixing with the surrounding ocean until a force is applied. If a slick was viewed on a calm day, it would remain a single body, spreading radially. In the presence of stronger local forces, the slick would elongate and break apart into a number of disjoint slicks, which may eventually recombine. With very strong local effects, the body may disperse 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, which spread, advect and disperse in the same manner as the underlying collection of slick bodies. The task now is to recover the shape of the oil slicks, and determine an approximation of 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 described by its position P, volume V and radius r. A regular mesh of p x q samples is generated over a rectangle which contains some collection of parcels. At each point in the sampling region a 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 mesh points that will receive a contribution from a given parcel partially or wholly contained  Chapter 7. The System^  85  in the rectangular region will satisfy:  iPi; — P1 < 7- 2 0 < i < p—1  0<j<q—1 where r is the radius of the parcel P. At each mesh point the contribution of the parcel to its thickness will be accumulated. The thickness at a point is determined by the distribution 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— V  irr 2  If the distribution of petroleum over the parcel was conical:  h = h o (r — AO ho = 3  V irr2  Or = iPi.i — Pi After the parcels have been traversed each location of the mesh has the accumulated slick thickness recorded. This information is encoded into the p x q image. The thickness is encoded into a single byte, using a logarithmic encoding scheme. This schema was developed to preserve the large changes in thickness over the integrated slicks (Figure 7.6). The sampling rectangle by default is determined by the intersection of the plane z=0 and the four sided pyramid emanating for the viewers eye position (see figure 7.5). The size 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 b meters where b > 0 and a < 0. The maximum thickness difference is b — a orders of  Chapter 7. The System  ^  86  Viewpoint  Figure 7.5: Viewpoint and Sampling  ^  q  Region Figure 7.6: Sampling of p x q Grid magnitude. If the thickness t = 0, then 0 is placed in the image. If t > 0 then the pixel is encoded into an 8 bit byte as Wlog io (t10 a) where W = (2 8 — 1)/(b — a). The -  thickness can be recovered by reversing the operation. There will be some loss minor loss in precision, but the range of thicknesses will be preserved. For the system, b = 1 and a = —5, on the assumption that the slick will never be thicker than a meter, nor thinner that 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 the surface and as a visual representation of the slick in the form of a grey scale image where intensity is proportional to thickness. The boundaries of the radial parcels can in some cases be quite pronounced. To remove this effect a low pass filter can be passed over the thickness map before encoding. 7.4.2 Interpolation of the Slick Between Frames Given two frames of an oil slick sampled at different times, it is necessary to be able to interpolate an intermediate state. This can be done either by performing a metamorphosis  Chapter 7. The System^  87  between two images, or by direct interpolation of the script files from which the images were generated. To attempt a metamorphosis from one thickness map to another, extraction of the boundary 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 by frame basis. To accomplish this, methods of surface reconstruction were considered. Contours of equal value (in this case slick thickness) are determined from two thickness map separated by some time step, and an attempt is made to match them. From this match, an intermediate contour can be determined. Repeating this process over some range of contours on both thickness maps, an intermediate state can be interpolated and reassembled. These methods are discussed in [Slo87]. The main problem is the matching of points from different contours. The problem reduces to a search on a class of graphs called toroidal 4 . Contours were extracted from the slick image. The zero contour (the boundary between no pollution and some pollution) was extracted. Each pixel in the image was set to 0 if it was void of pollution, and 1 if any trace of the pollutant existed. This provides a noise 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 or edge detection. Contour tracing, such as discussed in [Pav82], involves determining the chain of pixels which form the boundary of a patch of given intensity in an image. Edge detection 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 are dealing with a binary image, the two operations can be seen as more or less equivalent. 4 A toroidal graph consists of two collections of vertices Ci and C2 where the edges between vertices within a given set C1 or C2 form only cycles. These cycles form the contours. Edges between the sets C1 and C2 form the match between the contours  Chapter 7. The System^  88  The Laboratory for Computational Intelligence has a developed set of image processing tools called VISTA. This software was used as an alternative to developing software to perform image processing tasks. The grey scale image was generalized into a binary image, with non-zero regions denoting regions contaminated by oil. The binary image was converted into the LCI format, and the components of VISTA used to extract the boundaries. As the image was binary and synthetic, no noise was present. As a result low pass filtering was not required. The edges were found by application of a Canny edge detector [Tor82], followed by a routine which linked up edge pixels. These links were converted into an generic ASCII format. The recovery of edges was quite good, with very few breaks occurring in the recovered boundaries. Unfortunately, attempts to enforce correspondence on a frame by frame basis could not be accomplished without human intervention. This was due to the fact that the number of components from frame to frame was not fixed. Applying the methods of surface reconstruction from collections of contours with this property tends to result in self intersecting surfaces. This results in a many to many relationship from N to M contours, where N M. Further, as the slick was the result of the reinitegration of a set of circular parcels, each one could potentially represent an independent slick (in practice very few parcels were found to be distinct from a main body). Therefore, effort was redirected to a direct interpolation of the slick scripts. Over time, the slick scripts will consist of seaborne and grounded parcels. No interpolation is required with grounded parcels, as their positions do not change with time. The unique identifier 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 the environment, 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 is  Chapter 7. The System^  89  the identifier of P, then if a parcel 13 ' exists on Lk+1 such that .13.„ is its identifier, then  P and .13 ' match. The position, volume and radius of the interpolated parcel .P" can be interpolated at time t k II where tk < t k " < tk+1 from: tk+1 — 4+1 - tk sPk + (1 — s)Pk +1  V " = s Vk + ( 1 — s )Vk + i r = srk + (1 — s)rk +1 If no point with the identifier Pu exists, then the parcel vanished in going from frame k to frame k + 1. This parcel will be interpolated at t k " from: S=  tk-} i — t k' tk+1^tk  13 " = Pk V = SVk  r = srk This process traverses both lists, removing matches from both lists. Parcels in Lk which do not match in Lk +i are also removed. When Lk is exhausted, parcels may remain on Lk +1 . These parcels are the result of venting from the oil source which occurred between  tk and 4+1, and are not interpolated. This interpolation is demonstrated in figure 7.7. 7.4.3 Other Rendering Types The tidal model presented an opportunity to render in detail certain properties of the surface currents and tidal elevation. These are specified with the particular rendering options in the script files. The tidal elevation over a domain can be rendered, using colour as the visual cue. In a manner similar to the thickness map, samples are laid over the extent of the domain  Chapter 7. The System^  90  t=o  t = t`  t=1 3  Figure 7.7: Interpolation of an Oil Parcel Over Time with each sample point corresponding to a pixel in the resulting image. A tidal elevation is sampled at each location. This is mapped into a colour by using the HSV colour model. The colour subspace for V = 1 corresponds to a circles with the hues along the circumference. A specific colours is specified by its angle 7 , with complementary colours located 180 degrees opposite of each other. The value of S, or saturation, is a ratio which is 1 at the circumference of the circle and 0 at the center. Varying S has the effect of interpolating 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 transformation would not be as straight-forward as HSV. If the tidal elevation is 77 then the colouring is based upon:  H = flifi7>0 H = 180 — ,8 if ri < 0 Ti S =^ qmax  1n the HSV model, V=1 corresponds to maximum intensity Alternately defined as hue 7 With red located at 0 degrees 5  6  Chapter 7. The System^  91  17 = 1 where # is the current hue. Using a single hue, or the more familiar grey scale, to represent a thickness or range is a relatively intuitive approach. Using the HSV system allows 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 height of the tide, and black or no intensity to represent the minimum tidal elevation, then midtide or zero tidal elevation would be represented by the middle grey value. This cannot be readily discerned. A function to map HSV to RGB completes the colouring. The effect 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 is approximately mean tide. In an animation, this shows the rise and fall of the tide if an adequate time step and number of frames is used. However, as the change in elevation is very smooth over time, no local spatial variation in elevation on a frame by frame basis is discernable. While little variation in elevation is evident, this is not true with the horizontal currents (figures 7.8 and 7.9). For the given domain a sample of the horizontal flow at each point corresponding to a pixels is determined. The maximum flow in the region must be determined or estimated for the region. At present a grey scale value is determined at a location based upon: pp  = ( 28 1) Tfmax — Tf Tim=  The resulting image displays areas of large current flow as bright white, with areas of relative calm as a dark shade of grey, or black in the case of no flow at all. In the context of an animation, the effects of channeling and shallow water can be seen clearly. This will be expanded further in chapter 8. The ability to output a representation of an oil spill into the Optik ray tracing system  Chapter 7. The System  ^  92  Figure 7.8: Current Flow in Howe ^Figure 7.9: Current Flow in Vancouver Sound Harbour [Ama90] was provided. For every parcel in a slick, a corresponding Optik object is created. At present these objects are either particles or spheres. The object is translated to its model coordinate position, and the radius of the object 8 is set to the parcel's radius (assuming the parcel is a sphere). An oil spill is represented by a collection of intersecting spheres or disks. A sphere can be scaled into an ellipsoid with its two major axis 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, with  z used to differentiate between time slices. 7.4.4 Texture Mapping an Oil Spill One of the intended purposes of generating the thickness map in the form of an image was to use it as a texture map in a pre-existing rendering system, such as Optik. The slick at a given point will produce a colour dependent on its thickness, state of degradation, 8 With  respect to a particle or sphere.  Chapter 7. The System^  93  and the orientation of the surface at that point. Ideally for a given spill scenario, an image of the scene should be rendered which shows both the position and the general optical properties of the oil slick. Texture mapping is the geometric mapping of a two dimensional 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 a region, and the surface is the ocean's surface and the local foreshore which correspond to the spatial location of the slick. In our case, the surface is a plane. However, given an approximation of the surface of the ocean resulting from wave action, such as developed in [Fou86], the scene can be made more consistent with reality. The oil slick is not a texture in the strict sense. When the slick is mapped onto the surface, at each point where oil is present we would like the corresponding surface colour in the scene to be affected. However, at each point where there is no contamination, then the radiance of the corresponding region should remain unchanged. Most graphic systems support texture maps in some form. Unfortunately, two problems surface on closer examination. These problems are the registration between the texture map and the underlying surface, and compositing the thickness map onto the surface. In Vertigo, Wavefront and Alias, it was not possible to bring the three dimensional model into the system without performing scaling, rotation and translation operations. As a result, the series of transformations that 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 locations generated in the oil slick modelling system. The Optik rendering system had no such constraint. This texture mapping will be discussed further in chapter 8. The various models and custodial functions which comprise the oil spill simulation were introduced in this chapter. Integrating the components developed during modelling yields a system which provides a means for defining and visualizing a marine oil spill.  Chapter 8  Experiments and Examples  8.1 Standard Slick Formulation For the testing stage, a fixed set of spatial, temporal, and physical parameters for the simulation are set as a fixed model, from which the performance of adjusted models will be measured. The parameters for this fixed model are listed in table 8.1. The resulting oil spill is shown in figure 8.1. Parameter  t P l'h a crT  P. Pw  At Vo fr  Value July 17, 1992, 6 am Fixed at test time 5.0 0.025 1.00 0.88 1.0275 Time Step 106 10  Comments Date of Sampling for Interpolated Wind Data Spatial Location of Source Horizontal Dispersion Coefficient Wind Coupling Coefficient Tide Coefficient Specific Gravity of Test Petroleum Specific Gravity of Ocean Water Adaptive or Fixed Liters Liters per second  Table 8.1: Default Slick Model  8.2 Testing for Sensitivity to Space and Time The first collection of tests will determine the sensitivity of the simulation to initial conditions in space and time. This will be a good measure of the quality of the overall integrated simulation. For a fixed set of slick parameters, executed at time t and position 94  Chapter 8. Experiments and Examples ^  95  Figure 8.1: The Standard Oil Spill  P, the simulation after some period of time TD will produce a result. Visually, it will be measured by the final distribution achieved by the parcels within the spatial domain. The spread, size, orientation, number of components, and the amount of shoreline affected can all be measured qualitatively by visual inspection of the rendered scene. If the slick 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 outcome of the simulation for time TD should be similar if the system is not chaotic. If the result is unacceptably different, then the simulation is exhibiting chaotic behavior, or a sensitivity to initial spatial and temporal conditions. If the system does not have this chaotic property, then small changes in position and time will result in a small changes in the distribution of the particles over time and space.  Chapter 8. Experiments and Examples ^  96  The means of determining this similarity can be qualitative or quantitative. The visual inspection of the two results is the most direct route. An intervening human in most cases will be able to determine if, in general, the simulation is behaving in a stable manner, and that patterns occurring in one simulation occur in the other. Quantitatively, there are various means by which the stability of the system could be determined. The analysis of the system on a parcel by parcel basis would be difficult. The random processes involved in the dispersion and grounding processes remove the possibility of retaining sensible correspondence between parcels from different simulations, despite the fixed spill parameters. Rather, generalized statistics would need to be derived, such as the center of mass of the respective parcel clusters, the difference in slick area, and the volume remaining 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 two simulations could be performed, and the results accepted or rejected based on comparison of these statistics.  Figure 8.2: Position 200m West ^Figure 8.3: Position 200m East  Chapter 8. Experiments and Examples^  97  A collection of simulations were run on the system and it was visually confirmed that small 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 position and the simulation executed for 120 minutes. A change in the order of 500 meters usually resulted 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 vary spatially. The starting time was varied between 1 and 60 minutes, with changes over 25 minutes usually producing vastly different results. This would be expected, as the tide is changing with time, as are the winds. Figures 8.2, 8.3, 8.4, and 8.5 show the distribution patterns achieved by varying space and time.  Figure 8.4: Time Set 15m Back  Figure 8.5: Time Set 15m Ahead  8.3 Testing for Spreading and Advection The spreading model was tested by altering the dispersion coefficient -yh and alternately turning the physical spreading off and on. The dispersion coefficient was varied from  Chapter 8. Experiments and Examples ^  98  0.5 to 25 m 2 s -1 . The larger the dispersion coefficient, the greater the diffusion of the particles over the surface of the ocean. The results of these tests are shown in figures 8.8 and 8.9.  Figure 8.6: No Turbulent Diffusion Figure 8.7: No Physical Spreading The physical spreading was tested by executing an instantaneous spill. The rate of flow was set to a value sufficiently high as to allow for the entire volume to vent in one time 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 seen in figures 8.10, 8.6, and 8.7, which show the spreading resulting from just physical, just turbulent, 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 of tests were run with values ranging from 0.001 and 0.10. The influence of this parameter on the spread and advection patterns is very apparent, with the distance induced by drift proportional to this value. The aT parameter, while having no physical basis, was  Chapter 8. Experiments and Examples ^  Figure 8.8: -yh = 1.0m 2 s -1  99  Figure 8.9: yh = 25.0m 2 8 -1  intended to allow for varying the magnitude of the tidal velocity field. Spread patterns exclusively due to the wind, or to the tide, could then be observed. The results of this experimentation is shown in figures 8.11, and 8.12. Visualizing only the effects of advection of the particles was achieved by switching off the physical spread, and setting -yh to be very small. This would be analogous to spilling a collection of buoys rather than oil. It shows how the system can be adjusted to serve other purposes, such as tracking a drifting vessel. A result of this type is shown in figure 8.13. It was originally stated that the time step used in the simulation was adaptive. However, 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 adaptive time step, and a large fixed time step of 300 seconds. The results of this are shown in figure 8.14 and 8.15. They should be compared with figure 8.1. These images show how the time step greatly affects the result.  Chapter 8. Experiments and Examples ^  100  Figure 8.10: No Wind or Tide, Both Spreading Types On  8.4 Testing the Oil Source The effect of the density of the oil was tested. The difference was measured by visual inspection of the thickness maps generated by the extended visualizer. Key-frames of the simulation were captured involving petroleum with different specific gravities. The effects of petroleum density are manifested with respect to the amount of radial spread occurring at a particle level, and to the overall physical spreading of the oil body. Figures 8.16 and 8.17 show the effect of marked differences in the densities introduced to the model. One feature that is apparent is that the scale of diffusion of the high density slick is far too large (figure 8.16). The fact that many of the spread discs do not overlap shows  Chapter 8. Experiments and Examples^  101  Figure 8.11: a = 0.08 Figure 8.12: a j = 0.10 that the turbulent diffusion process is overwhelming the spreading process far to early in the simulation. This can be corrected by making -yh depend on the density of the petroleum.  8.5 Interpolation Between Frames Two methods were applied to interpolating between two key-frames of slick simulation data. The first method involved direct interpolation between the two thickness maps. This involves generating a slice of the thickness map by threshholding on thickness. This results in a high contrast binary image, which has an edge detector applied. All edges which 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. The result of interpolating a slick is shown in figure 8.20.  Chapter 8. Experiments and Examples ^  102  Figure 8.13: No Effects from Spreading or Diffusion 8.6 Altering Parameters Based on Tests The parameter l'it was determined by a function of both sea state and oil density. This was discussed in chapter 6. The magnitude of this parameter depends on the turbulence and the properties of the oil. -yh was parameterized over sea state and density by applying: 7h = (S + 1.0)L '10 where S is the sea state reported as a Beaufort number, P o is the density of the source petroleum, and -yo is the horizontal dispersion coefficient estimated for an oil with the same density as water on a flat calm sea under no influence from currents. While this formulation is empirical, it is constructed to be consistent with the processes in action. The term -yo/p o grows with decreasing density. The term S +1.0 applies a linear increase  Chapter 8. Experiments and Examples^  103  Figure 8.14: Adaptive Time Step Figure 8.15: A t = 300 Seconds in diffusion based on the surface turbulence, as estimated by the sea state model. The horizontal diffusion grows most rapidly when a low density petroleum is introduced into a heavy sea. Using the values for the standard slick set, the result of this experiment are shown in figure 8.21.  8.7 Texture Mapping a Slick into the 3 Dimensional Model -  A thickness map was converted into the Wavefront system to render a consistent image of an oil spill in the test domain. A second texture map to approximate the surface properties 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 rectangular colour map corresponding to the modelled region. Alternately, if satellite imagery was available, it could be used as the texture map. The 3-dimensional model was ported into Wavefront, and then both texture maps were registered to this model. The image of the oil distribution was coloured using a brownish-green shade, modelled after descriptions  Chapter 8. Experiments and Examples^  104  Figure 8.16: Density p o = 0.75 Figure 8.17: Density p o = 0.99 of 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 change in properties. The result of this work is shown in 8.22.  Chapter 8. Experiments and Examples^  Figure 8.18: Slick Thickness Map ^Figure 8.19: Resulting Edges  Figure 8.20: Results of Interpolation  105  Chapter 8. Experiments and Examples  Figure 8.21: Slick Produced Using Adaptive Spread Algorithm  106  Chapter 8. Experiments and Examples  ^  Figure 8.22: An Oil Slick in Vancouver Harbour  107  Chapter 9  Conclusion  9.1 Future Work A system of this ilk can rarely be termed complete. There is always one more development, or level of refinement which could be implemented. This chapter will briefly discuss the future of systems of this nature, and the short and long term changes that could be made to improve this particular system. This system operates strictly on the ocean's surface. While slick pollution is initially a surface problem, the chemical degradation extends the pollution into the atmosphere and down into the water column A system of this nature must eventually address this third dimension, not only with respect to the oil source, but also in relation to the environmental processes. Subsurface currents induced by wind and tides, the surface wind current, and the topography of the ocean surface, could be modelled, increasing the accuracy of the system. This extension would not only increase the consistency of the model, but would change the nature of the visualization that could be achieved. Volume rendering of the resulting subsurface and atmospheric plumes could be considered. The model would be improved by a more complete description and modelling of the weather process in action. This would require the collection of extensive data on a number of major petroleum crudes and distillates, and their chemical and physical properties. The inclusion of detailed weathering models would degrade the performance of the simulation, and could potentially render it non-interactive. It should also be  108  Chapter 9. Conclusion^  109  noted that weather systems and fluid dynamics remain a open area of research. Even the leading edge of this science would be insufficient to completely describe all weather related processes. The greatest weakness of this system is its lack of concrete validation. This is primarily a result of the lack of hard data on oil spills in the domain. While spill occurrence data is available, it does not provide any information on the spatial distribution of the slick over time, nor the breakup time of the spill. A possibility for validating the system would be to acquire the required data for modelling in a domain which does have detailed spill data from various events. Generation of the coastal terrain, tides and winds of the area at 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 the actual slick event. In terms of software, the system should be redeveloped to use an system independent interface. This could be done using the X windows paradigm [Nye90], or by developing a simple and portable interface. The component software could use a re-evaluation in terms of both general functionality, and time efficiency. With a careful redesign, the system could be redeveloped to run on most existing platforms (assuming adequate graphics capabilities). Further, it could be extended into a general library of particle functions and visualization. The generalization of the system into a generic particle system would require the addition of the third dimension in terms of both advective fields and topography. Vertical processes would then require more attention. With a system of this nature visualization of the flows under the ocean's surface and into the atmosphere could be realized. 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 should be developed. Barriers, such as oil spill booms and skimmers could be introduced to block the flow of the spill. Chemical dispersants, which increase the diffusion rate of  Chapter 9. Conclusion^  110  the spill, could be included by altering the existing chemical degradation models. The addition of these features would allow for comparison of the amount of damage caused by an oil spill with or without contemporary oil spill fighting equipment. This could be used as a measure of the relative effect of these techniques. Further, theoretical spill fighting methods could be tested within the framework of this system.  9.2 Conclusion This thesis set out to develop a consistent system for modelling and visualizing the dynamic process of an oil spill in the marine environment. The criteria for development was to produce a system which was interactive, both in terms of human interaction, and in the performance of the system, and possessed a generic mechanism for the input of source data, thus achieving domain independence. The system has achieved consistency at the level dictated by the assumptions and generalizations under which the models were developed. As has been stressed in previous discussions, the issue of modelling both the oil spill and the environment requires some generalizations. Attempts have been made to include all features that directly affect the process of oil spill pollution. While the inherent complexity of a system is not a valid reason for choosing a model, or for not including the process, many of the components of the system have been generalized to a point where further condensation would begin to remove validity. The visualization developed provide a base level tool for determining the spatial distributions of oil spill pollution over time. Additional functionality, in the form of flow visualization and in-depth rendering of particular environmental processes, allows the user of the system to analyze visually not only the resulting spill, but the major influences on the spill.  Chapter 9. Conclusion^  111  Over the last decade environmental awareness has spread and became a mainstream topic from the coffee shop to the board room. Whether or not this apparent concern will lead to a more sensible, rational and forward looking approach to development and the usage of ever-diminishing resources in the future is a question for philosophers and historians. Be that as it may, this system has attempted in a small way to contribute to the body of knowledge that can be applied to arriving at this lofty goal.  Appendix A  Appendix  A.0.1 Determining the Wind Drift To determine the direction and magnitude of the wind-induced drift current, the solution of the equations of motion with friction present must be solved. The system presented is due to [Pon83]. Two terms must be considered: the wind as a result of the horizontal pressure gradient, and the vertical friction. Note that this derivation is intended for the northern hemisphere. Minor adjustments are required for the southern hemisphere. The components are solved separately and added together, yielding: f v = f (v. + v E ) where f = 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 velocity vE = the Ekman velocity associated with vertical shear friction This separation is possible as the equations are linear. The wind induced drift current is determined by the superposition of these terms. Assumptions made include a homogeneous water body, a flat surface, an infinitely deep ocean avoiding bottom friction, and infinite horizontal extent to avoid associated boundary problems. Given this, and assuming the wind to be blowing in the y direction the Ekman equations becomes: r E z)eff u = lio cos(i + ±  z  112  Appendix A. Appendix^  113  i .-z)e* z v = Vo sin(i + f  where Vo 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 counterclockwise to the wind direction. To obtain a numerical relationship between the surface current Vo , observational experiments determined: Tn = PaCD  I W I, the wind stress magnitude,  where  p a 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 _ W  0.0127  VIC01)  yielding  DE =  4.3W  N/sin101)  From this if the wind in ms -1 is known, Vo and DE can be found at any depth, and  v,,  in particular on the surface. At the 45th parallel the ratio CF = i the wind coupling coefficient, has been estimated at 0.015.  sometimes called  Appendix A. Appendix^  114  A.0.2 System Specifics The system was developed under ANSI C on various Silicon Graphics platforms. The code is split into three groups, user interfaces, graphics and simulation-utility. The user interface was developed using the FORMS interfaces library. The graphics routines uses the Silicon Graphics GL embedded programming language for graphics initialization and primitives. The simulation-utility group was developed primarily by the author, with some 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 to find 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. The  database program does not require any support files, but will require that permission is set for read/write for the graphics group. To extract a dataset type database. This will bring up a graphic display of the domain, and an interactive panel. Specify the path and name of the dataset, and extract. When the extraction is complete the user can inspect the resultant directory. It will contain files with the names g100b, tide.; tide.u, tide.v which are respectively the terrain and tidal representation. The directory wind contains data files describing the wind. To run the oil simulator, the user must have a copy of the file parameters. dat. This file contains the defaults of the physical and simulation parameters used by the system. The default parameters.dat contains: 1. wind coupling coefficient, default 0.03 2. tide coupling coefficient, default 1.00 3. horizontal dispersion coefficient, default 1.0 4. vertical dispersion coefficient, default 0.001  Appendix A. Appendix^  5. tidal component flag 1, default 1 6. tidal component flag 2, default 1 7. tidal component flag 3, default 1 8. tidal component flag 4, default 1 9. tidal component flag 5, default 1 10. tidal component flag 6, default 1 11. tidal component flag 7, default 1 12. tidal component flag 8, default 1 13. tidal component flag 9, default 1 14. tidal component flag 10, default 1 15. tidal component flag 11, default 1 16. tidal component flag 12, default 1 17. tidal component flag 13, default 1 18. tidal component flag 14, default 1 19. hue for HSV, default 75.0 20. upper range for elevation, default 10.0 21. lower range for elevation, default 0.0 22. denominator and numerator of sea state coefficient, default 10,1  115  Appendix A. Appendix^  116  23. fixed drop point flag, default 0 24. 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 system type oilsim. The input panels are self describing.  Bibliography  [A1h75] S. W. Ahlstrom, A Mathematical Model for Predicting the Transport of Oil Slicks In Marine Waters, Battelle Pacific Northwest Laboratories, Richland, Washington, 1975. [A1R89] A. H. Al-Rabeh, H. M. Cekirge, and N. Gunay, A Stochastic Simulation of Oil Spill 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, Optik Users Manual, Version 2.1, 1990 [ASP83] American Society of Photogrammetry, Manual of Remote Sensing, 2nd ed., Falls Church, Va., 1983 [Ara82] K. P. Aravamudan, P. Raj, J. Ostlund, E. Newman, and W. Tucker, Break up of 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 Atmospheric 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, recent developments, and major issues, Simulation, Vol. 49, Number 3, pp 109-116, 1987. [Ben77] A. P. Bentz, Effects of Weathering on Oil, Workshop on Pattern Recognition Applied to Oil Identification, IEEE Cat. 76, 1977. [Bow83] K. F. Bowden, Physical Oceanography of Coastal Waters, Ellis Horwood Limited, 1983. [Bra92] P. Bragatto, N. Mazzino, and P. Palamidese, Animated Visualization of Scalar Fields, Second Eurographics Workshop on Animation and Simulation, Vienna, pp 115-127, 1992. [Bre77] J. A. Brewer, and D. C. Anderson, Visual Interaction with Overhauser Curves and Surfaces, Quarterly Report of Siggraph, Vol. 11, Number 2, pp 132-137, 1977.  117  Bibliography^  118  [Bri91] M. Briscolini, and P. Santangelo, Animation of computer simulations of twodimensional turbulence and three-dimensional flows, IBM Journal of Research and Development, Vol. 35, Number 1/2, pp 119-137, March 1991. [Bro90] F. P. Brooks Jr., Ming Ouh-Young and James J. Batter, Project GROPE - Haptic Displays for Scientific Visualizations, Proceedings ACM Siggraph, pp 177-185, 1990. [Buc90] J. Buchanan., MxMatrix Math Routines, Imager Laboratory, University of British Columbia, 1990. [Cek88] H. M. Cekirge, A. H. Al-Rabeh, and N. Gunay, Determining the Wind Driven Surface 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 Oceanographic Numerical Modelling, D. Riedel Publishing Company, 1986. [Cre88] P. B. Crean, T. S. Murty, and J. A. Stronach, Mathematical Modelling of Tides and Estuarine Circulation: The Coastal Seas of Southern British Columbia and Washington State, Springer-Verlag, Berlin, 1988. [EMR88] Digital Topographic Map Data, Topographical Mapping Division, Department of 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 of Oil, 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, Proceedings of the Joint Conferences on the Prevention and Control of Oil Spills, American Petroleum 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, pp 191-199, 1989. [GL90] Graphics Library Programming Guide, Silicon Graphics Inc., Document Number 007-1210-010, 1989 [Gun87] E.R. Gundlach, Oil-Holding Capacities and Removal Coefficients for Different Shoreline 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 Model for Scientific Visualization Systems, IEEE Computer Society Press Tutorial, Los Alamitos, pp 74-92, 1990. [Hec86] P. S. Heckbert, Survey of Texture Mapping, Proceedings Graphics Interface, pp 56-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, Academic Press, San Diego, 1979. [Hor72] B. Horstein, The Appearance and Visibility of Thin Oil Films on Water, EPAR2-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, A Mathematical Model for Transport and Mixing of Oil Slicks in Rivers, Proceedings of 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 the Marine Environment, Ann Arbor Science, 1980. [Kaw70] Kawahara and Ballinger, Characterization of Oil Slicks on Surface Waters, Industrial 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 residual currents and pollutant transport in the Arabian Gulf, Applied Mathematical Modelling, Vol. 12, pp 379-390, 1988. [Mac77] D. Mackay, and P. J. Leinonen, Mathematical Models of the Behavior of Oil Spills on Water with Natural and Chemical Dispersion, Economic and Technical Revue Report, University of Toronto, 1977. [Mac78] Mackay and Wilks, A Statistical Analysis of Oil Spills in Canada, Research Report No. 50, Government of Canada. June, 1978. [Mac80] D. Mackay, S. Paterson, and K. Trudel, A Mathematical Model of Oil Spill Behavior, Environmental Protection Service, Fisheries and Environment Canada, 1980. [Mas85] C. F. Mass, and D. P. Dempsey, A One-Level, Mesoscale Model for Diagnosing Surface Winds in Mountainous and Coastal Regions, Monthly Weather Review, July 1985. [Nae79] A. Naess, The Mixing of Oil Spills Into the Sea by Breaking Waves, Journal of Petroleum Technology, pp 1113-1122, 1980. [Nie89] G. M. Nielson, Visualization in Scientific Computing, Computer, pp 10-25, August 1989. [Nih84] C. J. Nihoul, A Non - Linear Mathematical Model for the Transport and Spreading of Oil Spills, Ecological Modelling, Number 22, pp 325-339, 1984. [Nye90] A. Nye, Volume 1 : Xlib Programming Manual for Version II, O'Reilly and Associates Inc, Sebastopol, California, April 1990. [Obr85] J. O'Brien, Advanced Physical Oceanographic Numerical Models, NATO ASI Series, 1985. [Ove91] M.H. Overmars, Forms : A Graphical User Interface Toolkit for Silicon Graphics Workstations, 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 to a Risk Estimation Model, Proceedings International Symposium Natural and Man Made Hazards, D. Reidel Publishing Company, 1986.  Bibliography^  121  [Pes91] R. L. Peskins, S. S. Walther, A. M. Froncioni, and T. I. Boubez, Interactive Quantitative 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 Over South Florida, Monthly Weather Review, Vol. 102, February, pp 115-134, 1974. [Pon83] S. Pond, and G. L. Pickard, Introductory Dynamical Oceanography, Pergamon Press, 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 Publishing Company, 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 Initial Results, Graphics Interface 87, 1987. [Ste88] D. G. Steyn, and I. G. McKendry, Quantitative and Qualitative Evaluation of a Three-Dimensional Mesoscale Numerical Model Simulation of a Sea Breeze in Complex 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?, Atmospheric 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 Fate Modelling, 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 Pattern Matching 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 Service, Part 1: Theory and Model Evaluation, Atmosphere and Ocean, Vol. 22, pp 93 - 108, 1987. [Ven90] S. Venkatesh, Model Simulations of the Drift and Spread of the Exxon Valdez Oil Spill, Atmosphere and Ocean, Vol. 28, pp 90 - 105, 1990. [Viz73] K. N. Vizy, Detecting and Monitoring Oil Slicks with Aerial Photographs, Eastman Kodak Company, 1973. [Yap89] P. D. Yapa, R. J. Thomas Jr., R. S. Rutherford, and H. T. Shen, Microcomputer Model for Oil Spill Simulation, Journal of Computing in Civil Engineering, Number 3, 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