UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Investigating heinrich events : a continuum model of iceberg drift and sedimentation La Prairie, Douglas I. 2000

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

Item Metadata


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

Full Text

I N V E S T I G A T I N G H E I N R I C H E V E N T S : A C O N T I N U U M M O D E L OF I C E B E R G D R I F T A N D S E D I M E N T A T I O N By Douglas I. La Prairie B. Eng. (Naval Architecture) Memorial University of Newfoundland A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF S C I E N C E in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October 1999 © Douglas I. La Prairie, 1999 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 Earth and Ocean Sciences The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract Records from the North Atlantic reveal high deposition rates of ice-rafted sediments from the Laurentide ice sheet during at least six distinct periods. Known as Heinrich events, these episodes of rapid deposition are attributed to massive iceberg calving promoted by unstable behavior of the ice sheet. The exact mechanisms causing these occurrences of high iceberg flux are not well understood. As a means of elucidating and exploring the freshwater and sediment fluxes associated with Heinrich events, a computer model has been developed which attempts to simulate iceberg drift during these periods of massive calving. This enables reconstruction of the marine sedimentation pattern and provenance-labeled marine core stratigraphy. Using finite-difference methods over a dis-cretized problem domain, the model treats the armada of icebergs as a continuum rather than tracking an ensemble of perhaps thousands of individual trajectories. In this manner the drift of the model icebergs in a multi-layer ocean is governed by advective, diffusive and Coriolis forces as well as being subject to melting. The stochastic nature of iceberg drift is dealt with using a tensorial diffusion coefficient to approximate the effects of heterogeneity, such as different shapes, which strongly influence iceberg drift. n Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement x List of Symbols xi 1 Glaciomarine Mechanics and Heinrich Events 1 1.1 Introduction 1 1.2 Iceberg Sources and Production 2 1.3 Ocean Drift Patterns 4 1.4 Drift Forces 6 1.5 Deterioration 7 1.6 A Curious Sediment Record 8 1.7 Thermohaline Circulation 12 1.8 Global Links 14 2 A Continuum Model of Iceberg Drift: Application to Heinrich Events 16 2.1 Iceberg Drift Models: State of the Art to Present 16 2.2 Physics of Drifting 19 iii 2.3 A New Approach: Icebergs as a Continuum 20 2.4 A New Model 22 2.4.1 Iceberg Calving Sources 23 2.4.2 Ocean and Atmosphere Fields 25 2.4.3 Drift Dynamics 28 2.4.4 Melting 28 2.4.5 Size Coupling 29 2.4.6 Sedimentation 29 2.4.7 Sea-Ice Gating 29 2.5 Model Runs 30 2.5.1 L G M 1,000 Year Run 32 2.5.2 LGM/Heinrich Event 1,000 Year Run 34 2.5.3 Modern 1,000 Year Run 40 2.6 Drill Core Prediction 42 2.7 Conclusions 48 2.8 Future Developments 49 References 50 Appendices 54 A Formulation of the Governing Equations 54 A . l The Continuous Model 54 A.1.1 Volume Density of Icebergs 54 A. 1.2 Volume Production by Sources 55 A.1.3 Volume Diffusion of Icebergs 55 A.1.4 Reynolds Transport: Volume Advection of Icebergs 56 iv A.1.5 Melting of Icebergs 57 A. 2 The Discretized Model 57 A.2.1 Local Forces and Acceleration: Determining Velocity Fields . . . . 61 A.2.2 Bathymetry 66 A. 2.3 Two-Dimensional Coordinates 66 B Stochastic Diffusion 69 B. l Theory 69 B. l . l Velocity Dependence 70 B.1.2 Coriolis Dependence 70 B. l .3 Sea-Ice Gating 71 B.1.4 Formulation of Diffusion Coefficients 71 B.2 Experimental Determination of the Diffusion Coefficients 73 B.2.1 Control Condition 73 B.2.2 Analyses 74 B.2.3 Results 79 C Iceberg Size Distribution at Sources 81 D Melting, Size Coupling and Sediment Deposition 86 D . l Melting 86 D . l . l Solar Radiation 86 D.1.2 Vertical Buoyant Convection 87 D.l .3 Air Convection 88 D.1.4 Wave Action 88 D.1.5 Total Melt and Numerical Formulation 89 D.2 Size Coupling 90 v D. 3 Sediment Deposition E Time Stepping and Sparse System Solution E. l Temporal Discretization E.2 Sparse Matrix Solution F Main Program and Routine Flow Charts v i List of Tables 2.1 Model run input parameters 30 2.2 Provenance dependent input parameters 32 B . l Iceberg control parameters 74 B.2 Atmosphere and ocean control parameters 74 B . 3 Variable drift parameters 75 C. l Greenlandic source iceberg flux distribution by size 83 C.2 Icelandic source iceberg flux distribution by size 84 vn List of Figures 1.1 Iceberg sources 3 1.2 Iceberg calving 4 1.3 Greenland S 180 record 11 1.4 Schematic Bond and Dansgaard-Oeschger cycles 11 1.5 Schematic thermohaline circulation 13 2.1 Ice sheet model input 24 2.2 Ocean model input 26 2.3 Ocean model input 27 2.4 Ocean model input 28 2.5 Ocean model input 31 2.6 L G M ice distribution 33 2.7 L G M meltwater distribution 33 2.8 Heinrich Event ice distribution 36 2.9 Heinrich Event meltwater distribution 36 2.10 Heinrich Event 1,000 year meltwater distribution by provenance 37 2.11 Heinrich Event 5 year meltwater distribution by provenance 38 2.12 Heinrich Event melwater flux 39 2.13 Modern ocean meltwater distribution by source 41 2.14 Modern ocean meltwater distribution by source 41 2.15 International Ice Patrol iceberg sightings 1994-1998 42 2.16 North Atlantic sediment distribution 45 vm 2.17 Drill core locations 45 2.18 Synthetic marine SU 90-08 46 2.19 Marine core SU 90-08 46 2.20 Synthetic marine core Me 69-19 47 2.21 Marine core Me 69-19 47 A . l Discretized North Atlantic problem domain 58 A.2 Five-point finite-difference cell scheme on a spherical grid 67 A. 3 Spherical grid: two-dimensional equivalent 68 B. l Drift variation after three days 77 B.2 Example drift trajectories 78 B. 3 Diffusion coefficient variabilty with perturbation 80 C. l Greenland source iceberg size distribution 82 C. 2 Iceland source iceberg size distribution 84 D. l Sediment content by size and provenance 92 F . l Main program 98 F.2 Iceberg velocity calculation 99 F.3 Generic routine for depth dependent parameters 100 ix Acknowledgement Foremost, I extend my gratitude to Garry Clarke for being a truly ideal supervisor. Garry was accommodating, more than could reasonably be expected, in allowing me to work independently according to a sometimes awkward schedule. Thanks are also in order to Dave Hildes, Jeff Kavanaugh, Gwenn Flowers and Shawn Marshall for their generosity with time and advice when needed.. This work was supported by funding from the University of British Columbia and the Natural Sciences and Engineering Research Council of Canada. x List of Symbols A general area, m 2 AA iceberg area above water projected in flow direction, m 2 Aw iceberg area below water projected in flow direction, m 2 C velocity dependence coefficient CDA drag coefficient in air CDW drag coefficient in water D tensorial diffusion coefficient, m m _ 1 s _ 1 = s - 1 di ocean layer I depth, m DA along-stream diffusion coefficient, m m - 1 s _ 1 = s _ 1 Di iceberg draft, m DT cross-stream diffusion coefficient, m m _ 1 s _ 1 = s _ 1 / Coriolis frequency, s - 1 fw wave frequency, s - 1 FA air drag force, N Fw water drag force, N J diffusion flux per unit area, m 3 m~ 2 s _ 1 = m s - 1 g gravitational acceleration, 9.81 m s - 2 h sea surface height, m H salinity, %o i spatial discretization index in x or i? direction j spatial discretization index in y or A direction K global sea-ice matrix xi List of Symbols I ocean layer index in vertical direction LWL iceberg waterline length, m rh volumetric melt rate, m 3 d - 1 Mi iceberg mass, tonnes n general exponent Niayer number of ocean layers Naize number of iceberg sizes P cell ice melted fraction, % Qy total cell volumetric flux, m 3 s _ 1 r general position in Cartesian or spherical coordinates R Earth radius, 6371 km R rotation matrix s iceberg size index t time, d T temperature, °C U general continuous velocity in x direction where U = U(x,y,t), U' general discrete velocity in x direction where U' = U'(t), m s - 1 U velocity magnitude y/U2 + V2, m s - 1 XJAI velocity magnitude of air with respect to ice in x direction, m s~ U# meridional angular velocity, ° W s _ 1 UAI velocity of air with respect to ice in x direction, m s - 1 UG geostrophic velocity of water in x direction, m s - 1 Utide tidal current velocity magnitude, m s _ 1 xii List of Symbols Utide tidal current velocity in x direction, m s - 1 Uw water velocity magnitude, m s - 1 Uwi velocity magnitude of water with respect to ice in x direction, m s - 1 UWJ velocity of water with respect to ice in x direction, m s _ 1 V general continuous velocity in y direction where V = U(x,y,t), ms V' general discrete velocity in y direction where V' = V'(t), m s - 1 V\ zonal angular velocity, °N s _ 1 VAI velocity of air with respect to ice in y direction, m s - 1 VG geostrophic velocity of water in y direction, m s _ 1 Vtide tidal current velocity in y direction, m s _ 1 Vwi velocity of water with respect to ice in y direction, m s - 1 W solar radiation factor x general map coordinate, m X general position coordinate where X = X(t), m y general map coordinate, m Y general position coordinate where Y — Y(t), m 3 ice sediment fraction, % 7 angle measured between U and U , °CCW A latitude, °W fi inverse iceberg volume, m - 3 A4 total melt rate at waterline, m d _ 1 A4a air convective melt rate at waterline, m d _ 1 A4C water convective melt rate at waterline, m d _ 1 Xlll List of Symbols Mr solar heating melt rate at waterline, m d M.w wave action melt rate at waterline, md p inverse cell area, m~2 PA air density, tonnes m~ 3 pi ice density, tonnes m - 3 pw water density, tonnes m~3 $ v volumetric flux per unit area, m 3 m~2 s ip cell ice volume, m 3 •& longitude, °N v ice volume per unit area, m 3 m - 2 = m V total ice volume, m 3 xiv Chapter 1 G L A C I O M A R I N E M E C H A N I C S A N D H E I N R I C H E V E N T S 1.1 Introduction Heinrich events are among of the more perplexing features of the last glacial period. The marine sedimentary record shows that at least six and possibly seven of these massive discharges of icebergs from the Laurentide ice sheet, through Hudson Strait, occurred at roughly 10,000 year intervals. Heinrich events are a topic of increasing interest largely owing to their possible climatic implications. Substantial evidence indicates that the last Heinrich event, which took place roughly 11,000 years B.P., triggered the Younger Dryas cold event. Theory holds that meltwater flux from decaying icebergs drastically disrupted, and even shut down, the thermohaline circulation of the North Atlantic. Because of this perturbation, heat transfer from the tropics to the mid-latitudes and north is believed to have been greatly curtailed. The following briefly summarizes some of the issues associated with Heinrich events. This entails a general introduction to iceberg production and drift as well as reviews of Heinrich event sedimentary records. This preamble sets the stage for the development of a numerical model to investigate Heinrich events in Chapter 2. For the reader interested in the mathematical details, the appendices, which can be read stand-alone from the main body, contain the theoretical and numerical formulations of this model. 1 Chapter 1. Glaciomarine Mechanics and Heinrich Events 2 1.2 Iceberg Sources and Production Icebergs are fragmented from glaciers and ice shelves by processes collectively known as calving. The most important modern sources of icebergs are Greenland and Antarctica. Estimates of iceberg production are based on glacier speed and cross-sectional area. In the North Atlantic, the region of primary concern to this work, present-day Greenland fluxes have been variously calculated at 225 GT/yr (250 km 3/yr) [Paterson, 1994] and between 205 and 265 GT/yr (227 and 294 km3/yr) [Robe, 1980]. The total production is almost evenly distributed between the east and west coasts of the island. While there are calving glaciers in the archipelago of the Canadian Arctic, they are minor and irregular sources; in recent memory the only notable iceberg producer has been the Ward Hunt Ice Shelf of Ellesmere Island. Yet the warm stable global climate of the Holocene period, from about 11,000 years B.P. to present, has seen iceberg production in the North Atlantic that is a pale reflection of past activity. During the Pleistocene (approximately 1,000-11 kyr B.P.) major calving glaciers were found in Greenland, Iceland, Svalbard, Scandinavia and continental North America. Locations of past and present iceberg sources are shown in Figure 1.1. Corresponding plots of modelled iceberg production rates from palaeoclimatic regimes to present, of which more will be said in Chapter 2, are shown in Figure 2.1. The external characteristics of an iceberg: shape, size and surface features, are func-tions of the morphology of the calving parent ice mass and the fracturing mechanism active in the calving process. For example, tidewater glaciers without floating termini will produce small icebergs. Constrained to sliding motion and subject to the resultant internal stresses, the glaciers calve along internal weaknesses, where weakening is dic-tated by heterogeneities and thermal stresses. Of icebergs produced in this manner, the greatest dimension is typically much less than the glacier thickness. These icebergs are mainly of local significance, rarely travelling beyond the calving bay. Chapter 1. Glaciomarine Mechanics and Heinrich Events 3 Past and Present North Atlantic Iceberg Sources Figure 1.1: Iceberg past (in brackets) and present sources. Modeled calving rates are provided in Figure 2.1. The arrows indicate ocean currents and, correspondingly, drift patterns. Ice shelves and the floating tongues of ice streams, on the other hand, behave as elastic cantilever beams acted on by flexural loads in addition to internal forces [Robe, 1980]. Depending on the aspect ratio or slenderness of a floating tongue and the conditions of its grounding line, failure may occur by lack of buoyancy causing the ice beam to bend under its own weight. Conversely, for ice shelves as well as tidewater glaciers and ice streams, fluctuations in water level by waves, tides and storm swell can lead to excessive buoyancy. Repeated wave action may also excite natural frequencies in the floating portion of the ice mass leading to fatigue failure. In these flexural modes calving is generally expected to occur at the support location or grounding line, though weak spots can cause failure Chapter 1. Glaciomarine Mechanics and Heinrich Events 4 I C E S H E L F G R O U N D I N G L I N E Figure 1.2: Iceberg calving. closer to the terminus. B y these buoyancy mechanisms, relatively long and shallow icebergs are made, so-called tabular icebergs. Icebergs calved from ice shelves can reach length scales exceeding 100 km. Calving of terminal glaciers and ice shelves is schematically illustrated in Fig-ure 1.2. 1.3 Ocean Drift Patterns In contrast to the dearth of knowledge about iceberg size distribution, there is an abun-dance of information about drift trajectories. Since the sinking of the Titanic in 1912, considerable interest has been focussed on understanding iceberg drift so as to avoid collisions in shipping lanes. Those efforts have been renewed in recent years with the advent of offshore gas and oil development in the North Atlantic . Chapter 1. Glaciomarine Mechanics and Heinrich Events 5 A brief synopsis of the drift of Arctic Ocean icebergs gives an indication of the mean-dering nature of iceberg drift. In the polar basin of the Arctic, iceberg sightings are rare notwithstanding the lack of observers. Of those that have been seen, electronic tracking of Ellesmere Island icebergs, for instance, has found them to drift in the Arctic seas at an average speed of about 2 km/day, returning near their origin after a 10-12 year tour [Weeks and Campbell, 1973]. Though high-Arctic icebergs rarely stray from their polar enclave, a few are caught in the East Greenland Current and carried to the North Atlantic. Flowing southward between Svalbard and Greenland, the East Greenland Cur-rent is the major modern outlet for surface water from the Arctic Ocean; iceberg speeds have there been measured in the range of 7-22 km/day [Robe, 1980], reflecting the speed of the water surface. As an example, between 1960 and 1965 icebergs with waterline dimensions up to 600 m, believed to be of Ward Hunt Ice Shelf provenance, were sighted off the coast of Labrador and later observed as far south as 44° N [Robe, 1980]. Pack ice, the congregation of ice floes that form on the ocean surface, scarcely retreats above 70° N, which also delineates the southern limit of most Greenland iceberg produc-tion. Often held fast to the shore, pack ice can impede the drift of icebergs, confining icebergs to the sounds or fjords from which they calved. But, when not held to the shore, the entrapping pack ice itself drifts. Most Arctic basin icebergs that reach the open ocean drift south with pack ice on the East Greenland Current. Not only does the pack ice shelter icebergs from the effects of waves, but the cold surface waters largely prevent melting. Landfast ice may also prevent icebergs from running aground, thus promoting their longevity. Drifting with the East Greenland Current, some icebergs catch the warmer Irminger current, relatively temperate at 4° to 6° C, and travel along the northern coast of Iceland, east of which waves and high temperatures quickly cause decay. Rapid deterioration, by wind and waves, characterizes the general state of any iceberg below 70° N. Chapter 1. Glaciomarine Mechanics and Heinrich Events 6 An indeterminate, but small, number of icebergs from Eastern Greenland ultimately joins the West Greenland Current. Drifting north to Baffin Bay the icebergs successively join the Baffin Island Current and then the Labrador Current south of Hudson Strait. Southerly-bound by the Labrador Current, icebergs eventually reach the Grand Banks off Newfoundland. Diverging courses take the icebergs east to drift north of Flemish Cap or south between Grand Banks and Flemish Cap. The southern fringe of iceberg drift is generally the northern edge, within about 12°, of the North Atlantic current. Within the last century, icebergs have been occasionally spotted as far east as the Azores and as far south as Bermuda [Dyson, 1972]. According to the International Ice Patrol of the United States Coast Guard, the number of icebergs drifting south of 48° N in any year shows no discernable pattern [Robe, 1980]. 1.4 Drift Forces The forces acting on an iceberg, as with other floating bodies, arise from the relative motions with respect to air and water as well as inertia. In determining the drift traject-ory of an iceberg the principal forces are due to air and water drag, Coriolis acceleration and gravitational pull associated with sea-slope. In the hierarchy of forces, water drag is primary and air drag is typically an order of magnitude less. Coriolis and sea-slope forces generally cancel and are, in any case, of a lesser order than air drag. (See Appendix A, Section A.2.1 for a specific discussion of iceberg drifting dynamics.) A major limitation to conducting studies has been that field calculations of the forces acting on icebergs are difficult and expensive to collect. This is in part owing to uncer-tainties such as iceberg size, shape and surface roughness and ocean stratification. For example, drift predictions have most often been based on surface currents, which may flow counter to the drift of an iceberg depending on the velocities in the sub-surface layers Chapter 1. Glaciomarine Mechanics and Heinrich Events 7 in which the iceberg is immersed. 1.5 Deterioration Icebergs receive radiative energy from the sun, as well as heat conducted and convected by the surrounding air and water. Although all these fluxes contribute to melting, it has been suggested that direct and indirect radiation cannot account for even one percent of the energy needed to melt an iceberg [Robe, 1980]. As with accurate force measurements, field measurements of melting icebergs are virtually unobtainable. Icebergs are not only exceedingly unsafe to work near, but the melting process is remarkably complicated. As icebergs melt below their waterline, buoyant fresh meltwater mixes and rises in the surrounding saltwater in a process where turbulent heat transfer dominates [Huppert, 1980; Robe, 1980]. In short, waves remove the cold freshwater at the ice surface and deliver warm sea water in replacement. Furthermore, the freeboard portion of an iceberg exposed to air is also subject to turbulent heat exchange. Given typical iceberg dimensions, surface irregularities and wind speeds, the air flow over the surface of an iceberg is characterized by Reynolds numbers on the order of 106 to 107. The best information available on iceberg melting is highly anecdotal based on field observations or empirical based on in situ laboratory experiments [Huppert, 1980; Russell-Head, 1980]. For instance, the volume of an iceberg must be known in order to predict melt rates. Various rule-of-thumb methods have been proposed for predicting iceberg volume; volume estimates range from 3.35lwh to 4z.9lwh [Robe, 1980], where I, w and h are length, width and height, respectively. These seemingly haphazard approximations inescapably surface where accurate data are scarcely obtainable. Having thus qualified melt predictions, the reader is referred to Appendix D where an empirical scheme to predict iceberg melt flux is presented. Chapter 1. Glaciomarine Mechanics and Heinrich Events 8 As remnants of their glacial past, icebergs are strewn with cracks and crevasses. Breakup thus plays a significant role in deterioration. Sheltered from the open sea, icebergs are stable and stress free relative to their state in the open ocean. Fully exposed, waves cleave at planes of weakness: indentations, embayments and fractures. The form of an iceberg sometimes tells of the wave action that has taken place. Waves can cut through the top of an iceberg leaving the bottom intact; the resultant is distinctively marked by spires of ice above the water, separated and attached by an underwater portion, a so-called drydock iceberg. Icebergs that have experienced continuous rolling may have a smooth, rounded topography where waves gradually washed over the surface. Where the depth of an iceberg exceeds its width, rolling is unstable and will likely lead to over-turning. For tabular icebergs, with freeboards of perhaps 40 m, undercutting by wave action leads to promontories that fail in shear or flexure, depending on aspect ratio. Moreover, icebergs with waterline dimensions greater than the wavelength of the water waves may have susceptible resonant frequencies leading to flexural loading and fatigue. In the foreseeable future it appears that efforts to research iceberg drift and deteri-oration will continue, driven in large measure by commercial interests, not to mention for purely scientific purposes. As long as icebergs pose a hazard to shipping and ma-rine structures, quantifiable assessments of risk will need to be made based on the best available information. 1.6 A Curious Sediment Record The following introduces the principal features of the sedimentary record that first in-dicated the existence of Heinrich events and presents an overview of the ocean physics associated with them. The intention is to acquaint the reader with the relevant back-ground issues that motivate the model development in Chapter 2. Chapter 1. Glaciomarine Mechanics and Heinrich Events 9 Over a wide swathe of the mid-latitude North Atlantic (40° N to 55°N) sediment records show six distinctive horizons of large scale marine deposition. First reported in Heinrich [1988], these sediment layers are attributed to rapid and voluminous ice-rafting. Lead isotope analyses later found the anomalous deposits of ice-rafted detritus (IRD) to be primarily of Hudson Bay provenance [Gwiazda et al, 1996], lending credence to the notion of massive episodic calving of sediment laden icebergs from the Laurentide Ice Sheet (LIS) during what have come to be known as Heinrich events [Broecker et al, 1992]. Various estimates have these sediment layers being deposited over the course of a few hundred to around 2,000 years. Heinrich layers stand in marked contrast to adjacent layers from other, generically termed, detrital events which do not show such a significant Laurentide component [Elliot et al., 1998]. Rather, in those layers Greenlandic, Icelandic and Nordic ice-rafted sediment signals are prevalent. Drill core studies have focussed on the mid-latitudes where Heinrich IRD is clearly distinguishable from local detrital events near their sources. Attenuation of the Heinrich sediment signal in the distal lower latitudes imposes a southern limit on investigations. Heinrich sediment layers are starkly different from the background IRD cover in terms of grain size, age and mineral composition. Though perhaps the most salient feature is a notable lack of fossilized planktonic foraminifera N. Pachyderma [Bond et al, 1992]. The foraminiferal concentrations in the Heinrich layers are in fact of an order or two less than usual Pleistocene values and represent a greater depletion than simple dilution in massive sediment allows. A cold ocean would have low foraminifera production. Reduced d>180 content in the fossils also indicates low salinity in the near surface water of the Labrador Sea [Elliot et al, 1998]. The local 6180 minima lag behind the initial Heinrich sediment signal, suggesting that the early depositions were suspended sediments churned by the advancing ice sheet. The gradual build-up to the peak IRD signal marks the progression Chapter 1. Glaciomarine Mechanics and Heinrich Events 10 and maximum extent of iceberg melting [Clarke et ai, 1999]. The evidence implies that Heinrich events occurred at the height of cooling cycles. Greenland cores and pollen records tell that during the Pleistocene two orders of cooling cycles operated. So-called Bond cycles enveloped a bundle of abrupt temperature oscillations termed Dansgaard-Oeschger (D-O) cycles, with periods of around 10,000 years and 500 to 2,000 years, respectively. While the exact relationship between D-0 and Bond cycles remains unclear, both seem to have roots in the North Atlantic. Heinrich events appear to be in phase with Bond cycles, occurring at their culmination during the cooling phase of the last and coldest D-0 cycle in a series. These are followed by warming and a series of successively colder D-0 cycles until the peak of the Bond cycle and then another Heinrich event. Figure 1.3 shows a Greenland <5180 record with the timing of Heinrich events noted. The <5180 record indicates that Heinrich events started while the ocean was relatively cold. Local SlsO minima correspond to surface salinity minima or meltwater peaks. The temporal relation between Bond and D-0 cycles and Heinrich events is schematically illustrated in Figure 1.4. The repeated D-0 fluctuations are attributed to a series of abrupt changes in ocean circulation perturbing poleward heat transport [Broecker, 1994]. Before the onset of the Holocene, the final Heinrich event at the end of the Pleistocene ushered in a brief final glacial period. Known as the Younger Dryas, this cold episode was most noticeably felt in the North Atlantic and adjacent land masses. Coincidental timing leads to speculation that a Heinrich event had some role in disrupting the ocean circulation during the Younger Dryas. The unanswered question of how the Younger Dryas event and last Heinrich event are related calls for further investigation. Chapter 1. Glaciomarine Mechanics and Heinrich Events 11 Summit Greenland GRIP oxygen-isotope record 0 (0 20 30 « 50 fO 70 10 T i m e ( k y r B P ) Figure 1.3: Greenland S180 record indicating Heinrich events [Dansgaard et ai, 1993]. The S180 concentration is plotted as the deviation from the standard mean ocean value. S180 minima indicate peak cooling. D a n s g a a r d -O e s c h g e r C y c l e Bond Cycle T i m e ^ ~ 7 to 10 kyr Figure 1.4: Schematic relation between Bond, Dansgaard-Oeschger cycles and Heinrich events. Chapter 1. Glaciomarine Mechanics and Heinrich Events 12 1.7 Thermohaline Circulation The ocean circulation comprises wind- and thermohaline-driven components. Surface-imposed shear stresses make wind-driven circulation dominant in the upper ocean layers, while thermohaline flows prevail at depth. Energy from the Sun unevenly heats the ocean and atmosphere and is ultimately responsible for both types of ocean circulation. Thermohaline circulation is driven by changes in water density as governed by tem-perature and salinity fluctuations with climate. The cool winter of high latitudes causes water to lose heat and increases density. Sinkage of the dense surface displaces bottom water which, in turn, migrates toward the lower density regions near the equator, thus es-tablishing a conveyor system. One segment of the global thermohaline circulation system active in the Atlantic is known as the North Atlantic Conveyor; it delivers heat poleward from the warmer latitudes to moderate the climate of mid- to high-latitude Europe and North America. The effect of surface cooling is enhanced by the formation of pack-ice, which increases surface salinity and maintains low temperatures. The strength of the thermohaline circu-lation is not a function of absolute density, rather, it depends on gravity potential or the density gradient. This gradient is not only controlled by events at the ocean surface, but also by density changes in the lower layers. Bottom water buoyancy is increased through the eddy diffusion of heat from solar radiation, removal of salt and, to a minor degree, ge-othermal heat flux. Continuous decay and upwelling of the bottom water accommodates thermohaline circulation; it would otherwise cease. Essentially, thermohaline circulation originates as a vertical gravity flow in the high-latitudes leading to a meridonial flow at depth, satisfying volume continuity. Of the direct effect of differential heating between the poles and equator, its contribution to ocean circulation is thought to be minimal [Pickard, 1964]. In warm climates where the contrast between surface and deep water Chapter 1. Glaciomarine Mechanics and Heinrich Events 13 Figure 1.5: Schematic thermohaline circulation. temperatures is small, thermohaline circulation is sluggish. In the cool waters of high latitudes where warm poleward currents strongly influence climate, thermohaline circu-lation is relatively vigorous. Confirmed by numerical models [e.g. Manabe and Stouffer, 1993; Weaver et al, 1996], the prevailing wisdom holds that the strength of the thermo-haline circulation is highly sensitive to fresh water input near the polar regions. Because the vertical deep-water formation occurs over a small areal extent at high latitudes, it is believed that the thermohaline circulation is acutely vulnerable to interruption. The findings of Heinrich [1988] have lead to the suggestion of a link between massive glacial calving and disruption of the North Atlantic Conveyor. The proposition of vast iceberg armadas is supported by retrieval of nearly identical cores from other drilling programs [Broecker et al., 1992]. Iceberg melting is a plausible mechanism for switching Chapter 1. Glaciomarine Mechanics and Heinrich Events 14 the North Atlantic thermohaline circulation between interglacial and glacial modes. This is consistent with evidence from Greenland ice cores and North Atlantic marine cores that abrupt climate changes are associated with changes in ocean circulation [Broecker, 1994; Fanning and Weaver, 1996]. 1.8 Global Links Mounting evidence indicates that Heinrich events were strongly linked to a variety of other global scale climatic phenomena and that they occurred at major climatic boundaries [Broecker, 1994]. Greenland cores reveal troughs and peaks in methane levels that follow the coldest D-0 cycles and Heinrich events [Chappellaz, 1993]. Presumably, the tropical wetlands were the prime source of atmospheric methane. Pollen records from Florida hint at a correlation between the last Heinrich event and a high precipitation period in what is now the southern United States, where a switch in the North Atlantic Conveyor would not be directly consequential [Grimm, 1992]. It is suggested that this approximately 1,000 year wet period had features of an El Nino event; though the link with the North Atlantic thermohaline circulation is unresolved. It is conjectured that a change in the North Atlantic Conveyor could be manifest in the tropics by an alteration in moisture convection to the atmosphere. Further afield, evidence indicates synchronicity with ice advances in the Andes and Fennoscandia, among other examples [Broecker, 1994]. Starting as mere oddities in the sedimentary record, the evidence suggests that Hein-rich events had major climatic implications. The investigation of Heinrich events is an exercise in unravelling the connections between the oceans, atmosphere, cryosphere and global climate where the divisions separating the causes from the responses are murky. Chapter 1. Glaciomarine Mechanics and Heinrich Events 15 The present interest in Heinrich events underscores timely issues, for instance, those re-lating to modern ocean circulation, El Nino events and the warming of polar regions. Understanding the role of the meltwater fluxes associated with Heinrich events might shed light on the possible implications of the melting of the present polar ice caps, as an example. Having established the motivation for further investigation, the next chapter presents a model to examine Heinrich event iceberg drift, meltwater flux and sedimentation. Chapter 2 A CONTINUUM MODEL OF ICEBERG DRIFT: APPLICATION TO HEINRICH EVENTS 2.1 Iceberg Drift Models: State of the Art to Present Until the latter half of the nineteenth century, notions of iceberg drift were generally ill-conceived, simply because icebergs scarcely interfered with human activity. Growing cargo and passenger traffic along ice-infested trans-Atlantic shipping routes at the turn of the century meant that iceberg avoidance became an increasing problem. Driven by the sinking of the Titanic in 1912, the newly formed International Ice Patrol (IIP) of the Customs Cutter Service, forebear of U.S. Coast Guard, introduced a formal approach to predict short-term iceberg drift. While knowledge of the original method has been lost through the years, iceberg trajectories were essentially hand calculated by vector addition of shipboard measured wind and water current velocities. Though the IIP vectoral method evolved somewhat, manual calculations prevailed until the advent of computational methods in 1971. Although not the first automated technique adopted by the IIP, the model of Moun-tain [1980] lias been successfully used for trajectory predictions to present [Armacost, 1995]. This model comprehensively includes wind drag, water drag, Coriolis and gravi-tational (sea-slope) forces to formulate a set of differential equations of iceberg motion. It is used to predict the short-term, mainly 12 hour, drift of certain icebergs threatening shipping lanes. Real-time data, namely air and multi-layered water current velocities 16 Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 17 collected by buoys and drogues, are used in conjunction with iceberg mass, estimated based on above-water dimensions, to calculate forces. Similar models were independently developed, such as those of Napoleoni [1979] and Sodhi [1980]. Compared to drift dynamics, iceberg melting physics is poorly understood, in large part because icebergs are mostly inaccessible and the bulk of their mass is unobservable below the waterline. The flow regimes of heat-bearing air and water about the irregular surface of an iceberg is characteristically highly turbulent, quantified by Reynolds num-bers on the order of 106 to 107, as previously noted. Theoretical arguments to formulate a realistic description of iceberg melting do not exist, so the recourse for researchers has been to depend heavily on empirical derivations from laboratory experiments. The expected result has been that the various drift and deterioration models agree on the drifting physics, but diverge notably in their treatment of melting. Drift models to date have been very limited in scope; the computational expense of solving differential equations of motion and the sparcity of input data have confined model runs to single icebergs over short distances and times. For the purposes of avoiding collisions with ships and offshore structures, as dictated by the mandate of the IIP, for example, short-term trajectory predictions have sufficed. This approach, though, is inadequate for the investigation of large scale drifting, as for Heinrich events, and associated climatic phenomena. Of recent modelling endeavours with respect to iceberg drift, the model of Mat-sumoto [1996] is notable. The model predicts drift trajectories for multiple icebergs by the method of Mountain [1980], with the additional routines for predicting ice rafted depostion (IRD) and meltwater fluxes. Melting is governed by an empirical relation de-scribing iceberg life expectancy as a function of mass and geographical location. IRD release, in turn, varies with iceberg mass. Seasonality is included in the model by the use of observed monthly average air velocity fields. As well, the life expectancy equation Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 18 is a function of empirically determined terms that account for month. Water velocity fields, that with air velocities determine drifting dynamics, come from general circulation models (GCMs) for present-day conditions. While the work of Matsumoto [1996] seems to be the most rigorous among other contemporaneous drift/sedimentation models [e.g. Dowdeswell and Murray, 1990] it is, nonetheless, very limited in two obvious aspects. First, the numerics of Mountain [1980] are not well suited to modelling large numbers of icebergs. The method tracks individual trajectories, which for multiple simultaneous icebergs becomes cumbersome and compu-tationally expensive, especially if the modeler wants to predict iceberg drift over an area as extensive as the entire North Atlantic. This immediately curtails the scope of inves-tigations that can be reasonably carried out with the model of Matsumoto [1996] and, therefore, its ability to predict IRD and meltwater flux to their full extent. Secondly, even given that melting is a poorly grasped subject, iceberg decay in Matsumoto [1996] is crudely dealt with. In the model, melting is basically the product of iceberg mass and a month factor. Clearly, this is a rigid approach as the melting equation is only valid for the ocean region for which the empirical relation, based on observations, was formulated. Variables that are known to govern melting, among them temperature and salinity, are not explicity included. Actually incorporating the physics of heat transfer, albeit also empirically determined, would at least give the model generic appeal such that it could be applied to any ocean region. Thanks to GCMs, global ocean synthetic data are widely available. Along this vein, the Matsumoto [1996] model uses only present data, which are applied to model certain paleocean conditions which were believed to be similar to the present ocean. For model runs to investigate glacial period iceberg drift, which happens to be a topic of burgeoning interest, the Matsumoto [1996] model could not possibly apply. The object of the preceding criticism was not to deconstruct the admirable work Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 19 of Matsumoto [1996], but rather to demonstrate that the modelling of iceberg drift, sedimentation and freshwater flux could yet be improved. This is especially true where it is desired to examine ocean or global scale events. 2.2 Physics of Drifting As noted, iceberg drift is governed by wind drag, water drag, Coriolis and sea-slope forces. In a multi-layered ocean the discrete differential equations (ODEs) of motion of iceberg drift are written as dX dt = u (2.1) and dt (2-2) for velocities and as dU' dt fV - fVG + ^layer PACDAAA\CAI\UAI + PWCDW ^ AWL\CWil\Uwil 2MT (2.3) and dV' dt -fU +fUG+ N,a PACDAAA\UAI\VAI + PWCDW ^ ^ W , | U W J , | V V I , i=i 2MT (2-4) for accelerations in the x and y directions in Cartesian coordinates. The notation is listed at the beginning of the report on page xi. In a state of geostrophic equilibrium, the sea-slope pull and Coriolis accelerations ex-actly balance, thus sea-slope is expressed as the negative Coriolis acceleration of the local Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 20 geostrophic current, which gives the — JVG and fU~G terms in (2.3) and (2.4), respectively. This is further elaborated in Appendix A, section A.2.1. Mountain [1980], as typifies other models, solves (2.1) through (2.2) via a fourth-order Runge-Kutta (RK4) scheme. Melting may likewise be solved as an ODE describing the time rate of change in a characteristic dimension or mass. The IIP approach, for example, calculates the rate of change in iceberg length as a function F of, among other factors, temperature and present length. Generically, this is expressed as ^LWL = F(T,LWL,...). (2.5) (See Appendix D for details of the IIP melting routine.) For numerical stability, time stepping in the RK4 solution of (2.1) through (2.4) typ-ically ranges from seconds to a few minutes, in inverse proportion to iceberg velocity and acceleration. For the purposes of predicting gross meltwater and sedimentation patterns over, perhaps, a millenial time scale, modelling the discrete trajectories of icebergs is cumbersome. 2.3 A New Approach: Icebergs as a Continuum Iceberg drift models require for their operation input ocean and atmosphere data fields. For gross sediment and meltwater predictions, GCMs are an ideal source of these input data; real data are difficult and expensive to collect and, resultingly, relatively sparse. Moreover, GCMs can generate data for past ocean conditions, making feasible historical investigations. Given that the ability to couple with a GCM is a desirable feature of an iceberg drift model, the solution that naturally presents itself is to use the same modelling approach of GCMs, namely continuum modelling. Imagine an armada of icebergs with each member having the same mass. Assume that these icebergs are immersed in a uniform flow field, so that each individual experiences Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 21 the same air and water currents. Starting from rest, over time the icebergs will spread apart, although they will tend to drift in the same direction. In the uniform flow field all these equimass icebergs will eventually reach the same constant velocity. Until that time, however, they will variously accelerate according to their outward characteristics, namely shape and surface roughness, which effect the way the fluid flow acts on them. Inspection of (2.3) and (2.4) reveals that acceleration is linearly dependent on shape and surface parameters, in the form of drag coefficients. As the icebergs accelerate they will appear to diffuse, not in the classical sense of seeking a lower potential, but rather by the variable forces acting on them due to their individual properties. For a range of drag coefficients, suggested to be uniformly distributed between 0.6 and 2.0 [Russell-Head, 1978], the trajectory of a hypothetical iceberg with the midpoint drag coefficient represents the mean of all the iceberg trajectories within the armada. From this basis, the trajectories of the icebergs in the armada could be said to have two parts: (1) a determistic component, the trajectory of the hypothetical iceberg and (2) a probabilistic component, the stray from the mean trajectory, a function of the distribution of drag coefficients. Depending on their individual Coriolis accelerations the armada icebergs will deflect from straight line paths that the drag forces alone would determine. In a real ocean where flows are not uniform, icebergs in close proximity are variably subject to fluctuating air and water currents. The net result is that the icebergs will not only be scattered along the flow direction, but also in the cross-stream direction. In a mathematical description, the armada could be formulated as an advective-diffusive continuum of equal mass as the sum of all the individual iceberg masses. The advective component would be calculated as the constant velocity that the icebergs will achieve. Using the constant velocity may seem inaccurate; it is an approximation ignoring the acceleration period. But, in fact, icebergs generally accelerate over relatively short time frames. Over a 2° by 2° resolution, characteristic of which GCMs operate, whether Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 22 an iceberg is considered to start at rest or with the constant velocity has negligible influence on the calculated trajectory after as little as a couple days. (See Appendix A, section A.2.1 for more explanation.) The diffusive component, on the other hand, would comprise diffusion coefficients de-scribing the rate of drift along- and across-stream from the mean trajectory. Stochastic diffusion coefficients could be experimentally determined through the analyses of com-puted trajectories found by the solution of the equations of motion, (2.1) to (2.4). If icebergs generally follow ocean currents, then the local dispersion is probably minimal. Nonetheless, diffusion necessarily has to be included to produce a realistic picture of the spread of ice-rafted sediment. Also, while the local wayward drift may be small, as the icebergs experience new flow regimes they will increasingly spread. In spherical coordi-nates with the addition of melting and source terms, the advective-diffusive formulation of iceberg drift is written as Please refer to Appendices A and B for full derivations and descriptions of the continuous model and stochastic diffusion. 2.4 A New Model The following presents a model of iceberg drift and sedimentation based on the con-tinuum formulation established in the previous section. The model solves (2.6) over a discretized problem domain using a five-point finite-difference scheme that is implicit in time with upwind differencing for the advection terms. Time stepping is determined by the Courant-Friedrichs-Lewy (CFL) condition for numerical stability based on local iceberg drift velocity and cell spacing. Systems of equations are solved by the method of (2.6) R? sin A [d& Vsin \ d-d) d\\ 1 \6_ ( Ds dv\ d_ ( Dx sin A — ) = §v + rhpv. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 23 successive or simultaneous over-relaxation (SOR) described in Appendix E. The discret-ization and grid spacing, which may be non-uniform, are taken as user input. Details are provided in the appendices. The model requires input iceberg sources and ocean fields, for which a continental ice sheet model model (ISM) [Marshall, 1996] and a coupled ocean-atmosphere climate model (CCM) [Weaver et al., 1998] are used to generate data. It is informative to introduce the synthetic ice sheets and ocean upon which the iceberg drift model depends. 2.4.1 Iceberg Calving Sources Iceberg source data are gathered from a three-dimensional, thermomechanical ISM. The thorough treatment given to ice sheets includes viscous creep, bed decoupling and sedi-ment deformation and, importantly from the perspective of calving, fast flow instabilities. For the North Atlantic, with which our investigation of Heinrich events is concerned, the ISM supplies calving rates in terms of the annual volumetric fluxes from 10 sources of four provenances: North America (NA), Greenland (GL), Fennoscandia (FS) and Iceland (IS). Plots of these flux rates sampled every 100 years from 120,000 years B.P. to present are shown in Figure 2.1. In the iceberg drift model, each ice source cell on the discretized problem domain is ascribed a provenance label according to its modelled geographical location; indices 1 to 4 imply NA, GL, FS and IS, respectively. Each provenance in the model has a unique probability density function (PDF) associated with it describing the distribution of ice-berg sizes at its calving sources. This allows that the Greenland ice sheet would tend to produce much larger icebergs than would be manufactured at Icelandic sources, for example. The area under each PDF is divided according to size category; size produc-tion is then calculated as the product of the PDF size category area and the total source production rate. The number of different sized icebergs in the model is taken as user Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 24 /a£>rcrcfc>/-North American sources — SoOnSor 10 Icelandic sources S K. l u J j •JU. JN 120 100 80 60 40 Time (kyr B.P.) 20 Figure 2.1: Input model ice sheet calving sources [S. Marshall, unpublished]. The data are plotted every 100 years. Source locations are indicated on the map of Figure 1.1. input. Waterline length, or some other characteristic dimension, is considered the meas-ure of size. Assuming some standard shaped iceberg in the model, which is mitigated by applying diffusion, iceberg volume is then a function of the characteristic dimension. (Refer to Appendix C for specific details.) Heinrich fluxes for Hudson Strait are not included in the time series of Figure 2.1. For Heinrich event runs, the Hudson Strait flux is augmented by 200 km 3/yr, a rough estimate of the mid-range Heinrich flux value from Marshall and Clarke [1997]. Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 25 2.4.2 Ocean and Atmosphere Fields The iceberg drift model uses a variety of atmosphere and ocean data: atmosphere and ocean temperature and velocity fields; ocean salinity fields and an ocean bathymetry field. The data are output by the CCM for 19 ocean layers and a single atmosphere layer. The largest icebergs in the drift model, and in nature, will only sink to a depth of around 500 m or the eighth model ocean layer. The CCM atmosphere is rudimentary, simply consisting of energy and moisture balance equations for heat and fresh water. The ocean component comprises a three-dimensional, 19 layer, ocean GCM (OGCM). The OGCM was originally developed as the Princeton Geophysical Fluid Dynamics Laboratory -Modular Ocean Model, which is detailed in Weaver and Hughes [1996]. The C C M data, as supplied for the drift model, was obtained on a model grid res-olution of 3.75° zonally and 1.8555° meridionally. For use in the iceberg drift model runs discussed in the next section, these data were linearly interpolated to give a 1.875° zonal resolution. CCM runs produced data for the Last Glacial Maximum (LGM), around 21,000 years B.P., and for modern, around 10,000 years B.P. to present, ocean conditions. The data are in terms of their mean annual values. LGM model conditions are applied to runs throughout the last glacial period. In this respect, the drift model produces first-order iceberg drift predictions. An ideal model would have continual feedback between melting icebergs and the ocean. For runs that cross the glacial to present ocean date threshold, the model automatically switches ocean conditions. Figures 2.2, 2.3 and 2.4 illustrate the CCM input fields used by the iceberg drift model for the runs to be discussed. The data are plotted on a 60x30 grid using the iceberg drift model resolution already noted. Ocean data are shown for the surface layer only Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 26 Ocean Currents: LGM Ocean Salinity: LGM 20°W Ocean Salinity: Present Ocean Currents: Present -10 - 8 - 6 - 4 - 2 Salinity (ppt - 35) Figure 2.2: Input model ocean fields [Weaver et al, 1998]. Ocean surface layer velocity and salinity fields are illustrated on a 60x30 grid with 1.875° zonal and 1.8555° merid-ional resolution. Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 27 - 2 0 -10 0 10 Temperature (Degrees Celsius) Figure 2.3: Input model ocean fields [Weaver et al., 1998]. Ocean surface layer and atmosphere temperature fields are illustrated on a 60x30 grid with 1.875° zonal and 1.8555° meridional resolution. Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 28 O c e a n S t r a t i f i c a t i o n : L G M t o P r e s e n t 0 2 4 6 8 10 12 14 16 18 Number of ocean layers Figure 2.4: Input model ocean fields [Weaver et ai, 1998]. Ocean bathymetry field illus-trated on a 60x30 grid with 1.875° zonal and 1.8555° meridional resolution. 2.4.3 Drift Dynamics Air and water velocity fields determine iceberg drift velocity through the solution of the equations of motion. For each iceberg size at each cell in the model the velocity is calculated using (2.1) through (2.4) in an RK4 scheme (see Appendix A, section A.2.1). 2.4.4 Melting Melting in the drift model, using a variant of the the IIP method explained in Appendix D, is driven by temperature and salinity. The effective values of these parameters that an iceberg experiences are taken as the weighted vertical averages of the atmosphere and ocean layers in which the iceberg sits. These average parameters are calculated at each problem cell for each iceberg size. In the IIP melting scheme air velocity also plays a minor role in melting. For variables that govern melting that are not CCM output data: solar radiation and wave height and frequency, assumed conservative values are uniformly Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 29 applied over the model ocean. 2.4.5 Size Coupling The model accounts for the fact that as icebergs melt the fields that act on them change. The drift model keeps a running tally of the mass history of the ice as it moves from cell to cell. As melting is active, the portion of ice that is melted is gradually moved to the next lower size category, experiencing new drifting dynamics. Details of this procedure are given in Appendix D. 2.4.6 Sedimentation Sediment content of the model icebergs is a function of iceberg size and provenance. Within each size category, the sediment is assumed to be uniformly distributed through-out the icebergs. Through the coupling routine, as icebergs move to lower sizes their sediment content also changes. This attempts to mimic the sediment distribution of real icebergs in which sediment content decreases towards the centres. It seems crude to distribute sediment uniformly over surfaces. But, considering that icebergs may roll repeatedly as they drift, it is plausible that the effective sediment distribution does seem to be somewhat uniform viewed from the ocean bottom. For more information, please refer to Appendix D. 2.4.7 Sea-Ice Gating The iceberg drift model makes provisions for the presence of sea-ice, which can confine icebergs or impede their drift. This is accomplished in the diffusion routine by merely slowing, or possibly shutting down, the diffusion flux at a cell where sea-ice is present. Unfortunately, for the model runs discussed in the next section sea-ice data were not Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 30 available and therefore excluded. Sea-ice gating and other salient model features are graphically illustrated in Figure 2.5. 2.5 Model Runs Three 1,000 year iceberg drift model runs were conducted: (1) for L G M background conditions; (2) for L G M conditions with Heinrich event ice flux from Hudson Strait and (3) for modern ocean conditions. Lists of input parameters for these runs are provided in Tables 2.1 and 2.2. Parameter Units Default L G M L G M w/ Heinrich Modern start yr B.P. 21,000 21,000 2,000 duration yr 1,000 day 10 10 5 nx 60 ny 30 Ai? ° W 1.875 A A °N 1.8555 u mtn ° W -75 •d "max ° W 37.5 Amzn °N 28.76 Amaa; °N 84.425 1 * size 5 DA day- 1 0.005 DT day" 1 0.004 PI T m "3 0.900 PA T m " 3 0.001 Pw T m " 3 1.026 R km 6731 Table 2.1: Model run input parameters relating to the numerics and iceberg dynamics. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 31 Figure 2.5: Iceberg drift and sedimentation model features. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 32 Parameter Units Default NA GL FS IS LwLjmax m 500 LwL,avg m 90 90 75 75 Dist n(LWL) Rayleigh fimax % 1.0 Distn{(3) Linear Table 2.2: Model run input parameters relating to iceberg size and sediment composition by provenance. 2.5.1 L G M 1,000 Year Run Plotted in Figures 2.6 and 2.7 are the results for a model run under L G M conditions starting 21,000 years B.P. with a 1,000 year duration. Figure 2.6 shows the unmelted iceberg distribution, expressed as the equivalent ice column depth over the model cell. Barring an abrupt change in ice flux, steady-state ice distribution is generally reached after around 25 years. Figure 2.6 shows meltwater distribution, indicating the extent of iceberg drift. The meltwater distribution is also given as the equivalent ice column thickness. It is clear from Figure 2.6 that modern drift expectations do not hold for the modelled L G M results. Referring to Figures 2.2 and 2.3, the CCM predicts a significant weakening of the thermohaline circulation during the LGM through mechanisms including orbital forcing and continental freshwater run-off [Fanning and Weaver, 1996; Weaver et al., 1998]. The weakened thermohaline circulation with relatively cool, low salinity surface waters extending to the mid-latitudes would have promoted iceberg longevity. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 33 Equivalent Ice Thickness: 1,000 yr. LGM 0 1 2 3 4 8 6 7 Equivalent ice column depth (m) Figure 2.6: Unmelted ice distribution for LGM conditions. The results shown are for a 1,000 year model run beginning 21,000 yr B.P. Meltwater/Sediment Distribution: 1,000 yr, LGM 0 100 200 300 400 500 600 700 Equivalent ice column depth (m) Figure 2.7: Melted ice distribution for LGM conditions. The results shown are for a 1,000 year model run beginning 21,000 yr B.P. Chapter 2. A Continuum Model of Iceberg Drift: AppHcation to Heinrich Events 34 2.5.2 L G M / H e i n r i c h Event 1,000 Year R u n Figures 2.8 and 2.9 show the plotted results for a model run again using L G M conditions, but with a 200 km 3/yr Heinrich flux added at Hudson Strait. This is supposed to simulate Heinrich event H2. These results were also obtained for a run starting 21,000 years B.P. with a 1,000 year duration. The unmelted ice distribution Figure 2.8 is very similar to that of Figure 2.6 and also reaches steady-state after around 25 years. Figure 2.8 also does not show a pronounced difference from the meltwater distribution pattern of Figure 2.6, despite the difference in meltwater quantity. Noting that a 200 km 3/yr addition to Hudson Strait only represents an increase of about 30% in the North American sources, that the difference between the LGM and Heinrich event results is not more marked is unremarkable. Figure 2.10 shows meltwater distribution according to iceberg provenance. Surprising is the extent of Fennoscandian and Icelandic iceberg drift. But, inspection of the L G M surface currents in Figure 2.2 shows that for the input CCM ocean Nordic icebergs would have taken a relatively quick, direct route to the lower latitudes of the North Atlantic. The modelled North Atlantic gyre appears to reach somewhat more northerly than for the real modern ocean. To indicate the true extent of modelled meltwater distribution, which is obscured by the strong North American signal using the scale of Figure 2.9, Figure 2.11 shows LGM with Heinrich flux results by provenance for a five year model run. Mentioned in Chapter 1, it is believed that the North Atlantic Conveyor could have been switched off by the freshwater flux associated with Heinrich event icebergs. Estim-ates place the required increase in flux for shutdown from as little as 0.03 Sv to 0.06 Sv [Fanning and Weaver, 1996]; at 200 km 3/yr (0.006 Sv) the modelled Heinrich flux is well below the expected requirement. Given that the model ocean is static, it is advisable to Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 35 use Heinrich fluxes that would not likely upset the model ocean circulation if the circu-lation and meltwater were coupled. In this light, the assumed Heinrich fluxes seem to be well chosen. In a form more amenable to an oceanographic description, Figure 2.12 indicates the increase in freshwater flux, in Sv/km 2 , averaged over the first five years of the modelled Heinrich event. For modelled Heinrich fluxes other than used for these runs the pattern of Figure 2.12 would not change while the meltwater quantity would scale in proportion to the ratio of the fluxes after steady-state had been reached. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 36 Equivalent Ice Thickness: 1,000 yr, L G M w/ Heinrich Equivalent ice column depth (m) Figure 2.8: Unmelted ice distribution for L GM conditions with Heinrich flux from Hudson Strait. The results shown are for a 1,000 year model run beginning 21,000 yr B.P. Meltwater/Sediment Distribution: 1,000 yr, L G M w/ Heinrich Equivalent ice column depth (m) Figure 2.9: Melted ice distribution for LGM conditions with Heinrich flux from Hudson Strait. The results shown are for a 1,000 year model run beginning 21,000 yr B.P. Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 37 Figure 2.10: Melted ice distribution by provenance for LGM conditions with Heinrich flux from Hudson Strait. The results shown are for a 1,000 year model run beginning 21,000 yr B.P. The scale of Figure 2.9 applies. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 38 Figure 2.11: Melted ice distribution by provenance for LGM conditions with Heinrich flux from Hudson Strait. The results shown are for a 5 year model run beginning 21,000 yr B.P. The illustration is meant to indicate distribution patterns, so a scale is not shown. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 39 Increased Meltwater Flux, 5 year average : Heinrich - LGM Figure 2.12: Increased meltwater flux due to a Heinrich Event. The results shown are the average for a five year model run at the start of a Heinrich Event 21,000 yr B.P. Chapter 2. A Continuum Model of Iceberg Drift: AppUcation to Heinrich Events 40 2.5.3 Modern 1,000 Year Run Plotted in Figures 2.13 and 2.14 are the ice and meltwater distribution for modern con-ditions produced by a model run starting 2,000 years B.P. and lasting 1,000 years. With a relatively vigorous ocean circulation and warm temperatures, modelled modern ice-berg drift is far less extensive than for the L G M ocean. Note from Table 2.2 that the currents of the modelled present ocean require a decrease in the time step, by the CFL criterion, from the LGM runs. Pockets of meltwater near Greenland and Iceland indicate the relatively warm ocean and atmosphere compared to LGM conditions. The modern model run should give an indication of how well the model works, by comparison with, present day iceberg sightings. Plotted in Figure 2.15 are the locations of approximately 600 icebergs sighted by the IIP between 1994 and 1998. The plotted results of Figure 2.14 appear to predict surprisingly well the extent of modern iceberg drift, inspiring confidence in the continuum drift model and the CCM data. A series of other modern runs was conducted to test the role of the various melting mechanisms. It was found, in accordance with [Robe, 1980], that the removal of solar radiative melting indeed had practically no discernable effect on the meltwater pattern of Figure 2.14. In fact, using only the IIP turbulent melting term in the formulation of melting (see Appendix D) yielded almost the same meltwater distribution as when all the melting mechanisms were active. This somewhat substantiates laboratory findings [Huppert, 1980; Russell-Head, 1980] that turbulent processes are dominant in iceberg melting. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 41 Equivalent Ice Thickness: 1,000 yr, Present 0 0 . 0 5 0.1 0 . 1 5 0.2 0 . 2 5 0 3 0 . 3 5 QA Equivalent ice column depth (m) Figure 2.13: Unmelted ice distribution for modern ocean conditions. The results shown are for a 1,000 year model run beginning 2,000 yr B.P. Meltwater/Sediment Distribution: 1,000 yr, Present 5 1 0 1 5 2 0 2 5 3 0 3 5 Equivalent ice column depth (m) Figure 2.14: Melted ice distribution for modern ocean conditions. The results shown are for a 1,000 year model run beginning 2,000 yr B.P. Chapter 2. A Continuum Model of Iceberg Drift: Apphcation to Heinrich Events 42 International Ice Patrol Iceberg Sightings: 1994-1998 20°W Figure 2.15: U.S. Coast Guard-International Ice Patrol iceberg sightings from 1994 1° 1998. The locations of roughly 600 icebergs are shown. Files of the plotted data were obtained from the U.S. Coast Guard [U.S. C.G, 1999]. 2.6 Dri l l Core Prediction The iceberg drift model incorporates a dri l l core prediction routine for ice rafted sediment. The model predicts dri l l core depth and provenance labeled fractional composition. Core stratification is output at user specified sampling intervals. Dr i l l locations for gathering the synthetic cores are also taken as user input; any number of locations may be sampled. (See Appendix D.) Dri l l ing programs have yielded useful information on iceberg sediment content and Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 43 sedimentation rates which can be used to tune the model core prediction routine. For example, Figure 2.16 shows contours of Heinrich event H2 sediment layer thickness based on the evidence of more than 50 marine cores [Dowdeswell et al., 1995]; the H2 sediment distribution is similar to that of HI. The iceberg sediment content inferred in Dowdeswell et al. [1995] is based on the sediment being held mainly in a thin zone at the bottom of the iceberg, which does not translate directly for use in the continuum drift model. Nonetheless, Figure 2.16 gives an indication of the orders of sediment depths the model should predict and the general distribution pattern. For a 12,500 year model run beginning 23,000 years B.P., two synthetic drill cores were produced. Model sampling was carried out at 10 year intervals with input sediment content as listed in Table 2.2. The model drill locations are shown in Figure 2.17. The first synthetic core, shown in Figure 2.18, corresponds to Elliot et al., [1998] core SU 90-08, the record of which is shown in Figure 2.19. The second synthetic core, given in Figure 2.20, was taken at the location of Heinrich [1988] core Me 69-19, for which Figure 2.21 shows the record. The synthetic cores are presented in terms of the fractional composition of provenance-labeled sediment, while analyses of real marine cores are standardly given in terms of mineral or organic material composition. However, direct links can be drawn between abrupt increases in the proportion of North American sediment and Heinrich events IRD peaks in the marine core record. The general shape of the North American sediment plot of the synthetic core shown in Figure 2.18, indicating HI and H2, is broadly in accord with the coarse sediment peaks of core SU 90-08 seen in Figure 2.19. Moreover, Elliot et al. [1998] suggest an average sedimentation rate for core SU 90-08 of 6 cm/kyr, which agrees well with the predicted core. The model core shown in Figure 2.20 over predicts, by about a factor of two, the sediment depth between Heinrich events HI and H2 shown in Figure 2.21. But again, the synthetic North American sediment content and actual Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 44 IRD records show the same general shape. Also, referring to the HI peak in Figure 2.20, the North American sediment makes up approximately 30% of the total sediment content over a layer thickness of about 1 m. If the Heinrich iceberg flux, at 200 km 3/yr, makes up about 30% of the total North American flux, then the Heinrich event HI layer is roughly 9 cm deep. For drill location 2 shown in Figure 2.17 this agrees well with the distribution shown in Figure 2.16. This suggests that an abrupt increase in the modelled North American sediment content is valid as a proxy for Heinrich event IRD. Cores were also predicted for locations near the southern tip of Greenland. These synthetic cores, which are not shown, compared poorly with real drill records, showing very little, if any, North American sediment content. Figure 2.16 and the current lit-erature suggest this is not the actual case. This discrepancy could be attributed to a melting scheme in the model that is too conservative. It is also very likely that proximal to Hudson Strait the deposited sediments were stirred up at the ice stream terminus and transported while suspended in water, rather than being ice-rafted. Though, in general, the modelled Heinrich event H2 meltwater distribution shown in Figure 2.9 successfully predicts the observed sediment distribution of Figure 2.16. Referring to the core of Figures 2.18 and 2.20 as well as the provenance-labeled melt-water distribution of Figure 2.10, the model can readily be adapted to solve for sources within each provenance. If we wish to examine, say, the Hudson Bay and Strait sedi-mentary signal it would be desirable in the model results to distinguish Hudson sediment from that of the rest of North America. However, adding provenances or sub-provenaces to the model linearly increases the model execution time. This caveat likewise applies to the addition of iceberg size categories. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 45 80°N ' ! 1 - i 60°W 40"W 20°W Figure 2.16: Heinrich event H2 sediment distribution [Dowdeswell et ai, 1995]. The contours represent the depth in centimetres of the H2 layer. Marine Core Sample Locations Figure 2.17: Drill locations for the results of Figures 2.18 and 2.20. Location 1 roughly corresponds to core SU 90-08 (43° 31% 30°24'W) from Elliot et al. [1998], while location 2 roughly corresponds to Heinrich [1988] core Me 69-19 (47*19% 19°41'W). Chapter 2. A Continuum Model of Iceberg Drift: AppHcation to Heinrich Events 46 0 . 5 5 , — O S \ 0 . 4 5 -c= 0.4 -g. 0.35 <| 0.3 -o — . tt a 0 . 2 5 L l _ 0.2 -0 . 1 5 -O . l -C o r e s a m p l e at 1 ( 4 2 . 6 8 ' N , 3 0 . 9 4 0.2 W ) D e p t h ( m ) 0 . 4 o .e f s l . A m e r i c a G r e e n l a n d F e n n o s c a n d l a I c e l a n d H 1 O . O S l 16 18 20 T i m e of D e p o s t l o n (kyr B . P . ) Figure 2.18: Synthetic marine core corresponding to Elliot et al. [1998] core SU 90-08 (43°31'N, 30°24'W). The core location is indicated in Figure 2.17. Abo indicated on the figure are the North American sediment peaks corresponding to Heinrich events HO, HI and H2. The core was generated by a 12,500 year run starting at 23,000 years B.P. with a 10 year sampling period. Figure 2.19: Marine core SU 90-08 (43°31% 30°24'W) [Elliot et al, 1998] with Heinrich events HI and H2 indicated. Refer to Figure 2.17 for the drill location. The coarse sediment content, indicated by the number oflithic grains per g, is a measure of ice-rafted sediment content. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 47 0.55 0.5 0.45 0.4 0.35 0.3 ._! 0-25 0.2 O. IS O.I 0.05 Core sample al 2 ( 48.24'N, 1 9.69'W) Depth(m) C2 0.4 0.6 O.B I M . A m e r i c a G r e e n l a n d F e n n o s c a n d l a I c e l a n d ES 16 18 20 Time of Depostlon (kyr B.R.) Figure 2.20: Synthetic marine core corresponding to Heinrich [1988] core Me 69-19 (4T19'N, 19°41'W). The core location is indicated in Figure 2.17. Also indicated on the figure are the North American sediment peaks corresponding to Heinrich events HO, HI and H2. The core was generated by a 12,500 year run starting at 23,000 years B.P. with a 10 year sampling period. C o r e IVIe 6 9 - 1 9 D e p t h <rr>) Figure 2.21: Marine core Me 69-19 (4T19% 19°41'W) [Heinrich, 1988] with Heinrich events HI to H5 indicated. Refer to Figure 2.17 for the drill location. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 48 2.7 Conclusions The pattern of modern iceberg drift as predicted by the continuum model agrees well with observations. More notable, though, are the extents of modelled L G M and Heinrich event iceberg sediment deposits. These model results are in accord with marine records that show Laurentide ice-rafted depositions in the proximity of the Iberian Peninsula [Heinrich, 1988; Dowdeswell et al, 1995 ]. The extensive drift pattern of L G M icebergs is attributed to a relatively cool ocean and atmosphere, promoting iceberg longevity, compared to modern conditions. For the results shown, the iceberg drift and ocean-atmosphere models were not truly coupled; but, still we can speculate about the outcome had there been feedback between these models. From the perspective of location, Fennoscandian and Icelandic sources are best poised to disrupt the downward vertical flow of the North Atlantic Conveyor. However, as contributors to the total North Atlantic meltwater flux these Nordic sources are relatively minor. Conversely, with estimates of Heinrich event fluxes of more than 2000 km 3/yr (0.06 Sv) [MacAyeal, 1993], North American sources, in particular that of Hudson Strait, could have had the requisite calving rates to drastically alter the ocean circulation. Modeled LGM distribution patterns show, however, that North American icebergs did not significantly melt in the high latitudes where the North Atlantic conveyor would be most sensitive. However, the above inferences are drawn from what is known of the modern ther-mohaline circulation. There is reason to believe that just as the cold conditions reached further south for the LGM ocean than at present, as the input ocean temperature fields indicate, the downward flow of the North Atlantic Conveyor might have also shifted south [Lyle, 1997]. In this scenario, Heinrich event icebergs would have been well placed to diminish or shut down the North Atlantic Conveyor. Chapter 2. A Continuum Model of Iceberg Drift: Application to Heinrich Events 49 2.8 Future Developments The foremost task that should be carried out is the coupling of the iceberg drift model with ice sheet and ocean-atmosphere circulation models. Icebergs are essentially the ligature connecting the cryosphere and ocean-atmosphere systems. The development of model coupling these systems would almost certainly constitute a tremendous leap in climate modelling. The triggering of the Younger Dryas by iceberg melt, for instance, could then tridy be investigated. As for modern climatic phenomenon, such a model would allow the investigation of, say, possible links between the modern cryosphere and E l Nino events. With the addition of the continuum iceberg drift model just introduced, all the components to create such an encompassing climate model currently exist. The continuum iceberg drift model can be used over any region of the global ocean, or even over the whole ocean. Given, as an example, synthetic or actual calving data for Antarctic ice sources, the continuum drift model could readily be used in the investigation of iceberg drift and sedimentation in the southern ocean. References Armacost, R. L., 1995. Analysis of the IIP Iceberg Drift Deterioration Model, Re-port CG-D-28-95, U.S. Coast Guard Development Research Center, Groton, Connecticut Bintanja, R. and J. Oerlemans, 1996. The effect of reduced ocean overturning on the climate of the last glacial maximum, Clim. Dyn., 12, 523-533. Bond, G., and R. Lotti., 1995. Iceberg discharges into the North Atlantic on mil-lenial timescales during the last glaciation, Science, 267, 1005-1010. Broecker, W. S., 1994. Massive iceberg discharges as triggers for global climate change, Nature, 372, 421-424. Broecker, W. S., G. Bond, J. Mc Manus, M . Klas and E. Clark, 1992. Origin of the North Atlantic's Heinrich events, Clim. Dyn., 6, 265-273. Budd, W. F., T. H. Jacka and V. I. Morgan, 1980. Antarctic iceberg melt rates derived from size distributions and movement rates, Ann. Glaciol., 1, 103-112. Chappellaz, J., T. Blunier, D. Raynaud, J. M. Barnola, J. Schwander and B. Stauffer, 1993. Synchronous changes in atmospheric C H 4 and Greenland cli-mate between 40 and 80 kyr BP, Nature, 366, 443-445. Clarke, G. K. C , S. J. Marshall, C. Hillaire-Marcel, G. Bilodeau and C. Veiga-Pires, 1999. A glaciological perspective on Heinrich events, In Keigwin, L., R. Webb and P. U. Clark (eds) Mechanisms of millenial-scale global climate change, American Geophysical Union, Washington D.C Dansgaard W., S.J. Johnsen, H.B. Clausen et al., 1993. Evidence for general in-stability of past climate from a 250-kyr ice-core record, Nature, 364, 218-220. Dowdeswell, J . A., 1995. Comment on: "Deep Pleistocene iceberg plowmarks on the Yermak Plateau: sidescan and 3.5 kHz evidence for thick calving ice fronts and a possible marine ice sheet in the Arctic Ocean", Geology, 23, 476-477. Dowdeswell, J. A., A. Elverfi0i and R. Spielhagen, 1998. Glacimarine sedimentary processes and facies on the North Atlantic Margins, Quat. Sci. Rev., 17, 243-272. Dowdeswell, J . A., M . A. Maslin, J. T.Andrews and I. N. McCave, 1995. Iceberg production, debris rafting and the extent and thicknesses of "Heinrich layers" (H-l, H-2) in North Atlantic sediments, Geology 23, 301-304. Dowdeswell, J. A. and T. Murray, 1990. Modelling rates of sedimentation from ice-bergs, Glacimarine Environments: Processes and Sediments. Geological Society London, Special Publication 53, 121-137. 50 References 51 Dowdeswell, J. A., R. J. Whittington and R. Hodgkins, 1992. The sizes, frequen-cies and freeboards of East Greenland icebergs observed using ship radar and sextant, J. Geophys. Res., 97(C8), 3515-3528. Duxbury, A. C , 1971. The Earth and its Oceans, Addison-Wesley, Reading, Mass-achusetts Dyson, J. L., 1972. The World of Ice, Alfred A. Knopf, New York Elliot, M. , L. Labeyrie, G. Bond, E. Cortijo, J. L. Turon, N. Tisnerat and J. C. Duplessy, 1998. Millenial-scale discharges in the Irminger Basin during the last glacial period: Relationship eith the Heinrich events and environmental settings, Paleoceanography, 13, 443-446. Fanning, A. F. and A. J. Weaver. 1996. Temporal-geographical meltwater influences on the North Atlanic conveyor: implications for the Younger-Dryas, Paleocean-ography, submitted Flato, G. M. and W. D. Hibler III, 1992. Modeling Pack Ice as a Cavitating Fluid, J. Phys. Oceanogr., 22, 626-651. Grimm, E. C., G. L. Jacobson Jr., W. A. Watts, B. C. S. Hanson and K. A. Maasch, 1992. A 50,000-Year Record of Climate Oscillations from Florida and Its Temporal Correlation with the Heinrich Events, Science, 261, 198-200. Gwiazda, R. H., S. R. Heimming and W. S. Broecker, 1996. Tracking the sources of icebergs with lead-isotopes: The provenance of ice-rafted debris in Heinrich layer 2, Paleoceanography, 11, 77-93. Haupt, B. J., C. Schafer-Neth and K. Statteger, 1994. Modeling sediment drifts: A coupled oceanic circulation-sedimentation model of the northern North At-lantic, Paleoceanography, 9, 897-916. Heinrich, H., 1988. Origin and Consequence of Cyclic Ice Rafting in the Northeast Atlantic during the Past 130,00 Years, Quat. Res., 29, 142-152. Hibler III, W. D., 1979. A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815-846. Hibler III, W. D. and K. Bryan, 1987. A Diagnostic Ice-Ocean Model, J. Phys. Oceanogr., 17, 987-1015. Huppert, H. E. and J. S. Turner, 1978. On melting icebergs, Nature, 271, 46-48. Huppert, H. E. and E. G. Josberger, 1980. The Melting of Ice in Cold Stratified Water, J. Phys. Oceanogr., 10, 953-960. Josberger, E. G. and S. Neshyba, 1980. Iceberg melt-driven converction inferred from field measurements of temperature, Ann. Glaciol., 1, 113-117. Lyle, M. , 1997.Could early Cenozoic thermohaline circulation have warmed the poles?, Paleoceanography, 12, 161-167. MacAyeal, D. R., 1993. Binge/purge oscillations of the Laurentide Ice Sheet as a cause of the North Atlantic's Heinrich Events, Paleoceanography, 8, 767-773. References 52 Manabe, S. and R. J. Stouffer, 1993. Century-scale effects of increased atmospheric C O 2 on the ocean atmosphere system, Nature, 364, 215-218. Marshall, S. J., 1996. Modelling Laurentide ice stream thermomechanics, Ph.D. thesis, Univ. of B.C., Vancouver, B.C., Canada Marshall, S. J. and G. K. C. Clarke, 1997. A continuum model of ice stream ther-momechanics in the Laurentide Ice Sheet 2. Application to the Hudson Ice Stream, J. Geophys. Res., 102(B9), 20,615-20,637. Matsumoto, K., 1996. An iceberg drift and decay model to compute the ice-rafted debris and iceberg meltwater flux: Application to the North Atlantic, Pale-oceanography, 11, 729-742. Matsumoto, K., 1997. Modeled glacial North Atlantic ice-rafted debris pattern and its sensitivity to various boundary conditions, Paleoceanography, 12, 271-280. Mountain, D. G., 1980. On predicting iceberg drift, Cold Reg. Sci. and Tech., 1, 273-282. Napoleoni, J. G. P., 1979. The dynamics of iceberg drift, M.Sc. thesis, Univ. of B.C., Vancouver, B.C., Canada Parkinson, C. L. and W. M. Washington, 1979. A large-scale numerical model of sea ice, J. Geophys. Res., 84(C1), 311-337. Patankar, S. V., 1980. Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York Paterson, W. S. B., 1994. The Physics of Glaciers, 3rd ed., Elsevier Sci., New York Pfirman, S. L., R. Colony, D. Niirnberg, H. Eicken and I. Rigor, 1997. Recon-structing the origin and trajectory of drifting Arctic sea ice, J. Geophys. Res., 102(C6), 12,575-12,586. Pickard, G. L., 1964. Descriptive Physical Oceanography, Pergamon Press, New York Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, 1986. Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York Robe, R. Q., 1980. Iceberg drift and deterioration, Dynamics of Snow and Ice, Academic Press Inc., San Diego, 211-259. Russell-Head, D. S., 1980. The melting of free-drifting icebergs, Annals of Glacio-logy, 1, 119-122. Sodhi, D. S. and M. El-Tahan, 1980. Prediction of an iceberg trajectory during a storm, Ann. of Glacio., 1, 77-82. Thorndike, A. S. and R. Colony, 1982. Sea Ice Motion in Response to Geostrophic Winds, J. Geophys. Res., 87(C8), 5845-5852. Thorndike, A. S., D. A. Rothrock, G. A. Maykut and R. Colony, 1975. The Thick-ness Distribution of Sea Ice, J. Geophys. Res., 80(33), 4501-4513. References 53 U.S.C.G., 1999. www.uscg.mil/lantrea/iip, United States Coast Guard - Interna-tional Ice Patrol web page Wang, H. F. and M. P. Anderson, 1982. Introduction to groundwater modeling, W. H. Freeman and Company, San Francisco Weaver, A. J., M . Eby, A. F. Fanning and E. Wiebe, 1998. Simulated influence of carbon dioxide, orbital forcing and ice sheets on the climate of the Last Glacial Maximum, Nature, 394, 847-853. Weaver, A. J. and T. M. C. Hughes, 1996. On the incompatibility of ocean and atmosphere models and the need for flux adjustments, Glim. Dyn., 1 2 , 141-170. Weeks, W. F. and W. J. Campbell, 1973. CRREL Research Report 200, U.S. Army Cold Regions Engineering Research Laboratory, Hanover, New Hampshire Appendix A FORMULATION OF THE GOVERNING EQUATIONS The following presents the physical arguments leading to the formulation of a continuum discrete finite-difference formulation. A . l The Continuous Model There are two means of treating the icebergs in the model; they can be dealt with in terms of their actual number or equivalently expressed as the corresponding volume of ice. The model discussed herein uses the latter approach largely because of practical concerns, as glacial ice flux data is given in volumetric terms. As well, since this is essentially a probabilistic model, or perhaps quasi-deterministic, it stands to reason that the mathematical definition of the natural phenomenon should readily describe its spread. To this end, working with ice volume is is more suitable than defining the number of individuals. A.1.1 Volume Density of Icebergs Given an ocean surface area A, the volume of icebergs present as a function of time is given by model of iceberg drift and deterioration which, in turn, constitutes the basis for the V(t) 54 Appendix A. Formulation of the Governing Equations 55 where v is the local volume of icebergs per unit area. In a spherical coordinate system r denotes position (R, i?,A), where R RS 6371km is the Earth's radius and i? and A are longitude and co-latitude, respectively. . Clearly, variously-sized icebergs will drift differentially according to the forcings acting on them determined by, for instance, the number of ocean layers in which they sink. The model, therefore, must account for the spectrum of iceberg sizes, such that Jo where s is a size parameter such as volume, sail height or waterline length for example. If it is assumed that there is an finite upper limit Smax on iceberg size, ice volume is then J A Jo A.1.2 Volume Production by Sources Let $y represent iceberg volumetric flux such that $y is source production, describing calving at a floating glacier terminus, and $y can be viewed as representing iceberg breakup, or conversely, production of smaller bergs by breakup. Without a realistic method to model iceberg breakup, only production by glacier calving will be considered. In general A.1.3 Volume Diffusion of Icebergs A viable continuum model of iceberg drift requires that a probabilistic, or stochastic, component be incorporated. The object is, in an efficient and practical manner, to with which the addition of size dependence, as in (A.l), extends to (A.2) Appendix A. Formulation of the Governing Equations 56 account for various individual iceberg properties, like shape, that affect drift. If 3v(s, ri t) denotes volume diffusion flux of icebergs of size s, then following a simple diffusion law 3v(s,r,t) = -I)jk(s,T,t)Vv(s,T,t) where ~Djk is a tensorial diffusion coefficient such that in a local frame of reference diffusion differs along and transverse to the iceberg trajectory D DA 0 0 DT Determination of the stochastic diffusion coefficients will be specifically addressed in Appendix B. Iceberg volume flux as a function of diffusion is then dVi dt JA JO V • 3v(s,r,t)ds dA. (A.3) A.1.4 Reynolds Transport: Volume Advection of Icebergs Combining (A.2) and (A.3) yields an expression for global volume balance -= / f dt JA JO $+(.*, r, t) - V • t)jk(s, r, t)Vv(s, r, t) ds dA. (AA) Furthermore, Reynolds transport theorem allows the left-hand side of (A.4) to be written as dV ~dt fl flSmax £} r -| = J J -v(s,r,t) + V-[iJ{s,r,tWs,T,t)\dsdA. (A.5) where U = (U^,V\) is a velocity vector. Thus the local volume balance expression, by equating (A.4) and (A.5), is 8 dt v(s,r,t) + V- lJ(s,T,t)v(s,T,t) =fy{s,T,t) + V- D jk{s, T , t)Vv(s, r, t) . (A.6) Appendix A. Formulation of the Governing Equations 57 A.1.5 Melting of Icebergs By inspection, (A.6) is somewhat generic in that it could describe a host of natural phe-nomena; soil contaminant transport or general heat and fluid flow, for example. Icebergs, though, are also subject to volume loss by melting, requiring that another term be added to the local balance equation. Given an average volumetric melt rate — rh(s,r,t) of an iceberg of size s and its inverse specific volume p(s) (Volume -1) the global expression for melting is dV„, dt IA lo rh(s, r, t)p(s)v(s, r, t) ds dA Incorporating (A.7) into (A.6) gives a final equation dt" v(s,r,t) + V- V(s,T,t)v(s,r,t) -m(8,T,t)n(s)v(s,r,t) l + ( a , r, t) + V • [Bjkis, r, t)V v(s, r, t) which constitutes the general continuum description of iceberg drift. (A.7) (A.8) A.2 The Discretized Model Applied to a spherical coordinate system, rearranging and expanding of (A.8) gives 3 . | 1 dt Rsin X dv ' d ( D* dv\ d , . ^ a* l^d* +d\\DxsmXdx = $ y + rhpv (A.9) R2 sin A Obviously, d/dR terms vanish and are excluded. An approximate solution for (A.9) is found using finite-difference methods over a discretized North Atlantic Ocean on the surface of a two dimensional spherical grid, as shown in Figure A . l Appendix A. Formulation of the Governing Equations 58 Figure A . l : Discretized North Atlantic problem domain. The discretization uses a 60 x 30 spherical grid with 1.875° zonal and 1.855° meridional spacing. The problem domain spans the area from 75° W to 75°E and between 28° N to 85° N. The methods of finite-difference calculus will not be dealt with at great length here because they are well described elsewhere. For the uninitiated reader, the principles of finite-difference techniques are outlined in Patankar [1980] and Press [1986]. Before delving into the finite-difference equations it is instructive to rewrite (A.9) in a form that is more amenable to the description of probabilistic iceberg drift. Using L and T as general length and time units, consider a single cell on the discretized ocean of Figure A . l ; let v = pip, where p is the inverse cell surface area in units (L~2) and ip is the actual volume of ice located within the cell. Remember that the diffusion is purely stochastic, with no relation to any change in potential. Its units are not as expected from Appendix A. Formulation of the Governing Equations 59 more usual flow problems, such as soil diffusivity with units of ( L 2 T _ 1 ) representing the flux across a control volume surface due to a change in potential. Rather, in this case diffusion represents drift of ice astray from a mean drift trajectory, so that the diffusion units are in terms of distance wayward drift per distance mean drift per unit time or (LL^T'1 = T - 1 ) . This leads to the reformulation of (A.9) as 1 o l i p 4 , ) + RsinX d ( D# dVA d ( . Bib = $+ + mp(pib). (A.10) R2 sin A which is simplified at a discrete cell by integrating over the surface area '+- f^i+-A = R2 3 " 2 sin A dtidX A . 1 JS._! to get {U#pib)# x - (U#pib)# AA,-+ R (sin A V\pip)\,+ i 7 J A s m A -(sin \Vxptp)\ - | Dx sin A — J 2 A 0tf sm sin A <9i? AA, = Qv+mpib. (A .11) where the square bracketed terms represent fluxes at the cell boundaries located at i? i ± ±. and A-. i . Referring to Figure A.2 on page 67, illustrated is the cell arrangement for establishing the five-point finite-difference scheme used by the model. The transient ice volume ib is defined at the cell centres while the velocities U$ and V\ and diffusion coefHcents, D^ and D\, are staggered by half a grid space in their respective flow directions, to be defined at the cell interfaces. Appendix A. Formulation of the Governing Equations 60 In the spatial finite-difference approximation of (A.11) an upwind differencing scheme is used, with backward difference first derivative approximations, for the advection terms. Upwind differencing for advection as used in this model is outlined in Patankar [1980]. In the model, as in nature, iceberg drift is deterministic based on the forcings from air and water currents and the diffusive nature is characterized by a very low Peclet number. The diffusion in the model is merely a statistical quantity to capture the gross drift pattern of numerous icebergs and not associated with a physical diffusion process, please refer to Appendix B. The diffusion terms are approximated by second-order centred differences. The finite-difference formulation of (A.11) for the centre cell of Figure A.2, summing the fluxes across the East, West, North and South interfaces and considering only icebergs of size 5, is then At + R + R sin Xj sin A,-AA, + sin s , t - l , j + D i sin AA,-i sin A, AA, AA,-Mi = Qt,ij + ms,i,jPsi>s,i,j-(A.12) where the subscripts (s,i,j) refer to Xj). The formulation (A.12) is implicit in time such that As the subscript 5 denotes, different sized icebergs experience different flow fields and have different melting, production and diffusion rates. Thus, (A.12) superposes for the Appendix A. Formulation of the Governing Equations 61 range of discrete iceberg sizes in the model so that: ffl size s=l The calculation of size dependent variables is discussed in Appendix C. The problem domain is bounded around its periphery by a single depth of buffer cells where the Dirichlet boundary conditions have iceberg volumes and ocean velocity fields permanently set to zero. Within the problem domain indices 0 and 1 are used to specify free (ocean) and no-flow (land) cells, respectively. Where a cell is specified as no-flow, with an index of 1, iceberg volumes, diffusion coefficients and velocities are also set to zero. A.2.1 Local Forces and Acceleration: Determining Velocity Fields The finite-difference scheme in the model uses velocity fields determined by iceberg ve-locity drifting dynamics. Governing the drift of icebergs are four main forces: wind drag, water drag, the Coriolis effect and gravitational pull due to sea slope. To start, the general dynamic water and wind drag force expressions are FA = ^ A ^ I U ^ I U X , ( A 1 3 ) 2 and PWCDWAW\^JWI\^WI t \ I A \ tw = (A. 14) where UAI a n d Uwi are the iceberg velocities relative to air and water with general x and y components U and V. CDA and Cow are drag coefficients, PA and pw are densities and AA and Aw are the iceberg projected surface areas with respect to air and water currents. The drag coefficients are constants that quantify the effects of individual Appendix A. Formulation of the Governing Equations 62 iceberg characteristics on the external flow field forcing. They characterize form, friction and iceberg inertia or Reynolds number. Generally they are empirically determined and combined in a single constant, as is the case here. In reality, however, iceberg drag coefficients are orientation dependent. In this model orientation dependence is neglected, as is the case for other models [e.g. Mountain, 1980; Matsumoto, 1996]. Incorporating orientation dependence in this continuum model would add a level of complication that is not justified by any gain in accuracy. Only where details of the motions of a discrete iceberg are being examined might one be interested in the variation of the drag coefficient with orientation. Particular attention is paid to water drag owing to its dominance in governing the drift of icebergs. In a multi-layer ocean where currents at different levels act on the portion of the iceberg within those layers (A. 14) is written as ^ layer FW = ]T Awl\Vwh\\3wil (A.15) i=i where the subscript I denotes ocean layer. The primacy of water drag is not only con-firmed by historical observations of icebergs following ocean currents, but also by recent numerical models of individual iceberg drift such as Matsumoto [1996] showing that water drag is typically of an order of magnitude greater than wind drag. Of lesser importance than wind and water drag, but nonetheless significant in an accurate description, as objects that move slowly over long distances, icebergs are also influenced by Coriolis and sea slope forces. Examining the linear shallow water equations, the components of the acceleration of a particle over a flat bottom due to the Coriolis effect and sea slope are Appendix A. Formulation of the Governing Equations 63 and where / is the Coriolis frequency, h is the sea surface elevation and g is acceleration due to gravity. The input velocity fields in the model, as output by the coupled ocean-atmosphere model of Weaver [1998], are given as the mean annual East-West and North-South com-ponents for last glacial maximum (LGM) and present conditions. It is assumed that these relatively long term averages are essentially geostrophic values. The assumption is made that in a state of geostrophic equilibrium the surface slopes, dh/dx and dh/dy, are associated with the local geostrophic current U Q , of which UG and VG are the com-ponents. On this basis the surface slope pressure gradient force must exactly balance that due to the Coriolis acceleration of the surface geostrophic current. Recall that the Coriolis effect deflects positive x and y motions in the negative y and positive x direc-tions, respectively. Therefore, for positive accelerations due to sea slope we can write the equivalent expressions 9Yx = f V a ( A - 1 8 ) and 9 ~ = -fUG (A.19) Applied back to icebergs, (A.13), (A.15), (A.18) and (A.19) then compose to give dU' dt fV - fVG + er PACDAAA\UAI\UAI + PWCDW ^ AWi\Cwii\Uwii i=i 2Mi (A.20) Appendix A. Formulation of the Governing Equations 64 and dV' dt = -fU + fUG + er PACDAAA\UAI\VAI + pwCDw ^ Aw^wi^Vwii 1=1 2Mi (A.21) where Mj is iceberg mass. The above, (A.20) and (A.21), constitute the x and y equations for iceberg acceleration due to wind drag, water drag, the Coriolis effect and sea slope. In the absence of explicit x and y terms in (A.20) and (A.21), velocity components U' and V are functions of time t only. Given a range of typical values of water velocity and iceberg mass, 0.1ms - 1 and 5 X 106 Tonnes are representative, dimensional analysis on the dominant water drag term shows that the time scale for acceleration of the iceberg from rest to the surrounding fluid velocity is on the order of a few hours at most. This is also confirmed by model runs using individual icebergs to determine the stochastic diffusion coefficients, as dis-cussed in Appendix B. With model time steps of around five days, determined by the Courant-Friedrichs-Lewy (or CFL) criterion for numerical stability, the contribution of the acceleration period towards the net iceberg trajectory during each time step is rel-atively insignificant. This is especially true given that the icebergs in the model are generally not starting from rest; the ocean flow fields being continuous are not expected to yield abrupt iceberg velocity changes between cells. Thus the loss in accuracy in as-suming constant iceberg velocity within each model cell is negligible. Using (A.20) and (A.21) in addition with ~ = U' (A.22) dt v ; and — = V' (A.23) dt Appendix A. Formulation of the Governing Equations 65 to form a complete set of equations of iceberg motion, a fourth-order Runge-Kutta routine uses the input multi-layer ocean and atmosphere velocity fields to establish a velocity field for each size of iceberg in the model. The numerical formulation of (A.20) through (A.21) uses cell-centred input multi-layered ocean velocity fields and calculated Coriolis para-meters in the differential equations to determine cell-centred iceberg velocity fields. Over a cell of the discretized grid for size s icebergs these are written as dU, dt sL=f,<V.iJ-VGiJ,i) + N layer PAAAs[(Us.tj - UGi,jfi)2 + (Vs-. - VG^oflHU + Pwf2 Aw.^,i[{U'Sti- - UGiJtl)2 + (V;tiJ - VGi,j:l)2}Hu'3:itj - UGiij,t) UGi,j,o) 1 2 M 7 (A.24) and dV. dt P A A A S [([/;.. - uGi,jfiy + (I/;.,,. - vGitj,0)2} Hv;tij - vGiJt0) Nlay er pw 2M 7 (A.25) where the subscripts 0 and 1 on the geostrophic velocities imply air and surface water layers, respectively, and the drag coefficients are taken to be one. The assumed drag coefficients are corrected by incorporating diffusion into the model. To ensure continuity the iceberg velocities in (A.12) are defined in the continuous model at the cell interfaces through taking the harmonic average of adjacent cell-centred discrete iceberg velocities as follows U-<^ = K , ^ W + U : M J A A • ( A - 2 6 ) To determine the cell constant velocities at the start of a run the icebergs are initially taken to be at rest. When ocean conditions switch from LGM to present in a model run, iceberg velocity fields are calculated anew. Appendix A. Formulation of the Governing Equations 66 A.2.2 Bathymetry In the model, iceberg draft (depth below water) is a function of iceberg volume, as discussed in Appendix C. Since the local ocean depths from the model of Weaver [1998] are input fields this allows the model to include depth restrictions on iceberg drift. Where an iceberg of size s has a greater draft than the local depth its velocities USiij and Vs,i,j and diffusion coefficients DStij are set to zero. A.2.3 Two-Dimensional Coordinates It is somewhat desirable to work in two-dimensional coordinates owing to their relative simplicity compared to spherical coordinates. Since iceberg drift occurs on the surface of a sphere, the switch to two-dimensional, quasi-Cartesian, surface coordinates is readily facilitated by the conversions dx = Rsin Xd-d and dy = RdX which allow (A.12) to be written as Ax, Ayj + D-i Aa + Ax-Ay; Ay; = Qt.ij + in-wPd'.id-Ax D fps,i,j -A % (A.27) Please refer to Figure A.3 for details on the Cartesian grid corresponding to the spherical grid of Figure A.2. The time discretization and solution of the system of equations is deferred until Appendix E. Appendix A. Formulation of the Governing Equations Figure A.2: Five-point finite-difference cell scheme on a spherical grid. Appendix A. Formulation of the Governing Equations Figure A.3: Spherical grid: two-dimensional surface coordinates equivalent. Appendix B S T O C H A S T I C D I F F U S I O N B . l Theory Observations tell us that icebergs initially in close proximity to each other will diverge over time when acted on by external forcings in the form of wind and water drag. Sup-posing that the icebergs are close enough to each other such that they experience the same currents, it may be assumed that the different reactions are due to the individual iceberg parameters affecting drift, namely size and shape. Iceberg tracking data and previous numerical models show that icebergs closely follow the ocean currents because the forces driving them are dominated by water drag. Scale analyses of the local drift equations also reveal that in steady flows, icebergs of the sizes considered in this model reach constant velocity in a matter of a few hours, generally much less. In the finite-difference model we wish to investigate the long term, gross drift and sediment deposition patterns using constant geostrophic ocean and air velocities at each grid block with time steps of several days. The grid block dimensions are roughly 200 km North to South and typically around 100 km East to West. Under these conditions it is computationally expensive and impractical to consider the small time scale drifting physics, it negates the expedience of using a continuum model, especially if the local accelerations are of relatively minor importance in the total drift across the grid blocks over the time scales considered. Nonetheless, it would clearly be negligent to disregard the local drift physics. It is the variable forcings on the individual icebergs that cause them 69 Appendix B. Stochastic Diffusion 70 to spread from one another, playing an important role in determining the gross pattern of ice rafted deposition. As a means of accounting for the variable drift in each grid block the model includes a stochastic component describing the dispersion of icebergs through the range of expected sizes and shapes due to the variation of individual characteristics. B . l . l Velocity Dependence Iceberg scatter is a function of the forcings they experience; where they are not subject to external forces, in the form of velocity fields, they will not scatter. To account for this the diffusion law is modified by (B.l) where U and V are the local cell iceberg velocity components, U m a x is the largest ex-pected iceberg velocity giving maximum diffusion and n is an experimentally determined parameter describing how the diffusion varies with iceberg velocity. B.l .2 Coriolis Dependence The effect of Coriolis acceleration is to deflect icebergs from a straight-line that wind and water drag forces in equilibrium alone would determine. As the Coriolis frequency increases with latitude, it is conjectured that the rate of diffusion also thus varies. Drifting objects at the equator, where the Coriolis force vanishes, would not experience any scatter due to the Coriolis effect. Taking Coriolis acceleration into consideration, the diffusion law is modified by A = G(Xj) (B.2) where G is a function of the discretized latitude Xj. Appendix B. Stochastic Diffusion 71 B . l . 3 Sea-Ice Gating The drift of icebergs is constrained by the presence of level sea-ice inasmuch as it acts a barrier. To capture the effects of sea-ice as it hems in icebergs the model uses a global diffusion matrix of the factor in each cell by which iceberg scatter is impeded. For example, for an M by N cell problem domain, such as illustrated in Figure A . l on page 58, the global diffusion is f hi ••• kiN \ K (B.3) \ &M1 • ' • kMN J where 0 < kij < 1. Quite simply, each cell in the discretized ocean is assigned a user input value between zero and one which specifies the factor by which diffusion is affected by the presence of sea-ice; in a cell unencumbered by sea-ice hj — 1, where sea-ice completely prevents iceberg scattering hj — 0. B . l . 4 Formulation of Diffusion Coefficients To incorporate the effects of velocity and sea-ice gating into the mathematical description of diffusion, first recall the diffusion tensor from Appendix A, Section 1.3, defined in a local coordinate system ( DA 0 \ D = (B.4) V 0 DT ) where DA and DT are along and cross-stream diffusion coefficients. Let 7 be the rota-tion angle of the local coordinate system, measured counter-clockwise, from the global coordinate system where U 7 = cos 1 (U2 + V2)* Appendix B. Stochastic Diffusion 72 and define the rotation matrix R as . cos 7 sin 7 . R = • (B-5) — sin 7 cos 7 Using (B.4) and (B.5) in the identity D ' = R - D R = X X (B.6) 0 Dyy transforms D into the global coordinate system. At cell accounting for the local iceberg velocity dependence from (B.l), Coriolis dependence from (B.2) and the local sea-ice gating from (B.3) gives the final diffusion coefficients (D'L o \ D" = % A C D ' = (B.7) V 0 Dvy) In the model the input local along- and cross-stream diffusion coefficients of (B.4) are cell-centred. Operating on (B.4) using the rotation matrix (B.5) and accounting for local iceberg velocity and sea-ice gating as in (B.7) gives the final cell-centred x and y diffusion coefficients, hereafter written as Dx and Dy. To ensure continuity and numerical stability in the diffusive flow the diffusion coefficients are staggered to be defined at the cell interfaces, as in Figure A.2. This is accomplished through taking the harmonic average of the diffusion coefficients of adjoining cells, accounting for variabilty in cell spacing. For instance, given adjacent cells (i,j) and {i + l,j) with x direction cell-centred diffusion coefficients DXiij and D X | ; + i j , the cell-interface diffusion coefficient is calculated as n _ {Dx,s,i,jDXiSii+ltj)(Axiij + Axi+lij) U'>i+12>i ~~ _ n • . A T - , -4- D • - A T - • ^ ' As discussed in Appendix A and Section B . l , each size of iceberg in the model is dy-namically unique. Accordingly, the diffusion process in the model is also size dependent, denoted by subscript s, as well as a function of position and time. Appendix B. Stochastic Diffusion 73 B .2 Experimental Determination of the Diffusion Coefficients The following outlines the experimental method to determine the tensorial diffusion co-efficients used in the stochastic component of the iceberg drift model. A model of discrete, or individual, iceberg drift has been developed using the equations of motion (A.20), (A.21) (A.22) and (A.23) in a fourth-order Runge-Kutta routine with variable time stepping. Having chosen average iceberg mass, drag coefficients and a uniform ocean velocity field this model is used to establish a mean iceberg trajectory over the continuum model time step. The mean trajectory is the control condition for the experiment to follow. The force acting on an iceberg is a function of the external velocity field, mass and drag coefficients. Surface area is a factor, but is a function of mass in the numerical model and therefore implicitly included. By the estimate of Russell-Head [1978] the drag coefficients range between 0.6 and 2.0, which are assumed to be random and uniformly distributed. Performing multiple runs of the discrete drift model, holding either the velocity field or mass constant and varying the other over its expected range of values allows comparison against the control run to determine the role of the individual variables in determining drift trajectory. The results of the analyses are used to determine the mean daily drift astray, along and traverse, from the control trajectory and variation with iceberg mass and velocity. B .2.1 Control Condition Tables B . l and B.2 list the iceberg and ocean parameters used to establish, the control trajectory in the experiment. The control values were chosen to represent those charac-teristic of icebergs, the atmosphere and ocean. The experimental flow field comprised single ocean and atmosphere layers. Appendix B. Stochastic Diffusion 74 Parameter Units Value CDW , GDA 1.3 MI tonnes 5 x 106 Shape cylindrical height/width 1 Table B . l : Control trajectory iceberg parameters. Parameter Units Value A °N 55 U W m s- 1 (km day - 1) 0.1157 (10) VW m s- 1 (km day - 1) -0.1157 (-10) uA m s"1 5 VA m s _ 1 -5 Table B.2: Control trajectory flow field parameters. B.2.2 Analyses Using random and uniformly distributed drag coefficients, multiple model iceberg tra-jectories were generated to gauge the effect of mass, Coriolis frequency and flow field velocity on the dispersion of the model icebergs. The drag coefficients in air and water were taken to be the same and ranged between 0.6 and 2.0. For each model run, random perturbations were added to the flow fields to simulate the variability in real ocean and atmosphere currents. The meridional and zonal flow components of water velocities were each allowed to fluctuate within 10% of their mean values over six hour periods. The air velocity components were allowed to fluctuate over a range of 50%. As well, a diurnal tidal current was superposed on all the runs, including the control run. The tidal current velocity magnitude is given by Duxbury [1971] as Uude = (B.9) Appendix B. Stochastic Diffusion 75 Parameter Units Range N •L' runs Duration tonnes Ixl0 6 -200xl0 6 10 1,000 5d f s-1 /(30°N)-/(85°N) 10 1,000 5d U w m s- 1 0.01-0.5 10 1,000 5d Table B.3: Variation in drift parameters. Nruns denotes the number of model trajectories generated in each set (i.e. for each parameter 10,000 runs were performed). where A is the tidal amplitude, g is gravitational acceleration and d is the average water depth. For the purpose at hand, A and d were assumed as 0.5 m and 3,000 m, respectively. The x and y components the diurnal tidal current velocity are written as Ayfgd cos 2~xt Util de (B.10) anc Vtide Ay/gd sin 2wt ( B . l l ) where t is time in days. Holding, for instance, the Coriolis frequency and flow field velocity constant, multiple sets of runs were conducted for various iceberg masses. For example, 1,000 trajectories were generated for a given iceberg mass and then diffusion coefficients were calculated; without changing any other parameters, the mass was incrementally increased and an-other 1,000 trajectories generated and diffusion coefficients were again calculated. In this way, the effect of mass, or another parameter, on the dispersion of icebergs was examined. Table B.3 indicates the ranges over which the parameters varied in the experiment. The ranges chosen represent those expected for real icebergs. Figure B . l shows example positions after three days for 100 model icebergs drifting under the control conditions with the addition of random perturbations and tidal cur-rents. Incidentally, the average of the net trajectories for the multiple perturbed runs Appendix B. Stochastic Diffusion 76 and the net trajectory of the control run were negligibly different. Taking the final posi-tions of each of the probabilistic trajectories and that of the control or mean trajectory, diffusion coefficients for each iceberg were calculated. Referring to Figure B.2, showing sample three day trajectories, the along-stream diffusion coefficients were calculated as where .if OB- OA e = cos \OB\\OA The cross-stream diffusion coefficient was calculated as Dxrun = - — j ^ j j — • (B.13) Note that as formulated in (B.12) and (B.13) the diffusion coefficients are dimensionless, whereas in the governing equation (A.9) from Appendix A the diffusion coefficients are in terms of T _ 1 , where T is general time. The drift coefficients describe the wayward drift as a fraction of the net mean drift. To say, as an example, that an iceberg strays from the mean trajectory at a rate of 5% per day is the same as stating the wayward drift as 5% per week or month. So, the diffusion coefficients as calculated in (B.12) and (B.13) can generically be applied to any time unit. In the classical definition of diffusion this is a invalid approach; diffusion rates of 5 m 3 per day and 5 m 3 per week, for example, have widely different physical implications. For the multiple runs within each set the final or average along-stream diffusion coefficient was calculated as DA = 5 ] \DArun\ (B.14) which was likewise applied to the cross-stream coefficients. Appendix B. Stochastic Diffusion 77 Figure B . l : Iceberg drift variation after three days. The variability in the final positions shown is a function of the random distribution of drift coefficients and perturbations to the mean ocean and air currents with diurnal tidal currents superposed. Appendix B. Stochastic Diffusion 78 Figure B.2: Example drift trajectories to demonstrate the calculation of diffusion coeffi-cients. The line OA represents the net trajectory for the control condition while line OB represents the net trajectory of a typical iceberg for which the diffusion coefficients are calculated. Appendix B. Stochastic Diffusion 79 B.2.3 Results In determining the stochastic diffusion coefficients, the experimental results were some-what contrary to expectations. In the experiment, the iceberg diffusion process seemed remarkably robust in that the variation of parameters had very little bearing on the calculated diffusion rates. For iceberg mass and Coriolis frequency, their variations were found to negligibly, if at all, effect the calculated diffusion coefficients. The exception was over very short run durations, less than a half day, during spin-up from rest. During this initial period the icebergs dispersion varied linearly according to mass and roughly as the sine of latitude. However, after roughly a day the dispersion that took place during the spin-up from rest was found to insignificant contribute to the net trajectory. In any case, except at source grid blocks, the icebergs in the continuum model would not start from rest. As well, over the range of flow field velocities tested, the water velocity was found to have practically no effect except at very low velocities where the diffusion coefficients slightly increased. Even when the mean flow field velocity was near zero, tidal currents and perturbations still resulted in iceberg scattering. It is thought that the modelled increased dispersion at low flow velocity was due to Coriolis forces being relatively strong, deflecting the icebergs from their mostly water drag determined courses. The most salient finding from the experiment was that the along- and cross-stream diffusion coefficients vary nearly linearly as a function of the amount of perturbation to the flow field. In particular, perturbations in the water currents translate almost directly to the diffusion coefficients. The model icebergs proved to be only marginally sensitive to perturbations in the air currents; at velocities below about 15ms - 1 air current effects were almost unnoticable in the iceberg trajectories. Figure B.3 shows plots of the diffusion coefficients as a function of the mean water current perturbation. Based on what the Appendix B. Stochastic Diffusion 80 Di f fus ion coe f f i c ien ts a s a func t i on of p e r t u r b a t i o n 3 0 4 0 5 0 60 7 0 pe r tu rba t i on f r o m m e a n w a t e r ve loc i t y (%) Figure B.3: Iceberg diffusion coefficients as a function of the amount of random perturb-ation from the mean water velocity. modeler judges to be a reasonable amount of perturbation in the real ocean currents, the diffusion coefficients for use in the continuum model can be selected from Figure B.3. The experimental results conveniently allow the Coriolis and velocity dependence terms to be dropped from tensorial diffusion formulation (B.7). From this (B.7) is re-written as D" 0 D" yy (B.15) where, recalling (B.6) D' = R ^ D R Appendix C I C E B E R G SIZE DISTRIBUTION AT SOURCES Since calving sources are in general geographically remote and the bulk of an iceberg lies out of view, the size distribution of icebergs is not well understood nor comprehens-ively documented. Two things are clear: the size distribution of icebergs depends on the source and is skewed towards smaller sizes. A degree of flexibility in the model is necessary to accommodate the difficulties apparent in dealing with such a paucity of information. Provisions have been made for the iceberg size distribution at the sources to be readily modified to fit any appropriate continuous probability density function. For each provenance in the model (N. America, Greenland, Iceland and Fennoscandia) different moments of distribution, namely mean iceberg size and standard deviation, are taken as input parameters. As noted, the full extent of an iceberg is rarely observable. Since ice volume is not directly measurable, some other size parameter must be used as a proxy to infer volume. One such measurable iceberg parameter is waterline length. If in the numerical model some standard iceberg shape is assumed, as is the only viable option, such as a cube or cylinder, then the conversion from the observable quantity to volume is straightforward. To clarify, the following demonstrates how iceberg size distribution is dealt with in the model. As user input, the moments of distribution for each provenance, maximum iceberg size as well as the number of sizes present are supplied to the model. As well, the model reads an array of indices indicating the provenance of the sources at each grid block. 81 Appendix C. Iceberg Size Distribution at Sources 82 Consider, for example, iceberg flux from a Greenlandic source with iceberg waterline length LWL as the observable measure of size. At the risk of oversimplifying, assume a Rayleigh distribution 2aiexp(—x2/x2) p(< x) = ^ '-x for size production of icebergs at their calving sources, where x and x are a given size and the mean, respectively. The Rayleigh distribution does not require a standard deviation, of which field measurements are extremely poor. If Greenland icebergs have a mean waterline length of, say, LWL,GL = 150 m then we get a distribution curve as illustrated in Figure C . l below. ,-3 Rayleigh Distribution: mean = 150 Size Parameter: Waterl ine Length, (m) Figure C . l : Example: Greenland source iceberg size distribution. Given the user input number of size categories in the model, five for this example, the sizes are equally divided along the abissca between zero and a maximum, as shown in Figure C . l . The maximum size is a user input parameter specifying the upper limit Appendix C. Iceberg Size Distribution at Sources 83 of iceberg size to be considered in the model and applies to all provenances. The areas under the curve representing the probability of each size category are integrated using six-point Gaussian quadrature. The product of these areas and the total flux, in k m V 1 for instance, at the given source grid block are the production rates by calving for each size category. The model uses synthetic transient source rates, as output by the continental ice sheet model of Marshall [1996]. So, for the distribution and sizes of Figure C . l and a total flux of, say, Q = 100 km 3 a - 1 the associated probabilities of production and source rates for each size category are as presented in the following table. For the purposes of Index, s L W L (m) Area, As QS ( W a -1 ) 1 0< LWL <100 0.3588 35.88 2 100< L W L <200 0.4722 47.22 3 200< LWL <300 0.1507 15.07 4 300< LWL <400 0.0175 1.75 5 400< LWL <500 0.0008 0.08 Table C . l : Probabilities and fluxes according to iceberg size for a Greenlandic source. Total flux Q — 100km 3a - 1, mean waterline length LWL,GL = 150 m. determining the drift forcings and melt rates, discussed in Appendix D, the midpoint parameter values for each size category are used. To reiterate, using the same procedure with an Icelandic source, where the icebergs produced have a mean waterline length of, say, LWLJS = 90 m the probability density curve is as in Figure C.2. Appendix C. Iceberg Size Distribution at Sources 84 Rayleigh D is t r ibu t ion : mean = 90 Size Parameter: Waterline Length, L^ L (m) Figure C.2: Example: Iceland source iceberg size distribution. If this particular Icelandic source has a flux of Q — 100 km 3 a - 1 , then integrating multiplying the total flux gives the following tabulated source rates. Index, s LWL (m) Area, As Qs (krn^"1) 1 0< L W L <100 0.7090 70.9 2 100< L W L <200 0.2838 28.38 3 200< L W L <300 0.0072 0.72 4 300< LWL <400 1.5xl0"5 0.0 5 400< LWL <500 2 .6xl0 - 9 0.0 Table C.2: Probabilities and fluxes according to iceberg size for an Icelandic source. Total flux Q = 100 km 3 a - 1 , mean waterline length LWL,IS - 90 m. To summarize, each provenance has a unique probability density for iceberg s production. Using a finite number of size categories, specified by the user, the mo< calculates the probable production rate for each size group. Multiplying the total fl Appendix C. Iceberg Size Distribution at Sources 85 by the probability for each size yields the expected production rate for each size within each provenance. As the number of sizes increases the model is potentially more accurate with the caveat that execution time is directly proportional to the number of active sizes and provenances in the run. Appendix D MELTING, SIZE COUPLING AND SEDIMENT DEPOSITION D. l Melting Melting is responsible for ice-rafted sediment deposition and is the mechanism for size coupling. In the model, melting is calculated in accordance with the U.S. Coast Guard -International Ice Patrol (IIP) method to predict iceberg deterioration. In the IIP model, adopted by the continuum iceberg drift model, melting is governed by four mechanism: solar radiation, vertical buoyant convection, air convection and wave action. The fol-lowing presents the IIP empirical equations to predict iceberg melt rates in metres of waterline length melted per day. D . l . l Solar Radiation The daily melt rate A4T of an iceberg is taken as MT = 0.02W ( D . l ) where W is a factor to account for cloud cover, set to 1.0 for cloudy/foggy conditions and 2.0 for clear conditions. Based on qualitative observations of North Atlantic weather conditions, W is always set 1.0 in the IIP model. Since the ocean model input fields used in the continuum iceberg drift model do not include cloud cover the IIP assumption is also used here. 86 Appendix D. Melting, Size Coupling and Sediment Deposition 87 D.l .2 Vertical Buoyant Convection Vertical buoyant convection is the turbulent water-to-ice heat transfer that takes place as fresh melt water from the lower depths of the iceberg rises along its side mixing with the surrounding salt water. This daily melt rate of an iceberg Mc in metres is by this means calculated as Mc = 2.74 x 10-3(2.78T + 0.47T2) (D.2) where T is the difference between the temperatures of the surrounding sea and the iceberg. The temperature of the ice is taken as —1.63° C, the equilibrium temperature of ice in sea water at 30 parts per thousand (ppt). The IIP model uses sea surface temperature as the basis for its equations. The model detailed here uses the vertically-averaged sea temperature of the ocean layers in which the iceberg sits. Thus the sea water temperature experienced by an iceberg in the model depends on its size. Provisions have been made in the continuum model for the inclusion of salinity in the melt equations. However, no expressions have been found in the literature that specifically include salinity. Moreover, based on the findings of Russell-Head [1980] it is believed that, over the range of the concentrations in the model input salinity fields, between 25 and 37 ppt, the effect of salinity on melt rates is marginal. For the sake of thoroughness, though, it is proposed that the IIP buoyant convective term could be modified to account for variable salinity. Buoyant convective mixing at the surface of the melting iceberg is due to differential between fresh and salt and fresh water densities. If, to a first order, the variation in density is taken as a function of salinity, then the turbulent heat could also reasonably be taken a function of sea water salinity. Assuming linear dependence, as an approximation, then D.2 becomes 2.74 x HT3(2.78T + 0.47T2) (D.3) i f - 3 0 30 Appendix D. Melting, Size Coupling and Sediment Deposition 88 where H is salinity in parts per thousand. D.l .3 Air Convection The equation for convective melting in air A4a in meters per day is Ma = 0MFcT (D.4) where Fc is a convective melting factor. The factor Fc depends on waterline length LWL and the relative speed of the iceberg with respect to the local air | XJ^jr | as follows Fc = [0.934 - 0.202log(LWL)] \VAI\ (D.5) when | U . 4 j | < 25 cm/sec. When | U . 4 j | > 25 cm/sec the convective melt factor is Fc = [0.660 - 0A5Uog(LWL)](\VAI\ - 25) + 25 [0.934 - 0.202 log {LWL)] • (D.6) D . l .4 Wave Action Waves deteriorate icebergs by transferring heat to the ice, causing melting, and by in-ducing flexural stresses causing calving. Following the lead of the IIP, only melting is considered in the continuum model. Iceberg melting by wave action KAW is described by the equation r / o\°-2 1 (D.7) Mw = 2 \ 0 ' 2 0.13fe( H Tfw where h is the wave amplitude in centimetres and fw is the wave frequency. Though the model ocean input does not include any information on waves it would be remiss not to attempt to include them, for wave action strongly dominates as a melting mechanism. According to IIP calculations wave action typically accounts for around 80% of daily melt. In the absence of knowledge of wave parameters, representative values for h and fw could be used or for a conservative melt rate they could be set to zero. Appendix D. Melting, Size Coupling and Sediment Deposition 89 D.l .5 Total Melt and Numerical Formulation The total melt rate M, calculated for an iceberg size s in each grid block (i,j), is given by the sum of (D.l), (D.3), (D.4) and (D.7) d dt (LWL) = M = Mr + Mc + Ma + MI (D.S) which represents the daily shortening in metres of the waterline length LWL- Having assumed a standard iceberg shape such that volume is a function of LWL, as discussed in Appendix C, finding the daily volumetric melt rate of ice is trivial. In the model for an iceberg size of s at cell (i,j), the numerical formulation of (D.8) is Ms,i,j = 0.2Witj + Hs,i,j — 30 30 2.74 x 10"3 [2.78(f,,ij + 1.63) + 0.47(2^- + 1.63f 0.01Fci,j(T.,i,3 + 1-63) 0.13/i, h. i,3 0.2 (ftiitj + l.M)fwij (D.9) where T3iitj and HS}ij are the iceberg-size-dependent surrounding sea water temperature and salinity, respectively. This temperature, as with salinity, is calculated as d. E rpn ul iJtlDi. (D.10) subject to Nla E 1=0 where di is the depth of iceberg within ocean layer I, Drs is the size s iceberg draft, Titjii is the local water temperature in layer / and n is the power to which the iceberg melt rate depends. By inspection of (D.9), n is approximately one. For the time steps used in the Appendix D. Melting, Size Coupling and Sediment Deposition 90 model, (D.9) tends to become unstable when included as a term in the finite-difference equation (A.12). Equation (D.9) is therefore solved at each time step by a fourth-order Runge-Kutta routine with no effect on the results. D.2 Size Coupling By having the icebergs melt and thus move to lower size categories, the changing forces acting on the icebergs are realistically modelled. To facilitate this, the model simply keeps account of the fractional change in mass balance for each size of iceberg in each grid block at each time step. In other words, the model tracks the temporal and spatial variation in average iceberg size due to melting. The melting icebergs in each cell are gradually shifted to the next lower size category at each time step. Alternatively, when the cell average iceberg size, expressed as a fraction of its unmelted expected size, decreases below a specified threshold, all the ice in the cell could be shifted to a lower size category An outline of this simple procedure follows. Let Aibfl- and Aib'^ represent the net volume of icebergs of size s ice fluxed into cell over the time interval A i with melting and without melting, respectively. Written in expanded form these are and M ' i j = rtij-1>%? (D.12) The fraction of the ice volume belonging to size 5 icebergs melted over At is then A * * = (D.13) Appendix D. Melting, Size Coupling and Sediment Deposition 91 where Ps,i,j is the total fraction since time zero of the ice belonging to size s icebergs that is unmelted in cell Therefore A P is the fraction by which the average iceberg size has decreased in the cell over the time step. In the five-point finite-difference scheme used in this model the flux into cell is contributed to by the four adjacent cells as well as by its own sources, as outlined in Patankar [1980]. Thus Ps,i,j accounts for fractional sizes, not only of those icebergs already present, but also those from previous locations. In essence, the icebergs in a cell have a memory of where they have been and how much they have melted. At time t the fraction of ice present in the cell, compared to what would exist without melting, is then In general, where A P is zero at time zero, or the ice in the system is initially completely unmelted = ft1 - A P ^ ) (D.15) n=o which will reach a finite value when the iceberg and melt fluxes are at equilibrium such that A P ^ j is zero. The icebergs, or their equivalent ice volumes, are shifted to the next lower size category when for icebergs of size s in cell Ps.ij < ^ (D.16) where vs is the average volume for an iceberg of size category s. As soon as the icebergs in a cell have moved to size category (s — 1) the fraction Ps,i,j goes to one and the process continues. The model simultaneously shifts down all the icebergs in a cell, but could just as easily move portions of that ice for a smoother transition, although this seems to make no significant difference. Appendix D. Melting, Size Coupling and Sediment Deposition 92 D .3 Sediment Deposition Iceberg sediment content in the model is a user-specified input parameter for each size category within each provenance. A high degree of flexibility has been built into this aspect of the model owing to the relative scarcity of data at hand. The model reads input continuous functions for each provenance relating sediment content to iceberg size. As examples, Figure D . l illustrates hypothetical sediment/size curves for Greenland and Iceland. Example: Greenland Icebergs Sediment Content Example: Iceland Icebergs Sediment Iceberg Size Parameter, i.e. Waterline Length (m) teto9 s i i e P a r a m e , e r •'<*• Figure D . l : Example: Sediment distribution by iceberg provenance and size. The curves shown are for illustrative purposes only. To generate the curves shown the model merely requires an appropriate function and the sediment content for the maximum size iceberg. The model assumes that for each ice-berg size the sediment is evenly distributed throughout the iceberg volume. While there is a dearth of information on the actual sediment content of icebergs it is qualitatively clear that it is size dependent. Simply, large icebergs are more likely to have entrained sediment through grounding and contact with their glacier bed. The model allows the Appendix D. Melting, Size Coupling and Sediment Deposition 93 sediment/size function to be readily updated as better information becomes available. If at cell (i,j), for icebergs of size s, the daily volumetric melt rate is a function G of the daily change in waterline length AASiij, described by (D.8), then the daily melt water flux (pa,i,j is <b.iiij = G(M.,i,j)-^(l-8.) (D.17) Pw where pj and pw are ice and water densities and 0S is the the sediment content of a size s iceberg as a fraction of volume. The volume of sediment released daily a is just the product of the melted ice flux and the sediment fraction, as follows a t i i J = G{M.,itj)f3,. (D.18) It is implicit that the above apply to each provenance separately In the final analysis it is of more interest to look at the total melt water and sediment fluxes due to all the iceberg sizes from each provenance. Therefore, the final model output for each provenance is of the form N sizes fiij = ^ (D.19) 1 and Nsizes aitj = a * , * , j (D.20) I Volumetric quantities, such as melt water, deposited sediment and ice volume, are expressed as an equivalent column height at each grid block by dividing by the grid block surface area. Appendix E T I M E S T E P P I N G A N D SPARSE S Y S T E M S O L U T I O N E . l Temporal Discretization The temporal discretization of the model uses a backward difference, or fully implicit, scheme requiring that the ice volume over the problem grid be found at each time step through the solution of a system of linear equations. The fully implicit scheme is first-order accurate in time and implies that the best solution to the system comes through approximating spatial derivatives at the next time step. The size of the time step is determined by the C F L condition for numerical stability, as noted in Appendix A . Simply, the C F L criterion requires that ice does not entirely traverse any grid block in less than a single time step, or A t < min AXi min , min where U and V are x and y iceberg velocity components, respectively. E.2 Sparse Matrix Solution Written as a matrix, the discrete grid block equations of (A.12) results in a large, very sparse matrix, or to be precise, a tridiagonal matrix with fringes. Generically the problem is in the form [A] • [fl = [b] ( E . l ) 94 Appendix E. Time Stepping and Sparse System Solution 95 where at cell Aij = [AN, As, AE, AW, AC], the advection/diffusion coefficients asso-ciated with the Nor th , South, East, West and centre grid blocks in the five-point scheme. The other matrices contain the grid block ice volumes ipij and source/sink terms bij. The system solution is carried out via the method of successive or simultaneous over-relaxation (SOR) with Chebychev acceleration, following the method described in Press [1980]. S O R is an iterative, or relaxation, method that is a variant of the Gauss-Seidel scheme making use of the sparseness of the system of equations by only operating on the non-zero elements. A very general description of the method follows. The matr ix A of ( E . l ) is treated as having two parts; an approximation to A, [E], and a remainder [R] so that [A] = [E] - [R] (E.2) O n this basis ( E . l ) is written as [E] • [fl = [R] • [fl + [b] (E.3) Starting with an ini t ial guess of the unknown grid block parameter value tb^ the iterative scheme successively solves for ib^ in the manner [E] • feW] = [R] • [ ^ } + [b] (E.4) Expanding the remainder term of (E.4) at cell (i,j) at the nth iteration step yields A N 4 : I + + A^:l+A^+AcM (E-5) where looping through i and j the South and West, As and Aw, terms have already been evaluated at step (n + 1). The parameter is then adjusted to get a zero residual Appendix E. Time Stepping and Sparse System Solution 96 which constitutes a Gauss-Seidel iteration. The criterion for stopping is when residual R is below some specified minimum as determined by a suitable matrix norm. Convergence in the Gauss-Seidel method is highly oscillatory. In SOR monotonic convergence is more quickly achieved through the use of a relaxation parameter u>, really a damping coefficient, added to the residual term so that 4T1] = ^ - ^ (E.7) where 1 < co < 2. Again, the determination of OJ and the details of SOR are presented in Press [1980]. Advantages of SOR are that it is easy to program and relatively fast compared to other methods. As for memory, only a single storage vector is required. But unlike direct methods, such as Newton's, that give exact solutions using matrix methods, SOR is not guaranteed to converge. For ill-conditioned problems SOR rapidly becomes unstable and will not converge. In debugging the drift model this characteristic of SOR actually served usefully as a warning flag, for even minor errors in the problem formulation prevented convergence from the outset of a run. Moreover, it is argued that for the sake of efficiency approximate solutions are reasonable for approximate problem formulations. Using the above-discussed solution technique, the ice and meltwater mass balance in the model system well exceeds minimum expectations. Appendix F M A I N P R O G R A M A N D R O U T I N E F L O W C H A R T S 97 endix F. Main Program and Routine Flow Charts Start program Declare constants and allocatable arrays Input iceberg and sediment parameters Allocate arrays Input iceberg source time series Input problem grid and ocean fields Calculate advection-diffusion coefficients Solve cell ice content; SOR Calculate melt rate in each cell Solve melt flux; fourth -order Runge-Kutta Establish problem grid Read flow, no-flow indices Initialize varibles Calculate Coriolis parameters Read ocean depths Read air and water velocity fields Read air and water temperature fields Read salinity fields I Establish iceberg velocity fields Calculated effective temperature fields Calculated effective salinity fields Figure F . l : Main program flow chart. Appendix F. Main Program and Routine Flow Charts 9 9 Start of routine Calculate Coriolis acceleration Calculate iceberg areas above and below water T Calculate forces acting on iceberg and acceleration Solve equations of motion by fourth-order Runge-Kutta Calculate harmonic average iceberg velocities at cell faces End of routine Calculate air drag sumfract = 0, fraction = 0 oldfract = 0,1 = 0 1 = 1 + 1, where I indicates ocean layer fraction = depth(l) / draft sumfract = sumfract + fraction fraction = 1 - oldfract sumfract = 1 Calculate water drag acting on berg by water layer I Calculate acceleration due to Coriolis, water and air drag oldfract = sumfract Return to routine Figure F.2: Iceberg velocity calculation flow chart. Appendix F. Main Program and Routine Flow Charts 100 Start of routine Read multi-layered input field G sumfract = 0, fraction = 0 oldfract = 0,1 = 0 1 = 1 + 1, where I indicates ocean layer fraction = depth(l) / draft sumfract = sumfract + fraction fraction = 1 - oldfract sumfract = 1 Calculate effect of parameter on iceberg by layer I, G(l) Calculate harmonic average parameter at cell faces End of routine Figure F.3: General routine for calculating the effective salinity and temperature acting on an iceberg immersed in multiple ocean layers. 


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