N u m e r i c a l simulation of flow in a laboratory tank using a z-coordinate numerical m o d e l by Sergio Jaramillo Uribe B.Sc. (Physical Oceanography) Escuela Naval Almirante Padilla, 2001 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Oceanography) THE UNIVERSITY OF BRITISH COLUMBIA March 23 2005 Sergio Jaramillo, 2005 Abstract A z-coordinate numerical model is employed to study the dynamic response of an upwelling favorable flow in a rotating laboratory tank with a submarine canyon. A vertically uniform body force is exerted on the fluid inside the tank over a time of 30 s, equivalent to the duration of an upwelling event in the real ocean. The results show a flow pattern that is in qualitative agreement with previous numerical, theoretical and field studies. The agreement between laboratory and numerical results is an improvement on previous results obtained using a terrain-following numerical model. However consistent differences at the upstream rim of the canyon exist. A number of processes are investigated to understand the reason for these differences. Results show that the model is extremely sensitive to changes in the bottom friction and interior mixing coefficients. Although the model can be "tuned" to obtain better agreement, this approach is not used because it could mask other important problems. Instead, results from an increased resolution run using distributed processors show that the observed differences are not related to errors in the bottom boundary layer parameterization. Non-hydrostatic effects are estimated, showing a package of standing internal waves that is generated at the upstream side of the canyon. These internal waves are dominated by a second baroclinic vertical mode that is well predicted by the linear theory. However the non-hydrostatic effects are not strong enough to influence the main characteristics of the flow in the canyon. Different advection schemes for the active tracers are compared, showing that non-linear schemes with a flux limiting algorithm produce the best results from a physical point of view. Flux limiters are designed to regulate the flow in areas with steep gradients of the advected quantity; the influence of these algorithms is investigated and it is shown that they constrain the stretching of isopycnals at the upstream rim and thus they account at least in part for the observed errors in this area. Comparisons between stretching and relative vorticity show that the model results are overwhelmed by frictional effects and that the stretching and compressing of isopycnals is poorly represented. Through these modeling steps, it is concluded that the vorticity observed in the numerical model is more related to flow separation than to stretching vorticity effects. This also contributes to the observed differences between the numerical and laboratory results. ii Contents Abstract ii List of Tables v List of Figures vi Table of symbols x Acknowledgements xiii Chapter 1. Introduction 1 1.1 The coastal ocean and submarine canyons 1 1.2 Background on submarine canyons research 2 1.3 Objectives 11 1.4 Methodology 11 1.5 Structure of the thesis 13 Chapter 2. Massachusets Institute of Technology general circulation model 14 2.1 Model description 15 2.2 Advection schemes in MITgcm 30 2.3 Experiment configuration 39 Chapter 3. Results 53 3.1 Flow description 53 3.2 Comparison with laboratory results 70 3.3 High resolution results 79 iii Contents iv 3.4 Non-hydrostatic run 82 3.5 Comparison of advection schemes 3.4 Vorticity C h a p t e r 4 S u m m a r y and conclusions Bibliography , field 88 98 102 106 List of Tables 1. List of the symbols most used in the present work. x 2. Physical and numerical parameters for the coarse resolution experiment. 47 3. Initial and final depths on 72 floats. 4. Effects of the flux limiting algorithm for five points. v 96 List of Figures 1.2.1 Contoured sections of temperature, attenuation, and stretching vorticity across Astoria Canyon. 5 1.2.2 A schematic of flow over a canyon. 7 1.2.3 Numerical grid of a cross section of the submarine canyon in the laboratory tank. 2.0.1 10 Numerical grid of a cross section of the submarine canyon in the laboratory tank. 14 2.1.1 Sketch of the staggering of the variables using the Arakawa C-grid. 24 2.1.2 Staggering of horizontal grid descriptors. 25 2.1.3 Vertical grid in MITgcm. 26 2.1.4 Sketch of the x-z plane showing the location of the non-dimensional fractions h c and h . 27 2.3.1 Schematic of the bathymetry in the laboratory tank. 40 2.3.2 Top view of the horizontal grid over the canyon. 42 2.3.3 Top view of the bathymetric contours in the model domain. 42 2.3.4 Vertical section of a generic water column in contact with the bottom w topography. 44 2.3.5 Density profile used in the numerical simulation. 45 2.3.6 Potential temperature profile used in the numerical simulation. 46 2.3.7 (a) Buoyancy frequency profile, (b) Dispersion relation for rotating internal waves in the laboratory tank. 2.3.8 49 Shapes of the first three baroclinic modes calculated for the deepest depth in the tank. 50 vi LIST OF FIGURES vii 2.3.9 Top view of the high resolution horizontal grid over the canyon. 51 2.3.10 High resolution numerical grid of a vertical cross section of the submarine canyon in the laboratory tank. 51 3.0.1 Schematic view of the computational domain. 53 3.1.1 Plot of the velocity vectors at level 5 (-1.27 cm). 54 3.1.2 Plot of the velocity vectors at level 6 (—1.55 cm). 55 3.1.3 Plot of the velocity vectors at level 9 (-2.39 cm). 55 3.1.4 Plot of the velocity vectors at level 12 (-3.23 cm). 56 3.1.5 Vertical cross section of the along shore velocity across the shelf (Line C). 58 3.1.6 Vertical cross section of the radial velocity across the shelf (Line C). 59 3.1.7 Vertical cross section of the vertical velocity across the shelf (Line C). 59 3.1.8 Vertical cross section of the density across the shelf (Line C). 60 3.1.9 Vertical cross section of the azimuthal velocity at the canyon mouth (Line A). 61 3.1.10 Vertical cross section of the along shore velocity at the canyon mouth from Allen et al. (2003) using the model SCRUM. 61 3.1.11 Vertical cross section of the radial velocity at the mouth of the canyon (Line A). 62 3.1.12 A cross section of the radial velocity at the mouth of the canyon from Allen et al. (2003) using the S C R U M model. 3.1.13 62 Vertical cross section of the vertical velocity at the mouth of the canyon (Line A). 3.1.14 63 Vertical cross section of the vertical velocity at at the mouth of the canyon from Allen et al. (2003). _ 3.1.15 Vertical cross section of the density at the mouth of the canyon (Line A). 3.1.16 A cross section of the density at the mouth of the canyon (Line A) from Allen et al. (2003). 3.1.17 63 64 65 Vertical cross section of the azimuthal velocity at approximately the middle of the canyon (Line B). 65 LIST OF FIGURES 3.1.18 Vertical cross section of the radial velocity at approximately the middle of the canyon (Line B). 3.1.19 66 Vertical cross section of the vertical velocity at approximately the middle of the canyon (Line B). 3.1.20 viii 66 Vertical cross section of the density at approximately the middle of the canyon (Line B). 67 3.1.21 Vertical cross section of the azimuthal velocity along the canyon axis (Line D). 67 3.1.22 Vertical cross section of the radial velocity along the canyon axis (Line D). 68 3.1.23 Vertical cross section of the vertical velocity along the canyon axis (Line D) 69 3.1.24 ertical cross section of the density along the canyon axis (Line D). 69 3.2.1 Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm (filled blue squares), and S C R U M (unfilled black squares). 3.2.2 71 Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm reduced drag run (filled blue squares), and MITgcm full drag run (unfilled black squares). 3.2.3 75 (a) Vertical cross section of the radial velocity at the mouth of the canyon (Line A) with the drag reduced, (b) radial velocity with drag reduced minus original radial velocity. The view is toward the head of the canyon. 3.2.4 Vertical cross section of the density at the mouth of the canyon (Line A) for the reduced drag run. 3.2.5 76 77 Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm free slip run (filled blue squares), and MITgcm full drag run (unfilled black squares). 3.3.1 Results from the high resolution run (480x164x64 points) upstream from the canyon along Line C. 3.3.2 80 Results from the high resolution run (480x164x64 points) at the mouth of the canyon (Line A). 3.3.3 78 81 Results from the high resolution run (480x164x64 points) at the middle of the canyon (Line B). 82 LIST OF FIGURES 3.4.1 Differences between hydrostatic and non-hydrostatic runs at the mouth of the canyon (Line A). 3.4.2 ix 84 Top view of the differences in the vertical velocity between hydrostatic and non-hydrostatic runs at 1.2 cm depth, (a) t=10 s. (b) t—20 s. (c) t=30 s. (d) t=35 s. 3.4.3 86 Comparison between calculated baroclinic vertical mode shapes and vertical velocity profiles from the numerical model. 3.5.1 Temperature in the surface layer of the tank approximately across the head of the canyon (inshore from Line B) for different advection schemes. 3.5.2 87 89 Temperature anomaly at 1.2 cm depth (this is the depth of the maximum temperature anomaly on the shelf) calculated using different advection schemes. 92 3.5.3 Temperature along Line A above the rim depth (1.8 cm) for the 3-DST advection schemes with and without using a flux limiting algorithm for the low and high resolution runs. 3.5.4 94 Contour plot of the Sweby limiter for the vertical flux at 2.04 cm depth. (a)V»(r ). (b)V(r ). - + 98 3.6.1 Vorticity across canyon mouth. 100 3.6.2 Velocity vectors at 2.1 cm depth, a) t=10 s. b) t=20 s. c) t=30 s. d) t=35 s. 101 Table 1: List of the symbols most used in the present work. Symbol Explanation _ i f — (ft+l/2+Pi-l/ ) 2 2 A Horizontal eddy viscosity A Vertical eddy viscosity h v a Radius of the earth at Thermal expansion coefficient b a, b Surface and bottom control parameters in S C R U M Buoyancy Ps Salinity contraction coefficient B Burger number u Si 5iP = Pi+i/2 - Pi-1/2 5 X Parameter used to stabilize the Adams-Bashforth method Parameter controling the presence of non-hydrostatic terms in the equations V Free surface F Flux F Froude number r Forcing terms in the governing equations 9 Ge Gravitational acceleration Geometric parameter Summarize all the terms in the momentum equations except the temporal (*u,v,w,S,9 and surface pressure gradient terms Summarize all the terms in the momentum equations except the temporal, surface, and hydrostatic pressure gradient terms h e Parameter related to the area where more resolution is required in SCRUM X xi Table 1: List of the symbols most used in the present work. Symbol H s i, 3, k Explanation Depth of the shelf break Spatial index relative to x, y, z directions. Horizontal Laplacian diffusion coefficient K h Vertical Laplacian diffusion coefficient K z L Length of the canyon N Buoyancy frequency V Horizontal gradient operator :Vp = {S p, S p} x y Horizontal divergence operator : V- v 2 V • 17 = ^- {c\Ayu + SjAxv} Laplacian operator : S/ p = V • V p 2 Earth's rotation rate V Phyd Hydrostatic Pressure <t> Latitude IP Flux limiting function R Radius of curvature at the shelf break isobath n Linear bottom drag coefficient R Rossby number +- Slope ratio of an advected tracer 0 r Pressure P Density s Vertical coordinate used in S C R U M S Salinity a Linear function of z used in terrain-following models as the vertical coordinate. Table 1: List of the symbols most used in the present work. Symbol Explanation t Time T Stress e Potential temperature u Zonal velocity (azimuthal in the tank) U Velocity at the shelf break depth U, V, w, T Generic variable (u, v, w, S, 9) V Meridional velocity (radial in the tank) V, V, V , u Fluxes in the x, y, z directions v w V e Volume of the grid box surrounding the u, v, w, 9 points w Vertical velocity W Width of the canyon mid-way along its length w sfc Width of the canyon at shelf break X Zonal coordinate (azimuthal in the tank) y Meridional coordinate (radial in the tank) z Vertical coordinate Acknowledgements I want to thank my supervisor Dr. Susan E. Allen for her constant advice throughout these years at U B C , and especially for the opportunity that she has given me by involving me in the field of numerical modeling of submarine canyons. I also wish to thank Dr. William Hsieh for his agreement to examine my thesis in such short notice. For the same reason I thank Dr. Richard Pawlowicz, but also for allow me to stretch my sea legs by letting me participate in his cruises whenever I could. My gratitude is extensive to all the other students and post-docs in the Waterhole. A l l of them have helped me in some way along my masters, either with technical advice or by sharing a pint of beer. Especially Amy Waterhouse who, with the best disposition, helped me in issues ranging from laboratory canyons to shopping for groceries. I am grateful to the members of the Geophysical Disaster Computational Fluid Dynamics Centre at the University of British Columbia, for helping me with the use of their Monster Super-Cluster, especially with all the compilation problems that I had. Thanks to my mother,and sister who regardless of their own personal problems always worried about my well being. I particularly want to thank my father, because he taught me about the law of gravity. Finally, I want to say thanks to my wife Paola, who left her family, a good job, and the warm weather of Colombia to accompany me to Canada. She has always supported and encouraged me with both reason and love, especially in the darkest hours of these years. xiii CHAPTER 1 Introduction 1.1. The coastal ocean and submarine canyons The coastal ocean is one of the most interesting places for oceanography. Some of the most biologically productive areas can be found here, and it is home to intriguing physical processes that include coastally trapped waves, coastal upwelling, internal wave generation and dissipation, tides, and boundary and buoyancy driven currents. The main forcing mechanisms for coastal circulation are: winds, buoyancy due to fresh water input, and tides. When the wind blows parallel to the coast it generates an Ekman drift forming an angle of 90 degrees to the right (left) of the wind in the Northern (Southern) hemisphere, i.e. offshore or onshore. If the coast is on the left (right) of the wind in the Northern (Southern) hemisphere, the drift is offshore. This drift will displace the water in the upper layers creating a lower pressure center near the coast, which generates a geostrophic current along the coast. Water from deeper layers will have to move upward to replenish the water displaced by the offshore drift in a phenomenon known as coastal upwelling. This situation has an important impact on the productivity of coastal areas as it may bring nutrients into the euphotic zone making them available for biological uptake. Continental shelves bound practically every continent and many islands of the world, and represent a transition area from shore processes to deep oceanic processes. They vary in width from a few kilometers to more than a thousand kilometers in some parts of the world, with about 100 km in average. The edge of the shelfbreak occurs at a depth that ranges from 20 to 550 m, with an average of about 200 meters. The shelfbreak orientation typically mimics the coastal orientation but in many parts of the world (e.g. North - West coast of North America and the coast of the Mediterranean); the continental shelf is indented by submarine canyons. Hickey (1995) estimates that canyons occupy 20% of the shelf edge from Alaska to the Equator, and on the Vancouver Island - Oregon coast they cover almost 50% of the shelf edge. Observational (Freeland and Denman, 1982; Hickey, 1997; Allen et al., 2001), theoretical (Chen and Allen, 1996; Allen 2000) and numerical studies (Allen, 1996; Klinck,1996) suggest that in the 1 1.2. B A C K G R O U N D ON S U B M A R I N E C A N Y O N R E S E A R C H 2 vicinity of submarine canyons coastal upwelling is enhanced. Upwelling favorable shelf currents, (i.e. with the coast on their left in the Northern Hemisphere) produce unbalanced up-canyon pressure gradients in narrow canyons that can drive up-canyon flow enhancing the upwelling and thus the cross-shelf exchange of water from the open ocean to the coastal ocean. Kunze et al. (2002) suggest that, when canyons reach shallow water, enhanced turbulent mixing within the canyon may also supply nutrients into the euphotic zone, contributing even more to local productivity. For these reasons, if we want to understand the coastal ocean as a system, the study of the dynamics of submarine canyons is required. 1.2. Background on submarine canyon research Previous results from investigations about submarine canyons that have motivated the present study are summarized in this section. 1.2.1. T h e o r e t i c a l a n d field s t u d i e s . Oceanographic measurements are very difficult to ob- tain near submarine canyons due to the steep slopes and heavy fishing that occur near these areas. Hickey (1995) comments that historically most of the sampling in submarine canyons was conducted by geological oceanographers who were interested in turbidity currents; therefore most of the data collected was concentrated along the canyon axis (e.g. Shepard et al., 1974), and, unfortunately, this is not enough to form an over all picture of the variability of the flow. For these reasons there is only a handful of studies that have used field data to investigate the dynamics of submarine canyons (e.g. Preeland and Denman, 1982; Kinsella et al., 1987; Cherubin et al., 1990; Hickey, 1997; Allen et al. 2001; Cherubin et al., 2003). Even though the conclusions that are drawn from these studies are still the subject of discussion, most of them agree in suggesting that the basic mechanism exposed by Freeland and Denman (1982) is responsible for the enhanced upwelling observed near submarine canyons. This mechanism, based on observations obtained in their study near a small branch of the Juan de Fuca Canyon, can be summarized as follows: far from the canyon the flow is assumed to be in geostrophic balance, so that the Coriolis force that is acting on the along-slope flow is in equilibrium with the cross-shelf pressure gradient (which is directed in-shore under upwelling favorable conditions), over the canyon this situation is the same, but within the canyon the along-slope flow is reduced and thus the Coriolis force cannot balance the pressure gradient. This results in an acceleration of the in-shore flow inside the canyon, and therefore the upwelling is increased. 1.2. B A C K G R O U N D ON S U B M A R I N E C A N Y O N R E S E A R C H 3 However, the argument of Freeland and Denman (1982) neglected to consider the implications that the shape of the canyon may have on the behavior of the flow. Other studies expanded the theoretical understanding of the dynamics in submarine canyons to include such geometrical considerations, and have shown that the width of the canyon determines the strength of the cross canyon flow, and therefore, the strength of the effect of the canyon on the overlying coastal current, with the influence of the canyon becoming smaller as the canyon width decreases (Klinck, 1989). The length of the canyon is also of importance. For instance, for canyons such as Juan de Fuca Canyon (which are considered long, as it practically reaches the coast), the geostrophic solutions show that the flow is generated by first order geostrophic dynamics, while in canyons such as Astoria Canyon (considered short as it reaches approximately to the middle of the continental shelf), the in-canyon flow is generated by second-order dynamics, such as time dependence and non-linearities. This leads to strong and deep seasonal upwelling in Juan de Fuca Canyon, and weaker episodic upwelling in Astoria Canyon (Allen, 2000). One of the most distinctive features of the flow in a submarine canyon is the presence of a cyclonic eddy in the interior of the canyon. This eddy spans the whole depth of the canyon, and in some occasions is observed above the canyon rim. For example, Hickey (1997) presents observations of an eddy rising above the canyon rim in Astoria Canyon, and Allen et al. (2001) observed a closed cyclonic eddy above Barkley Canyon off the west coast of Vancouver Island. Different mechanisms have been suggested to be responsible for the vorticity generation in the canyons. First, one can consider that below the rim the vorticity pattern can be related to boundary separation of the slope flow that, when it reaches the strong turn at the upstream side of the canyon, is not able to follow the bathymetry (Hickey, 1997). Another mechanism that can explain the vorticity inside the canyon is related to the stretching of the water columns. Potential vorticity is a dynamic quantity that is conserved for inviscid, non-diffusive flows. This conservation principle requires that when the fluid layers are stretched or compressed, the relative vorticity is increased or decreased, respectively. Inside the canyon under upwelling conditions, there will be stretching of the water columns (and therefore increase in the cyclonic vorticity) if the upwelling is stronger for the layers near the rim and weaker for the deeper layers. Above the rim, if shelf water falls into the canyon, the water columns will be stretched leading to an increase in cyclonic vorticity on the upstream side of the canyon. Hickey (1997) calculates the stretching vorticity assuming conservation of potential vorticity and 1.2. B A C K G R O U N D O N SUBMARINE C A N Y O N R E S E A R C H 4 negligible upstream relative vorticity, as: V=(^-ljf (1.2.1) S where h is the initial thickness of a layer of water, h is the stretched thickness, / is the Coriolis c parameter, and SV is the stretching vorticity. In this way she can explain the order of magnitude of the vorticity observed in Astoria Canyon (Figure 1.2.1). Using similar arguments, Allen et al. (2001) show that the stretching vorticity generated over Barkley Canyon is strong enough to produce a closed cyclonic eddy over the rim of the canyon. Scaling analysis has been one of the most important tools for the understanding of the canyon dynamics. Allen and Hickey (in preparation) perform a detailed scaling analysis aimed to identify the non-dimensional numbers that are essential to represent the dynamics of the flow in a submarine canyon. The scaling is based on observations from Astoria Canyon, Barkley Canyon, Carson Canyon, Quinault Canyon, and a canyon off Tidra in North Africa. The flow considered for the scaling is nearsteady advection-driven upwelling of an inviscid stratified fluid. The canyon under consideration is a short canyon. They argue that friction can be neglected from the scaling, because the vorticity generated by stretching of the water columns is enough to account for the observed vorticity in canyons such as Barkley Canyon (Allen et al., 2001). These parameters are: • DH = fL/N, where / is the Coriolis parameter, L is the length of the canyon, and N is the buoyancy frequency at the shelf break depth, represents the vertical scale for the processes inside the canyon. • a Rossby number, R A = U / / i ? , where U is the incident velocity at the shelf-break depth, and R is the radius of curvature of the shelf break isobath at the upstream corner of the mouth of the canyon. This number determines the degree to which rotation is important in a flow, specifically for canyons, it represents the capacity of the flow to cross the isobaths and flow across the canyon. e a Rossby number, RL = U/fL. This number measures the ability of the flow to lift isopycnals (compared to DH). • a Burger number, B U = NH /fV\l s : where H S is the depth of the shelf break, and W is the width of the canyon at mid-way along its length. This numbers determines, generally, whether vertical motions are governed by stratification or rotation or both. For submarine 1.2. BACKGROUND ON SUBMARINE CANYON RESEARCH F I G U R E 1.2.1. Contoured sections of temperature, attenuation, and stretching vorticity across Astoria Canyon on 21 May/1983. The cross section is near the mouth of the canyon, and the view is toward the canyon head. Note the layer of positive vorticity at the upstream corner and the interior of the canyon that results from the stretching of water columns. Reproduced from Hickey (1997) with permission of Dr. B. M. Hickey. 5 1.2. B A C K G R O U N D ON S U B M A R I N E C A N Y O N R E S E A R C H 6 canyons the inverse of the Burger number, and the ratio of the canyon length to the radius of curvature, determines if the eddy at rim depth is present. • a Froude number, Fr = \J/NH . S This number measures the importance of stratification. In general if Fr < 1 the stratification effects are important; as Fr becomes smaller, the importance of stratification increases. • a geometric parameter, Ge = V\l/\N b, where s is the width of the canyon at shelf break depth. This parameter describes the horizontal shape of the canyon. The less tapered a canyon is, the more likely it is to develop an eddy at the rim depth. From these studies the conceptual picture of the flow in and around a short submarine canyon (taking Astoria Canyon as an example) during an upwelling event (southward flow in the northern hemisphere) can be divided in four parts (see schematic depicted in Figure 1.2.2, from Allen et al., 2001): (1) Near the surface, the flow is basically unaffected by the presence of the canyon. (2) Just above the canyon rim, the flow is strongly affected by non-linear advective effects. The flow falls from the shelf into the canyon causing stretching of the water column and subsequent cyclonic vorticity. This vorticity can be strong enough to generate a closed eddy. The flow is deflected up canyon crossing the downstream rim near the head of the canyon. After it crosses the downstream wall, the water columns are compressed generating anticyclonic vorticity and the flow turns offshore. (3) The flow below the shelf break enters the canyon following the bathymetry, but once inside the canyon, it is constricted by the canyon walls. The Coriolis force is reduced and the geostrophic balance breaks producing an unbalanced pressure gradient that drives the flow up the canyon and forces it to emerge over the downstream side of the canyon near the head. (4) Deep inside the canyon, the up-canyon flow is reduced, but as in the middle layers the upwelling continues, the water column is stretched generating a cyclonic eddy trapped within the canyon walls. Separation of the slope flow may also contribute to increase the cyclonic vorticity inside the canyon. 1.2.2. Laboratory and numerical studies. Numerical modeling of the flow inside and over submarine canyons has proven to be very challenging for several reasons such as the presence of open 1.2. B A C K G R O U N D ON S U B M A R I N E C A N Y O N R E S E A R C H 7 F I G U R E 1.2.2. A schematic of flow over a canyon. The near surface and deep layer are shown explicitly. The near surface flow-passes over the canyon with little deflection. Near the rim the flow "falls" into the canyon generating cyclonic vorticity. Below the rim the flow is constricted by the bathymetry and turns inshore and is upwelled onto the shelf. Flow in the deep layer turns cyclonically. Reproduced from Allen et al. (2001) with permission of Dr S. E. Allen. boundaries, steep slopes, and sudden changes in isobath orientation. In one of the first attempts to model submarine canyons, Allen (1996) used a three-layer model to study the spin-up (one-half to two day timescale), steady state (longer timescales), linear and non-linear dynamics of a rectangular canyon with vertical walls. Klinck (1996) used a semi-spectral primitive equation model (SPEM, Haidvogel et al., 1991) with the equivalent of 13 levels in the vertical, to investigate the steady state response of a rectangular canyon to upwelling and downwelling favorable incident flows using both strong and weak stratification. Both modelling studies showed good qualitative agreement with observations over Astoria Canyon (Hickey, 1997), especially in the vorticity patterns near the rim during upwelling, the cyclonic pattern inside the canyon, and the flow near the surface where the influence of the canyon is not observed. However, the eddy above the rim was not predicted. Due to the scarcity of comprehensive field data, the use of laboratory scale models of submarine canyons is regarded as an alternate approach to understanding the dynamics of submarine canyons and a useful tool to determine differences between laboratory and numerical modelling results. Following this philosophy, Perenne et al. (2001a, b) obtained good results for an idealized short canyon (less steep than most of the real canyons) using the spectral element ocean model (SEOM) when comparing it with a laboratory scale model on a rotating table. However, when analyzing the residual currents in the model, they found that the laboratory currents near the rim of the canyon easily separate from the shelf bottom and flow over the deeper fluid layers of the canyon, while 1.2. BACKGROUND ON SUBMARINE CANYON RESEARCH 8 the numerical velocities exhibit a greater tendency to follow the bottom slope. They attribute this misrepresentation of the vertical separation by the numerical model to the parameterization of the bottom boundary layer (which was not resolved due to the vertical resolution used) or the fact that the model is hydrostatic, and non-hydrostatic effects could be important in this area. Allen et al. (2003) use a terrain following coordinate numerical model to compare to a laboratory model of upwelling over a short submarine canyon, also with good results. However, they found differences in the flow pattern at the canyon rim both in the vertical and offshore components of the velocity. As an important part of the present work is devoted to the comparison with the results obtained by Allen et al. (2003), next I briefly discuss the classification of numerical models according to the vertical coordinate that they use, and I offer a short description of the numerical model used in that study, including the errors that they found. The numerical models of oceanic circulation are based on approximations of the Navier Stokes equations. These equations are normally written in geopotential coordinates (z or height), but can can also be written with different choices for the vertical coordinate. The most widely used alternatives are: • Terrain-following (a) coordinates: One of the main challenges of numerical modeling of the ocean is the accurate representation of the topography. With ^-coordinate models, an accurate representation of the bottom topography requires a very dense grid that can become computationally expensive, and the phenomena that occur near the bottom can be poorly represented (Haidvogel and Beckmann, 1999). In order to avoid these problems, some numerical models use the cr-coordinate system. The a-coordinate is a linear function of z which allows the mapping of the topography into a regular domain and simplifies the numerical computations. The lowermost coordinate level is assigned to the bottom of the ocean, and thus this coordinate system is useful to investigate processes with strong dependence on the bottom topography, like deep water circulation and turbidity currents. • Density (isopycnal) coordinates: Density models, also known as layered models, divide the water column in constant density layers. The physical reason for this division is to take advantage of the fact that most oceanic currents flow along isopycnals surfaces. Mathematically this is equivalent to mapping from the vertical coordinate z to a density coordinate system. In this coordinate system, moving constant-density layers are used, and each layer thickness is treated as a prognostic quantity. The advantages of this system include: ease of 1.2. B A C K G R O U N D ON S U B M A R I N E C A N Y O N R E S E A R C H 9 development, minimization of cross-isopycnal diffusion, elimination of pressure gradient errors, representation of baroclinic processes, and cost savings over a fully three-dimensional formulation (Douglas et al., 2003). This kind of models are useful to investigate phenomena such as thermohaline fronts, eddy generation, and wind-driven circulation. Because the discretization of these models introduces truncation errors that are different for each vertical coordinate system, the decision in favor of the use of one model depends on the physical problem to be studied, and the impact that the numerical errors particular to the vertical coordinate system used will have on the solution of the problem. 1.2.2.1. S-coordinate Rutgers University Model (SCRUM). The linear mapping of the topography done by a- coordinate models permits one to resolve the bottom topography very well, but does not allow the same resolution in the surface layers. In the development of S C R U M , Song and Haidvogel (1994) use a nonlinear function of z that it makes it possible to have high resolution in the upper ocean and at the same time retains the terrain-following advantages of the a-coordinate system. This s-coordinate in S C R U M is: (1.2.2) z = j] (1 + s) + h s + (h - h ) C (s), - 1 < s < 0, e e where, r ( B s C { s ) /, = { l - ^ sinh (as) b ) -^T + b sinh [a (s + 1/2)] - tanh ((1/2) a) 2 tanh ((1/2) a) where a (0 < a < 20), and b (0 < b < 1) are surface and bottom control parameters, respectively, that regulate the density distribution of layers in this way: the greater a the higher the number of layers near the surface, and the greater b the higher the number of layers near the bottom. The depth is h, the constant h is related to the area where more resolution is required (minimum depth, e surface or bottom boundary layer), r](x,y,t) is the free surface, where x and y are the horizontal coordinates, and t is the time. A cross section of the numerical grid in a laboratory tank with a submarine canyon illustrates the treatment of the vertical coordinate by S C R U M (Figure 1.2.3). Allen et al. (2003), reported truncation errors in S C R U M due to the abrupt change in the curvature of the bottom topography when comparing the results of the numerical model with laboratory experiments. The laboratory experiment considered a tank with a submarine canyon, forced 1.2. B A C K G R O U N D O N S U B M A R I N E C A N Y O N R E S E A R C H 10 Zonal Cross Section (Canyon) 0-01.. I I S.OLJ 30 I 1 L H-M-4.' I I I ' I ' I ' 4 - M - H I 1 . -.• n 23 f -!- I , V': 40 • • I I <5 I ' I 1 I J 1 I ' I R • • I 1 50 F I G U R E 1.2.3. Numerical grid of a cross section of the submarine canyon in the laboratory tank. The vertical scale is in cm. The horizontal scale is in grid points with horizontal grid spacing dx « 6.2 cm for this cross section. Reproduced from Allen et al. (2003) with permission of Dr. S.E. Allen. with an impulse-like body force for the equivalent duration of an upwelling event in the coastal ocean (the complete apparatus will be discussed in the next Chapter). In their results, the numerical vertical velocity showed small scale spatial features. Upward velocity on the upstream side of the canyon, and downward velocity on the downstream side were predicted, both of which are contrary to the current understanding of the dynamics of submarine canyons and to available observations of the isopycnal distribution over Astoria Canyon (Hickey, 1997). Comparison of particles tracks between the laboratory and the numerical model showed also that the model was consistently underestimating the off-shore component of the velocity field at the upstream side of the canyon. In S C R U M , the most important term for the calculation of the rate of change of the density is the vertical advection. For the numerical simulation only the temperature determines the density, and in this model, the advection of temperature is calculated using the Shchepetkin and McWilliams (1998) Gamma advection scheme. Using this scheme, Allen et al. (2003) calculates a positive 1.4. M E T H O D O L O G Y 11 truncation error in the vertical advection at the upstream r i m of the canyon of about 70%. This error is large enough and with the correct sign to account for the observed differences. They found that this error is due to the use of terrain following coordinates combined with a strong vertical gradient i n density, vertical shear i n the horizontal velocity, and topography with strong curvature. W i t h the exception of the problem related to the choice of a terrain-following model, all of these conditions are going to be present whenever the flow near a submarine canyon is modeled. For this reason, a z-coordinate model should be the next logical tool to be used i n the investigation of the dynamics of submarine canyons. 1.3. Objectives From the literature review presented there are still many questions about the dynamics of submarine canyons that remain unanswered. However, only some of those questions fit into the scope of this thesis, and they represent the objectives conceived for this research: (1) Would the use of a z-coordinate model (as opposed to a terrain-following model) reduce the magnitude of the errors in the vertical advection of density, and hence the discrepancies between numerical models and laboratory observations at the upstream r i m of the canyon? (2) What is the importance of friction for the circulation i n and around the canyon? (3) How important are the non-hydrostatic effects for the flow i n submarine canyons? (4) Considering a uniform incident flow, what would be the depth of the deepest upwelling in a submarine canyon? (5) W h y have numerical models not been able to predict the off-shore velocities at the upstream wall of the canyon observed in laboratory models? (6) W h a t is the contribution of the stretching vorticity to the relative vorticity in submarine canyons? 1.4. Methodology Different approaches can be taken to investigate the questions presented i n the previous section. Numerical models have been used widely to investigate the dynamics of the coastal ocean circulation, however the results from these models have to be validated first. This validation can be done through comparison with analytical solutions, but for the flow in a submarine canyon there are not analytical solutions available. Another approach would be to compare numerical solutions with field 1.4. M E T H O D O L O G Y 12 data, however in this kind of comparison it is very difficult to isolate possible sources of errors in the model, because there are many physical processes that have to be accounted for (i.e. turbulence, tides, fresh water input). The comparison of results between two numerical models could be used too, but in the case of a submarine canyon, there are already many questions regarding the performance of models used in previous studies to consider them a trustworthy frame for comparison. For these reasons, the alternative of using laboratory observations to verify the results of a numerical model was the methodology chosen for the present work. With the appropriate laboratory techniques and scaling analysis, laboratory models provide a good test case for numerical models. As there are no truncation errors, turbulence, or sub grid scale parameterizations to deal with, it is easier to isolate possible sources of errors for simplified geophysical problems under laboratory conditions. The flow in the laboratory, however, has not exactly the same characteristics as the flow in the real ocean because of scaling limitations (e.g. the necessity to exaggerate the vertical scale to reduce viscous effects). However, regardless of these limitations, numerical models should be able to simulate the laboratory experiments in a satisfactory way, and the discrepancies between the laboratory and the numerical results can help to understand any model shortcoming before the numerical model is applied to more realistic oceanic conditions. For the present work, the experiment that is going to be analyzed consists of an impulsively started flow in a rotating tank with a submarine canyon. The rotation rates, direction, and duration of the forcing are designed to simulate an upwelling event in the real ocean. The laboratory conditions are the same as described in Allen et al. (2003). The numerical modeling is done with the "Massachusetts Institute of Technology general circulation model (MITgcm)", which is configured to reproduce the conditions of the laboratory experiment. This model was chosen mainly because it meets the following requirements: • Uses a geopotential coordinate system (z-coordinate). This avoids the truncation errors inherent to terrain-following models, and therefore should give better results at the upstream wall of the canyon. • Is able to solve the non-hydrostatic Navier-Stokes equations. This is necessary to investigate objective 4, the role of non-hydrostatic effects. • Is written in MPI which allows it to run on a Beowulf cluster to give resolution high enough to resolve the bottom boundary layer. 1.5. S T R U C T U R E OF T H E THESIS 13 Some modifications were made to customize MITgcm for the comparison with laboratory data. These modifications are described in detail in Chapter 2, and include changes to the representation of the bottom topography in the calculation of vertical dissipation terms, and the inclusion of 3-dimensional advection of floats for runs with single and multiple processors. 1.5. Structure of the thesis The rest of the thesis is organized as follows: In Chapter 2, there is a description of the numerical model used in this work, of the laboratory and numerical model configurations, and of the parameters choices made for this study. Chapter 3 describes the results of the low resolution run for different areas of the domain and compares them against the laboratory observations available. It contains a description and discussion of the importance of the frictional effects in the numerical results by testing the sensitivity of the model to different values of the bottom drag, a description of the results from a high resolution run aimed to identify any errors related to the bottom boundary parameterization. Then, the results from a non-hydrostatic run are described and the presence of internal waves in the model is analyzed. A comparison between the different advection schemes available in the model is presented using the calculation of the depth of the deepest upwelling to identify the effects of the choice of scheme on the models results. The characteristics and performance of flux limiting algorithms are discussed with emphasis on the repercussions on the model results. The stretching and relative vorticity fields are compared for the low and high resolution runs. In Chapter 4 some conclusions and directions for future work are presented. CHAPTER 2 Massachusetts Institute of Technology general circulation model (MITgcm) The very large truncation errors described in section 1.2 are a direct consequence of the use of terrain-following coordinates when the topography changes abruptly. For this reason, in the subsequent study of the flow in a submarine canyon, a model that uses z-coordinates was selected instead. MITgcm employs what is called the intersecting boundary method, which allows the boundary to intersect a grid of cells modifying the shape of the intersected cells, thus creating a partial step representation of the bottom topography (Adcroft et al., 1997, Figure 2.0.1). -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 -0.045 0.16 0.18 0.2 0.22 0.24 0.26 FIGURE 2.0.1. Numerical grid of a cross section of the submarine canyon in the laboratory tank. The bold line represents the real bathymetry. Scales in meters. The representation of the bottom topography in this model can be in the form of full steps or partial steps. The former is a very crude representation of the topography that needs very high resolution to reasonably represent any problem with strong dependence on the bottom, as was mentioned in section 1.2. The second approach is more accurate and allows the bottom boundary to intersect the grid cells modifying the shape of the cells that are intersected. This partial step 14 2.1. MODEL DESCRIPTION 15 representation of the bottom topography produces a significant improvement in the results when compared to the full step approach, permitting at the same time a better approximation to the real topography with less vertical resolution (Adcroft et al., 1997). MITgcm also has the possibility of hydrostatic primitive equations (HPE), quasi-hydrostatic (QH), and non-hydrostatic (NH) treatment of the equations for oceanic and atmospheric circulation and allows the use of different curvilinear grids. It has many different options for physical parameterizations and can run on personal computers or parallel computers which reduces the constraints related to the resolution in the model. This section will describe the main aspects of the model, including the system of equations used, boundary conditions, method of solution, discretization, and parameterizations relevant to the present study. Making just a few adjustments, the model can be applied to both the ocean or the atmosphere, however the terminology that is going to be used is appropriate for the ocean. For a more detailed explanation see Marshall et al. (1997a, b), or the MITgcm web site (http://mitgcm.org/) . 2.1. Model Description 2.1.1. System of equations. Most of the numerical models (including SCRUM) use the approximate form of the Navier Stokes equations known as the "Hydrostatic Primitive Equations", where the vertical momentum equation only considers the hydrostatic balance. On large scales this approximation is usually very good, but when dealing with small scales, or when the possibility of having important non-hydrostatic processes, such as internal wave propagation, is open, then the assumptions made for the HPE's can produce difficulties. MITgcm offers the alternative of using a quasi-hydrostatic approximation, where the balance between gravity and pressure gradient is relaxed, and a non-hydrostatic approximation, that uses the incompressible Navier Stokes equations. These equations can be treated in polar spherical coordinates, Cartesian coordinates or general curvilinear coordinates. The state of a Boussinesq fluid at any time can be characterized by the distribution of the velocity (v ), potential temperature (6), salinity (5), pressure (p), and density (p) fields. Defining a kinematic pressure as where Sp is the deviation of the pressure from that of a resting hydrostatically balanced ocean and p f is a constant uniform reference density. The equations that control these fields in MITgcm are: re 2.1. M O D E L DESCRIPTION 16 • Zonal momentum: (2.1.1) ^ - {-2fiu sin 0 + 2fhu cos 0} - if — = -1? • V u - ^ dt (_a_ a J ' L + F u ox where ~v = (Vh,U>) = (u,v,w) is the velocity vector, u is the zonal velocity, v is the meridional velocity, w is the vertical velocity, <j> is the latitude, a is the radius of the Earth, fi is the Earth's rotation rate, and g is the gravitational acceleration. • Meridional momentum: (2.1.2) — = -1? -Vv- at < 1 } - 2fiu sin c/> + F J ay v a a Vertical momentum: (2 13) 1 ^ = J L ; where 5p — p - - ^ . V . = p ef r ( ^ j + I. 2 O n c o s 0 - ^ - | + a I Pre} + ^ OZ is perturbation density. a • Continuity: (2.1.4) V • 17 = 0 • Heat: an (2.1.5) | ! =-V • (1ft) 4- • Salt: (2-1.6) — = - V • (ITS) + F S • Equation of state: (2.1.7) p= p(6,S,z) On the right hand side of (2.1.1) and (2.1.2) the terms are the advection, earth curvature, Coriolis, pressure gradient, and forcing-dissipation terms, respectively. The terms are the same in (2.1.3), except the fourth term which is the buoyancy term. 2.1. M O D E L DESCRIPTION 17 The non-hydrostatic approximation uses all the terms in these equations, the quasi-hydrostatic ignores the terms underlined twice, and the hydrostatic approximation ignores all the underlined terms. 2.1.2. Boundary conditions. In the vertical, the boundary conditions are w = 0 at the bottom, and Dz f w = ——^ sur L Dt at the surface, where DZsurf dz f Dt is the total derivative, z f sur _ sur dt s = z + ?], z is the position of the resting surface of the ocean (zero 0 0 meters for the ocean), and ry is the deviation in meters of the surface of the ocean. In the horizontal, the boundary conditions are that of no flow through any solid boundary, ~v • ~n = 0, where ~n is the normal to the solid boundary, and the no-slip (VT = 0), or free slip (dvr/dn — 0) conditions can be used for the tangential velocity VT on all solid boundaries. 2.1.3. Forcing. The different forcing terms {F g s) on the right hand side of (2.1.1) to UtViWt t (2.1.6) can be introduced to the model using data fields of sea surface temperature, fresh water, salt or heat fluxes, surface pressure, wind velocity fields, and initial velocity fields. Different "physics packages" are also available for the ocean, including for example: • Gent/McWilliams/Redi sub-grid-scale (SGS) parameterization of geostrophic eddies (Redi, 1982; Cox, 1987; Gent and McWilliams, 1990; Danabasoglu and McWilliams, 1995). • A biogeochemical package to model the cycling of carbon, phosphorous and oxygen in the ocean (McKinley et al., 2000; Follows et a l , 2002). • Ocean vertical mixing - the nonlocal K-profile parameterization scheme K P P (Large et al., 94). • A thermodynamic sea ice package (Bitz and Lipscomb, 1999; Winton, 2000). 2.1.4. Dissipation. In general, the dissipation of momentum is done through Laplacian and biharmonic friction: dv + A -^ 2 (2.1.8) D = A V\v v h v +A V v A A h 18 2.1. M O D E L DESCRIPTION where Ah and A are (constant) horizontal and vertical viscosity coefficients, = d /dx 2 v 1 + d /dy 2 2 is the horizontal Laplacian operator, A4 is the horizontal coefficient for biharmonic friction, and = d /dx 4 i + d /dy A 4 is the biharmonic operator. These coefficients are the same for all velocity components, and can become variable in space and time when the K P P package is activated. Similarly, the mixing and dissipation of salinity and temperature is represented as: here r represents an arbitrary active tracer (salinity or temperature), Kh and K are the horizontal v and vertical diffusion coefficients for salt and temperature, which can be treated as variables (again, when using the K P P package) or used as constants. These coefficients are the same for all tracers. 2.1.5. M e t h o d of solution. Prognostic equations for flow velocity, u and v, salinity, S, and temperature, 9, are solved by stepping forward in time (2.1.1), (2.1.2), (2.1.5), and (2.1.6), respectively. Diagnostic equations are solved for vertical velocity, w, density, p, and for pressure, p, under H P E and Q H dynamics. The model solves a 2-D elliptic equation to find the surface pressure and the hydrostatic pressure at any level is computed from the weight of fluid above. Under N H dynamics a 3-D elliptic equation must be solved for the non-hydrostatic pressure before the horizontal momentum equations are stepped forward; here w is found diagnostically from the vertical momentum equation. 2.1.6. Discretization. To solve the prognostic equations while satisfying the constraints imposed by the diagnostic equations, the model has several different formulations. However, for brevity, only the formulation that was finally used is going to be described. This formulation uses a semi-implicit pressure method for the hydrostatic equations with an implicit linear free-surface and variables staggered in time. 2.1.6.1. Temporal discretization. The horizontal momentum equations (2.1.1), and (2.1.2) and the continuity equation (2.1.4), can be expressed as: (2.1.9) d u + gd r) - G.u (2.1.10) dv t t x + gdyT] - G v 19 2.1. M O D E L DESCRIPTION (2.1.11) du + dv + dw = 0 x y z where subscripts indicate partial derivatives, and the G terms summarize all the terms in the momentum equations except the temporal and surface pressure gradient terms. From the continuity equation integrated over the depth of the fluid, H, the linear free surface equation becomes: dtV + d Hu + dyHv — 0 x here Hu = J Hudz, and similarly for Hv. H The horizontal momentum equations and the free-surface equation are discretized as: (2.1.12) u (2.1.13) v (2.1.14) 7? n+1 + Atgd rj n+1 x n+1 + Atgd Tf + AtG ^ =v AtG^ n n +1 y n + 1 =u + M8 Hu^+ i x + + AtdyH^ { +1/2) +1/2) = rf where n represents the time level, and At is the time step. The implicit terms on the left hand side of (2.1.12) to (2.1.14) are the ones calculated at time level n + 1; the implicit backward time stepping scheme is being used. A l l other terms on the right hand side of the equations are explicit in time. Using the the superscript * to represent a prediction for the temporary flow variables at time level n + 1, let: (2.1.15) u* = u + AtG^+W v* = v + AtG^+W n and (2.1.16) n 20 2.1. M O D E L DESCRIPTION so (2.1.17) u = u* - n+1 &tgd 7f +l x and (2.1.18) v n+1 = v* - Mgd r) n+l y Then let (2.1.19) V* = V ~ Atd Hu* n x Atd Hv* y Substituting (2.1.17), (2.1.18) into (2.1.14) and then using (2.1.19) gives: (2.1.20) d gHd if +1 x x + d gHd ri +^ n+1 y y = Equations (2.1.15) to (2.1.20) are solved in the order (2.1.15), (2.1.16), (2.1.19), (2.1.20), (2.1.17) and (2.1.18). This sequence represents the algorithm with a backward implicit, linearized free surface. The implicit treatment of the free-surface allows the flow to be divergent and for the surface pressure to respond on a finite time-scale (as opposed to instantly). In MITgcm, the equations for the free surface are implicit in time, but the equations for the horizontal velocity (u, v), and the active tracers (0, S) use the scheme: (2.1.21) T* = T " + AtG j { + 1 / 2 ) where T is any of the variables u, v, 6, S. The term at level n + 1/2 is extrapolated using the quasi-second order Adams-Bashforth method (Hall and Watt, 1976) (2.1.22) 4 n+1/2) = (3/2 + e )G^ AB - (1/2 + where CAB is a small value that stabilizes the method. n-1 e )G^ AB 2.1. M O D E L DESCRIPTION 21 2.1.6.2. Algorithm. Due to the choice of advective scheme (to be discussed in the next chapter), the thermodynamic variables are staggered in time with respect to the flow variables. This is useful for well stratified problems where internal gravity waves can be the limiting factor in choosing a stable time step, because the staggering of the variables allows the gravity wave terms to leap-frog in time providing second order accuracy and more stability. Finally, the complete MITgcm algorithm for an implicit linear free-surface, with explicit and staggered Adams-Bashforth time stepping is as follows: (1) Integrate the hydrostatic pressure (note this does not include the surface pressure effect) (2.1.23) = j ° b (6 '\ p ^ n /2 n+1 h S ') dz n+1 2 where b=-^-(p(e,S,z)-p ) ref Pre} is the buoyancy. (2) Introduce the new variables Gl"-, G[™\ (and G ^ if non-hydrostatic) to represent all the terms in the momentum equations except the surface and hydrostatic pressure terms, and evaluate them from time level n (2-1.24) G^ = 2 ^ (ft?) (3) Extrapolate the tendencies in time (2.1.25) (2.1.26) G^ 1 / 2 = (3/2 + e ) ) G^ AB ^ = Itl + At ( G ^ V - (1/2 + _ e A B )fe y^ ) 2 (4) Calculate the right hand side of the surface pressure equation (2.1.27) - 77* = if AtV • JflP* - = - i & (5) Find the surface pressure (2.1.28) V - ^ + 1 ^ 1 2.1. M O D E L DESCRIPTION 22 (6) Update the horizontal velocities (2.1.29) 17£ = -v* - AtgVr) n+1 +1 h (7) Integrate the continuity equation for w (2.1.30) w n+1 =~ f {d u + dv ) n+1 dz n+l x y (8) Calculate the advection of the thermodynamic variables using the latest velocities (2.1.31) G^ + 1 / 2 ) = G ,s (v^\ e 0 " ( + 1 / 2 ) , S^ / ^ 1 2 (9) Extrapolate the tendencies forward in time (2.1.32) Gf n e (2.1.33) = (3/2 + e )G + n AB l/2 s ( > + / , 5"+ / ) - [e l , 3 2 3 2 n+1 2 - (1/2 + S /) n+1 2 e )G^~ 1/2 AB + AtG^ 1 (10) Update indices n to n + 1, and repeat the process. 2.1.6.3. Spatial discretization. The spatial discretization uses the finite volume method, where a second order centered difference scheme is used for the interior but allows boundaries to intersect a regular grid allowing a more accurate representation of the position of the boundary. The horizontal and vertical directions are treated separately. Notation. The notation used to describe the spatial discretization in the following pages is summarized here: General • Ax, Ay, Az: grid spacing in x, y, z directions. • A : Area of the face orthogonal to "o" direction (o = u, v, w...). Q • V , V , V , VQ: Volume of the grid box surrounding the u, v, w, 6 points. U V W • i, j, k: Current index relative to x, y, z directions. Operators • 8i: 5ip = Pi+i/2 • Pi-l/2 : l i _ (Pi+l/2 +P1-1/2) 2.1. M O D E L DESCRIPTION 23 • S: x S = ±8* x • V = horizontal gradient operator : Vp = {5 p, 5 p} x y • V - = horizontal divergence operator : V • ~v = {8iAyu + SjAxv} where A is the horizontal area of the corresponding cell. • V = Laplacian operator : V p = V • Vp 2 2 Grid. Figure 2.1.1 shows the space staggering of the flow variables using the Arakawa C-grid. The zonal velocity (u) points are located at the middle of the interface between continuity cells in the zonal direction, the meridional (v) and vertical (w) velocity components have the same configuration. The active tracer (9,S) points are located at the middle of the continuity (or tracer, T) cell. Horizontal grid. The model uses a horizontal coordinate system that is curvilinear orthogonal. In this system, the grid can be Cartesian, spherical polar or general curvilinear. Each cell in the horizontal grid is described by the length of its sides and the area of the cell. Figure 2.1.2 depicts the lengths and areas related to the different active variables in the model, and their orientation with respect to each other. The notation is as follows: • The "g" suffix indicates that the lengths are along the grid boundaries. • The "c" suffix associates the quantity with the cell centers. • The suffix associates points with the vorticity points. • The "v" suffix indicates that the length is measured between the w-points. • The " / " suffix indicates that the length is measured between the (tracer) cell faces. • The "to" suffix associates points with the u-points (w stands for west). • The "u" suffix indicates that the length is measured between the u-points. • The "s" suffix associates points with the w-points (s stands for south). Vertical Grid. As for the horizontal grid, the suffixes "c" and "f" indicate faces and centers. Figure 2.1.3 shows the default vertical grid used by the model. The descriptor Azj is the difference in meters between the faces (i.e. Az/ = — S^z where the minus sign appears due to the convention that the surface layer has index k = 1.); similarly Az is the distance in meters between cell centers. c 2.1. M O D E L DESCRIPTION 24 FIGURE 2.1.1. Sketch of the staggering of the variables using the Arakawa C-grid. Topography. As it was mentioned at the beginning of this chapter, MITgcm uses the intersecting boundary method to represent the bottom topography. This method allows the boundary to intersect a grid of cells modifying the shape of the cells that are intersected. The topography is then a piecewise constant (partial step) representation of the bottom boundary. Figure 2.1.4 shows a cell that is in contact with the bottom topography in the x-z plane. The physical thickness of a tracer cell is given by h (i,j,k)Azf c side is given by h (i,j,k)Azf w (or h (i,j,k)Azf s and the physical thickness of the open in the y-z plane). The non-dimensional factors h , c h , and h , represent the fraction cell depth that is "open" for fluid flow, these values range between s w 0 and 1. Continuity and vertical momentum equations. With the staggered C-grid discretization, the continuity and vertical momentum equations are summarized as 25 2.1. M O D E L DESCRIPTION FIGURE 2.1.2. Staggering of horizontal grid descriptors (lengths and areas). The shading indicates the cell boundaries and is the reference cell for all panels, a) The area of a tracer cell, A , is bordered by the lengths Ax and Ay . b) The area of a vorticity cell, A^, is bordered by the lengths Ax and Ay . c) The area of a u cell, A , is bordered by the lengths Ax and Ayf. d) The area of a v cell, A , is bordered by the lengths Axf and Ay . c g c w g c v s u (2-1.34) e nh (2.1.35) ( d w + -^—5 p ) = e G V Azc ) t k A d r) + c\ c nh Ay Azfh u t g w k nh + b - -^—S p Az k w k c + Sj ^ Ax Azfh v g s —0 k where p h is the non-hydrostatic pressure. In the hydrostatic limit e h — 0, and G„; — 0, then n n (2.1.34) is reduced to the hydrostatic balance. Flux form of the momentum equations. The finite volume model is based on the Eulerian flux form of the momentum equations. The G (generation and dissipation) terms are divided into advective, Coriolis, horizontal dissipation, vertical dissipation and curvature terms: 26 2.1. M O D E L DESCRIPTION Az vv f \ w z W F I G U R E 2.1.3. Vertical grid in MITgcm. The tracer centers are placed at the middle point between cell interfaces. The vertical velocity points (w) are placed at the middle point of each cell interface. (2 1 36) G—* _|_ QCOTI _|_ Qh—diss QV—diss _|_ ^metric _|_ Qnh—metric For brevity only the first four terms which are the terms used in the experiments are will be described. For the description of the other terms the reader is referred to the MITgcm user's manual, or the MITgcm web site (http://mitgcm.org/). • Advection of momentum: The advective operator is second order accurate in space, and is expressed as (2.1.37) A Az h G (2.1.38) AsAzfhsGf" = <wV + 8jV v a dv w f w u = SiU^ + SjTv? + <S WV fc j j + 5 Wv j k k 2.1. M O D E L DESCRIPTION FIGURE 2.1.4. Sketch of the x-z plane showing the location of the non-dimensional fractions h and h . The thickness of a tracer cell is given by h (i,j,k)Azf and the physical thickness of the open side is given by h (i,j,k)Azf. The bold line represents the real bathymetry. The shaded region represents the model bathymetry. c w c w (2.1.39) A Az Gl dv c c = < W V + 8jV w k j where the volume fluxes U, V, and W are defined as (2.1.40) U= Ay Az h u (2.1.41) V = Ax Az h v (2.1.42) g f g f W = Aw c w s + 8Ww k k k 2.1. M O D E L DESCRIPTION 28 • Coriolis terms: The discretization of the Coriolis terms used is: (2.1.43) Cori _ f—j °™ = -fu " l G ^ ° " = fv- (2.1.44) 1 where the Coriolis parameter / , is set to a constant for the present work. • Lateral dissipation: The viscous stresses are discretized as (2.1.45) A Az h G^~ = 5iAy Azfh Tn + 8jAx Azfh^i (2.1.46) A Az h G^- = diAy Az h^T + 5jAx Az h T diss w f w f diss s f s c u f 21 v 2 f f C 22 and the lateral viscous stresses are discretized as (2.1.47) TH = A -^—8 (2.1.48) na = A ~-6 (2.1.49) rai = (2.1.50) r 2 2 h iU h jU Ah-^—SiV = A -^-SjV h Vf A Here the Laplacian viscosity coefficient, Ah, has units of m s 2 _ 1 and for the laboratory case is the molecular viscosity of water. The free-slip boundary condition sets the stress components against the wall to zero, while for the no-slip boundary condition (flow zero at the boundary) the boundary stresses are added as an additional body force of the form: (2.1.51) (2.1.52) sUe-dra G 9 = 4 _ Azf QSide-drag = ^ £AZf A * / Ay u ( J _ Ay„ LSX V 2.1. M O D E L DESCRIPTION 29 where the fraction of open water at the vorticity point, 7i£, is denned as: ) = fc m i n l n:{hj,k), h (i,j h - l,k),h (i,j, w k),h (i s s - k)} • Vertical dissipation : Vertical viscosity terms are discretized with only partial adherence to the variable grid lengths introduced by the finite volume formulation. This reduces the accuracy of these terms to first order directly adjacent to boundaries; where the bottom drag terms appear. A solution to this problem is proposed in section 2.3.2. The discrete form of the vertical dissipation terms is (2.1.53) Gr*- = -x^j-SkTu = (2.1.54) (2.1.55) ^-S r k 23 cr"" = t n ^ h m where the interior vertical stresses are (2.1.56) ri3 = A -^—5 u (2.1.57) T = A -^-5 v (2.1.58) r 2 3 3 3 = v k v k A -^—6 w v k As for the lateral viscous terms, the free-slip boundary condition is equivalent to simply setting the stress to zero on the boundaries. The no-slip boundary condition is implemented as an additional term acting on top of the interior and free-slip stresses. Bottom drag represents additional friction, in addition to that imposed by the no-slip condition at the bottom. The drag acts as a stress, and is expressed as a linear function of the mean flow in the layer above the topography: (2.1.59) bottom-drag T = (^A ^ v +r ) 6 U 2.2. A D V E C T I O N SCHEMES IN M I T G C M 30 (2.1.60) where r\, is the linear bottom drag coefficient with units of m-s The model has an alternate form that is available to treat the momentum equations. The use of it is recommended when working in curvilinear coordinates because it avoids problems generated by the earth curvature terms in (2.1.1), (2.1.2) and (2.1.3). This form is invariant under coordinate transformations so it can be applied in any orthogonal curvilinear coordinate system such as spherical coordinates, boundary following coordinates or the conformal spherical cube coordinate system. However, several attempts to use this form of the equations were performed and all of the results produced instabilities inside the canyon. Because of this, the vector invariant form of the momentum equations is not going to be described here. In the particular case of this study, the curvature terms should not represent any problem since the size of the domain is too small to be affected by the earth's curvarture, and also because the geometry information is already stored in the lengths and volumes used for the discretization of the equations. The discussion of the equations for the active tracers and the different advection schemes is crucial for the present work. For this reason, the advection schemes are going to be discussed more deeply in the next section. 2.2. Advection Schemes in MITgcm The active properties (temperature and salinity) are immersed in a loop in which they are transported by water currents, and in turn they affect the density distribution, which affects the pressure field, and this pressure field is used in the calculation of currents once again. For this reason, the accuracy of, and problems related to, the advection of active tracers have repercussions on the overall solution of the flow field. There are several advection schemes that can be used in MITgcm. Each one of the schemes could be more useful to investigate a particular problem or could be more computationally effective depending on the temporal and spatial resolution, and limitations particular to the specific task at hand. In this section, the basic discretization of the advection-diffusion part of the tracer equations and the various advection schemes available in MITgcm will be described. 2.2.1. Preliminary concepts. 2.2. A D V E C T I O N SCHEMES IN M I T G C M 31 2.2.1.1. Analytic solution of the advection equation. For illustration purposes, it is useful to describe the analytical solution of the basic advection equation without diffusion or forcing. This equation can be expressed as: (2.2.1) d 9 + V • {~v9) = 0 t where 9(x,y,z,t) is the potential temperature, but could be any arbitrary tracer, and ~v is the velocity vector. Simplifying (2.2.1) for the one-dimensional case we have (2.2.2) d 6 = -d u6 t x Equation (2.2.2) describes the variation of a tracer transported by a velocity field. For u constant, if 9(x , t — 0) = f(x ), the analytical solution of (2.2.2) has the form 6(x, t) = f(x — ut). 0 0 Substituting a wave-like solution 6(x,t) = 9(t) exp(ikx) into (2.2.2), yields the equation: dt9 = iku9, where the phase speed and the group speed are both u. This velocity is not dependent on the wave number k and therefore there is no dispersion in the solution of (2.2.2). Due to this, 9(x,t) should maintain its original shape, but displaced to x + ut after the tracer is advected for time t. 0 2.2.1.2. Problems related to the discretization . When the advection equation is discretized in time and space, several problems and approximation errors arise. The following is a brief description of some of these errors. Stability: Because of the similarity between the advection equation and the wave equation, the stability of the former can be investigated by analogy with results obtained for the latter. In this way, it can be demonstrated that explicit time integration schemes are generally stable if the grid spacing is greater than the distance that the fastest waves can travel in a model time step (i.e. Aa; > uAt). This condition is known as the C F L (Courant-Friedericks-Levy) criterion (Press et al., 1992). Phase and damping errors: Discretizing the the advection equation (2.2.2) using centered differences gives: nn+l _ an-l (an an \ 2.2. A D V E C T I O N SCHEMES IN M I T G C M 32 where j is the spatial index, and n is the time index. Using a wave-like solution to (2.2.2) of the form 0™ oc exp [i(kjAx - nAw)] where ALJ is the phase change of the approximate solution per time step, it can be shown that AUJ ~ /(sin(fcAa;)). The phase error is the difference between Au> and the true phase change ukAx, and it represents the lag between the finite difference solution and the true solution. Similarly, the damping error is the difference in the magnitude of the true solution and the finite difference solution. These errors are easy to avoid by the implementation of leap-frog based schemes (Haidvogel and Beckmann, 1999). As staggering in time means that the thermodynamic variables are effectively updated using a leap-frog scheme, numerical damping should not be a problem in the simulations. Dispersion error and false extrema. The dispersion error is closely related to the phase error, which is a function of the wave-number. Therefore waves that initially have the same phase but different wave number will gradually disperse. As a result ripples appear in the solution, producing false extrema (under and over shoots) in the advected quantity. Time-splitting errors. Time splitting errors appear when the advective scheme uses more than two time levels (i.e. more than n + 1 and n). In this cases, the difference equation can have two or more solutions, with different amplification factors (e.g. if the solution of a tracer equation is a function of exp(—kAt), and the other solution is a function of — exp(fcAi)). Note, that a negative solution (known as the computational mode) for quantities such as salinity or a passive tracer concentration gives an unrealistic solution. Also, depending on the advective scheme, this problem may bring divergence of the solution at odd and even time steps. Boundary condition errors. These errors occur when the difference equation requires boundary conditions different than the ones required by the partial differential equation. For example, the upwind differencing in space scheme requires the same boundary conditions as the original partial differential equation, while the centered in space schemes require special treatment at the downstream end. Also, higher order centered schemes often need extra data at the boundaries. Aliasing error. With non-linear advection equations, non linear terms of the form ud u6 appear. x When these terms are evaluated in a finite difference scheme, errors occur due to the inability of 2.2. A D V E C T I O N SCHEMES IN M I T G C M 33 the discrete grid to resolve wavelengths shorter than 2Az. That is, wavenumbers greater than kmax — K' I Ax Haidvogel and Beckmann (1999) demonstrate that for k > k , the wave will be misrepresented max as k* = 2k max — k. Errors of this sort are called aliasing errors. 2.2.2. T r a c e r e q u a t i o n s . In MITgcm, the default advection scheme is the centered difference second order method which requires a second order or quasi-second order time-stepping scheme to be stable. For this, the model uses the quasi-second order Adams-Bashforth method (ABII). The forced advection-diffusion equation is expressed as: (--) 2 2 ®t0 + G 3 e where G , G ^, and G ^ 8 e adv di o r c =G + G 9 adv e diff forc , summarize the terms due to advection, diffusion, and forcing, and are expressed as (2.2.4) (2.2.5) G e adv = d u6 + d v9 + 8 w6 - 0V • ~v x y G diff = z V-KV6 and the forcing can be some arbitrary function of state, time and space. The term, f?V • ~v , is the surface correction term that is needed to retain local conservation in the surface layer; this term is zero everywhere else as the interior flow is non-divergent. Subsequent discussion will ignore this term until multidimensional advection is considered in section 2.2.4.4. By default, the model uses the method of lines in which the space and time discretizations are treated separately. The tendencies are calculated at time levels n and n — 1 using the AdamsBashforth method described in section 2.1.6.1. However, the Adams-Bashforth method for the time stepping is only used when the standard second, third, and fourth order linear advection schemes are selected. For the rest of the advection schemes available, the model uses the forward Euler method for the time stepping. 2.2.3. L i n e a r A d v e c t i o n s c h e m e s . The centered second order, centered fourth order, first order upwind and upwind biased third order advection schemes are known as linear because the 2.2. A D V E C T I O N SCHEMES IN M I T G C M 34 coefficients for interpolation of the advected tracers are linear and are a function only of the flow, not the tracer field itself. 2.2.3.1. Centered second order. This discretization is designed to be consistent with the continuity equation to facilitate the conservation of properties. However, this scheme is known to be noisy, meaning that short waves pollute the solution. For this reason, a finite amount of diffusion is added using (2.2.5), and as the shortest waves are more effectively dissipated than longer waves, a better solution is produced. The advection operator is discretized as: (2.2.6) A Az h G a dv c f c 9 = SiF + 8jF + 5 F x y k z where the F terms are the area integrated fluxes: (2.2.7) F = U9 (2.2.8) F = V9 (2.2.9) F = W6 l x J y k z Adcroft et al. (1997) demonstrate that for non-divergent volume flux this discretization conserves the tracer locally and globally. 2.2.3.2. Third order upwind bias advection. Upwind schemes take the direction of the flow into account in order to approximate the spatial derivatives with more accuracy. These schemes are conditionally stable and monotonic, but the space-differencing errors that arise tend to damp the solution (Haidvogel and Beckmann, 1999). To reduce spurious damping and phase errors, the upwind biased third order advection scheme improves the difference formulation using points that are ±2Ax away, and then expanding the derivative in a Taylor series. In MITgcm, this scheme takes the form: (2.2.10) F = U X 6 2.2. A D V E C T I O N SCHEMES IN M I T G C M 35 (2.2.11) (2-2.12) F = w(e-^8 d^ z + kk -^\W\5 9 kkk . This scheme offers a good balance between smoothness and accuracy (Wicker and Skamarock, 1998), and although false extrema are permitted, the amplitude is much reduced compared to the centered second order schemes (Pietrzak, 1998). The boundary condition for this scheme is 6^6 = 0 at the boundaries. Currently, the model assumes that 5u, and 5m are also equal to zero near the boundaries. 2.2.3.3. Centered fourth order advection. Centered fourth order advection is the most accurate of the linear schemes implemented in MITgcm. Fourth order schemes improve the numerical quality of the solutions at low computational cost (Shchepetkin and McWilliams, 1997). However, the scheme is noisy and prone to show false extrema like the centered second order method, and so it must be used with some finite amount of diffusion. The centered fourth order fluxes are discretized as: (2.2.13) F = U[6-^8 8 (2.2.14) F = V[9-±5 6 (2.2.15) F = w(o-±6 e X y z U jj kk The boundary conditions are the same as for the third order scheme. 2.2.4. Non-linear advection schemes. The term "Non-linear advection schemes" is related to the use of non-linear coefficients for the interpolation of the advected tracer, not to the nonlinearity of the advection equation. As was described above, higher order schemes are in general more accurate, but do not preserve monotonicity (produce false extrema in the tracers), and can have problems with time-splitting errors (preservation of the sign of the tracer). To control or eliminate these problems, flux limitation algorithms are introduced. These are not advection schemes, but rather extensions to the schemes 2.2. A D V E C T I O N SCHEMES IN M I T G C M 36 that help to achieve or maintain certain properties of the numerically advected tracer, often playing a role similar to the addition of artificial diffusion in the linear schemes. The methods that use flux limitation algorithms are called flux limited advection schemes. The underlying schemes for the non-linear methods that are going to be described next, are the first order upwind, the second order centered difference Lax-Wendroff and the third order upwind biased. It must be remembered that when employing flux limited schemes, first order upwind or directspace-time methods the time-stepping of the tracers is switched to being forward in time instead of the Adams-Bashforth method. 2.2.4.1. Second order flux limiters. When using the second order flux limiter method in terms of a first order upwind flux and second order Lax-Wendroff flux, the limited flux is given as: (2.2.16) F = F+ ^(r)F 1 LW where ip(r)is the limiter function, F i = uO - i \u\ 8i0 (2.2.17) l is the first order upwind flux, \u\ (2.2.18) F LW =F + ^(l-c)5 e 1 is the Lax-Wendroff flux, and c = | u | At/Ax i is the C F L number. The limiter function ip(r), is a function of the slope ratio (2.2.19) r = ^~ 1 / 2 0i+l/2 (2.2.20) ~ ~ 0i _ Vu >0 3/2 "i-1/2 r = 7 " " Vu<0 "i+1/2 - "i-1/2 3/2 m / 2 The limiter function used is the Superbee limiter (Roe, 1985): (2.2.21) ib(r) = max {0, min [1,2r], min [2, r}} 2.2.4.2. Third order direct space time. When describing the tracer equation in section 2.2.2, it was stated that the model uses by default the "Method of Lines" in which space and time are treated 2.2. A D V E C T I O N SCHEMES IN M I T G C M 37 separately. The alternative is the direct-space-time (DST) method, which deals with space and time discretization together. Using the upwind biased third order DST scheme (Hundsdorfer and Trompert, 1994; Pietrzak, 1998), the flux is given by: (2.2.22) F = u (0;_ (2.2.23) F — u {6 1/2 i+1/2 + di (0 i+1/2 - 8i_ ) + d ( f l i _ i+1/2 - 9i_ ) - d (9 - di (6 1/2 1/2 2 2 1 / 2 - 0.-3/2)) V u > 0 - 9 )) i+3/2 i+1/2 V«<0 where (2.2.24) d = i (2 - |c|) (1 - |c|) (2.2.25) d2 = Ul- a \c\) (1 + \c\) 0 Note that when the C F L number, c, vanishes, d\ and d approach 1/3 and 1/6 respectively, 2 reducing (2.2.22), and (2.2.23) to the third order upwind method described by (2.2.10), (2.2.11), and (2.2.12). The third order direct space time scheme is stable for 0 < c < 1, and interestingly the accuracy increases with the C F L number. However, it still is susceptible to the production of false extrema. 2.2.4.3. Third order direct space time withfluxlimiting. The appearance of false extrema in the third order DST method can be controlled with the use of a flux limiter. The flux term becomes then: (2.2.26) F=\(u+ where \u\) [9i_ 1/2 + *l> (r+) (6 i+1/2 - 6i_ )) +\{u1/2 \u\) (8 i+1/2 - ip (r~) (6 i+1/2 - ^_ )) 1/2 2.2. ADVECTION SCHEMES IN MITGCM 38 and the limiter is the Sweby limiter (Sweby, 1984): (2.2.29) yj(r) = max JO, min 1-c min (1, di + d r), r 2 2.2.4.4. Multidimensional advection. When the advection equation has more than one dimension, advection schemes that were robust in one dimension can develop problems such as loss of monotonicity. In MITgcm these problems are fixed using a variant of the conventional flux splitting methods that allows the flux calculations to be implemented as if in one dimension: (2.2.30) 6 / n+1 = 6 -At 3 n (^S F (9 ) + ^ ^ . ^ (^SjFy (> / ) + 0 ^j«) (> ff^-l-SkW^j t n x (2.2.31) 6 l =9 ' - At (2.2.32) 6 =6 - At {^-JkF n+2 3 n+z/z n+l 3 n+2lz z +1 + 2 3 / ) 3 + n The way this method is incorporated into the general algorithm in MITgcm is by updating the advective tendency, G i , to n + 3/3, while the diffusive and forcing terms use the n time level. ac v This is: (2.2.33) Gl 1 / 2 = ± (<r«/3 - 0") and the complete time-stepping is: (2-2.34) c? n+1 = 6 + At (C?:+ + G n 1/2 diff (0") + G ) n forcing 2.2.5. Comment on the advection schemes. From this section it can be concluded that the most desirable advection scheme should include monotonicity, low implicit diffusion, accuracy, and low computational cost. However, these properties are not likely to be achieved at the same time. If a scheme offers monotonicity, it may lose accuracy; higher order approximations are accurate but more computationally expensive; non-linear schemes are more difficult to implement, and so on. Due to the nature of the present work, the advection scheme that is used is the multidimensional third order direct-space-time with flux limiting for the following reasons: 2.3. EXPERIMENT CONFIGURATION 39 • High resolution is needed due to the spatial and temporal variability of the physical processes in a rotating tank with a submarine canyon, and high order schemes offer good accuracy. • With non-linear schemes, as the CFL number increases the accuracy increases, and in consequence the limitation imposed by the advection on the time-step is reduced. • False extrema in the advected tracers is a non-physical process that should be avoided and this is main reason for the choice of the third order DST with flux limiting. However, a comparison between different advection schemes and their properties when applied to the model of the rotating tank is discussed later in Chapter 3. 2.3. Experiment configuration The configuration of the numerical experiments using MITgcm is mainly established by the details of the laboratory experiment where the flow is accelerated by changing the rotation rate of a tank filled with a stratified fluid that is initially in solid body rotation. This forcing is simulated in the numerical model using a body force exerted during one rotation period to match the time over which the rotation rate of the table is changed. The complete details about the configuration of the laboratory experiment to simulate an upwelling event in a rotating tank are described in Allen et al. (2003). The numerical experiments using MITgcm used the same length scales employed in the laboratory. However, when using the model SCRUM, Allen et al. (2003) increased the horizontal length scale by a factor of 10, as the slopes in the tank canyon were too steep for this model to accurately compute the horizontal pressure gradients. For this reason, in order to compare the MITgcm results with the results obtained with SCRUM, one experiment is conducted increasing the horizontal length 10-fold with a corresponding adjustment of parameters to maintain the relevant non-dimensional numbers. 2.3.1. Laboratory Configuration. The laboratory experiments were conducted on a computer controlled rotating table of 1 meter diameter, at the University of British Columbia. The tank representation of the coastal bathymetry included a flat abyssal plain, a steep continental slope (45°) and a continental shelf (slope of 5°). The continental shelf and slope extend inward from the outer wall of the tank. The shelf-break is indented by a canyon of 8 cm length (L) from the 2.3. E X P E R I M E N T CONFIGURATION 40 F I G U R E 2.3.1. Schematic of the bathymetry in the laboratory tank. The radius of tank is 50 cm. The "shelf slope is 5° and the "slope" slope is 45°. The numerical geometry is similar except that only 60% of the domain is modeled. Reproduced from Allen et al. (2003) with permission of Dr. S. E. Allen. head to the mouth, a width at the shelf break (W b) of 6.9 cm, and a width at mid-length (W) of a 2.4 cm (Figure 2.3.1). The tank was filled while it was rotating at / = 0.4 s , and then it was forced by accelerating _1 the table to / = 0.52 s linearly over 27.3 s. This is about one rotation period of the tank and is - 1 comparable to the duration of an upwelling event in the real ocean. With the change of rotation rate the fluid moves relative to the tank with the shelf to its left. The velocity at the shelf break is approximately 1.2 cm-s . -1 The flow near the canyon was visualized using neutrally buoyant particles. These particles were illuminated by two horizontal light sheets of 1 cm thickness, that covered depths from 1.2 cm to 2.2 cm (green light) and from 2.2 cm to 3.2 cm (red light). A video camera recorded the motion of the particles at intervals of 0.5 s, 2.7 s after the acceleration of the tank was complete (27.3 s) for a total duration of 5 seconds. This experiment was repeated 5-15 times. Comparing close tracks in different runs, Allen et al. (2003) calculate a velocity error for the laboratory experiments of 0.2 cm - s . -1 2.3. EXPERIMENT CONFIGURATION 41 2.3.2. Numerical model configuration. 2.3.2.1. Forcing. To match the acceleration of the flow inside the tank produced by the changes in rotation rate, a body force is introduced as an external forcing to the model's zonal momentum equation (F term on the right hand side of (2.1.1)). This body force is uniform in the vertical u direction, and increases radially toward the outer wall of the tank. The body force is applied during approximately one pendulum day (4TV//). The forcing increases from zero to its maximum over the first 0.6 s, then is held constant for 27.3 s, and then set back to zero again. The total time of the simulation is 35 s. The body force is discretized as: (2.3.1) Fij = a-t(l- (2.3.2) F = -t [l-°id (2.3.3) a 1 °- U™*-fi\ <ti 8 8 U m a x Jmax _, j ) J- it ) ,h<t<t 2 F^ = 0 , t < t < t 2 end where Uf re a — h-t— 2 is the slope of the forcing function , u f — 0.024 m -s is a reference velocity, t\ — 0.6 s, is the _1 re time of the initial increase in forcing, t = 26.7 s, is the time over which the forcing is held constant, 2 tend — 35 s is the total time of the simulation, j is the index in the radial direction, and j x ma is the maximum index in the same direction. The units of the F terms are the same as for the G terms, that is m-s~ . 2 2.3.2.2. Domain. For the initial simulations, the model domain included 240 points in the azimuthal direction, 82 points in the radial direction, and 32 levels in the vertical. This results in a spacing at the shelf break of approximately 0.5 cm in the azimuthal and radial directions (Figure 2.3.2), and of 0.28 cm in the vertical direction (Figure 2.0.1). The domain for the numerical simulation only covers 0.6 x 27T in the azimuthal direction in order to save computational time and memory requirements. This reduction does not cause significant problems because the model is set to run with periodic boundary conditions for the bottom topography, and the time required for a disturbance to travel around the domain is much bigger than the time of simulation. That is, a disturbance propagating with the flow at a velocity of 1.2 c m s - 1 2.3. E X P E R I M E N T C O N F I G U R A T I O N 42 0.34 0.32 03 0.28 0 26 0.24 022 02 0 18 0.16 0.1 FIGURE 2 . 3 . 2 . 0.15 0.2 0.25 0.3 T o p v i e w o f the h o r i z o n t a l g r i d over t h e c a n y o n . D i s t a n c e s are i n meters. T h e b a t h y m e t r i c c o n t o u r s are p l a c e d at 1 c m i n t e r v a l s . at the shelf b r e a k ( ~ 3 0 c m f r o m the center o f the t a n k ) , w o u l d p r o p a g a t e a r o u n d the n u m e r i c a l d o m a i n (0.6 x 27r) i n a p p r o x i m a t e l y 47 s, a n d the t o t a l d u r a t i o n o f t h e e x p e r i m e n t is 35 s. T h e m o d e l b a t h y m e t r y was p r o v i d e d b y M i k e D i n n i m a n f r o m O l d D o m i n i o n U n i v e r s i t y , a n d i t is o b t a i n e d f r o m d i g i t i z e d t a n k d e p t h c o n t o u r s ( F i g u r e 2.3.3). T h e m o d e l d o m a i n has a n i n n e r w a l l 10 c m f r o m the center o f the t a n k , a n d the o u t e r w a l l is 50 c m f r o m t h e center. I •Bo.01 I 0.04 -Tr~ 0 05 0.06 BSo.07 H o 08 ^0.09 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 FIGURE 2 . 3 . 3 . T o p v i e w o f the b a t h y m e t r i c c o n t o u r s i n the m o d e l d o m a i n . T h e u n i t s are i n m e t e r s . T h e d o m a i n covers 0.6 x 27T i n t h e a z i m u t h a l d i r e c t i o n i n o r d e r to save c o m p u t a t i o n t i m e a n d m e m o r y r e q u i r e m e n t s . T h e b o u n d a r y c o n d i t i o n s for the m o m e n t u m e q u a t i o n s are n o - s l i p ( S e c t i o n 2.1.2) b o t h at b o t t o m a n d at the l a t e r a l w a l l s . the 2.3. E X P E R I M E N T C O N F I G U R A T I O N 43 2.3.2.3. Dissipation. In the initial runs with MITgcm, the results obtained near the bottom topography presented sudden changes in the variable fields (mostly in the azimuthal, u, velocity field). By analyzing each term of the momentum equations and searching for inconsistencies, the source of these abrupt changes was found to be the vertical dissipation terms. As described in Section 2.1.6.3, when discussing the discretization of the interior vertical stresses and bottom drag terms, the vertical dissipation terms are discretized with only partial adherence to the variable grid lengths introduced by the finite volume formulation. For example, consider the zonal component of the interior stress from (2.1.56) and the water column in contact with the bottom topography as depicted in Figure 2.3.4. For the calculation of the stress between levels k, and k + 1, the expression = A,—— S u Tl3 Az k c is correct to second order accuracy, but for the calculation between levels k + 1, and k + 2 , Az c does not represent the correct spacing. This should be instead Az (l + h) 2 c w which leaves the interior vertical stress as (2-3.4) r 1 3 =A Su Az (1 + h ) 2 v k c w similarly for the meridional component of the interior stress. The bottom stress needed a similar adjustment. Considering only the no-slip condition without additional drag, the bottom drag term in the zonal direction in MITgcm is ^° - s ttom dra T _ (2A /Az ) V c u, because the bottom drag terms are only active in the cells that are in contact with topography, this expression is not completely correct. Instead we modify it to take into account the reduced thickness of bottom cells to obtain: (2.3.5) bottom-drag T = V 1 t\z • h c v The same is applicable to the meridional term: (2.3.6) bottom-dra T g = ^ 1 Az • h. r 44 2.3. E X P E R I M E N T C O N F I G U R A T I O N w—a Az f ..... ^ ....i y 2 w > ).s | ^ \ S '• VfcAz,h w 1 V FIGURE 2.3.4. Vertical section of a generic water column in contact with the bottom topography. The thickness of the deepest cell is Az • h . c w With these modifications the results were significantly smoother, improving the adaptability of the model to changing topography. 2.3.2.4. Stratification and equation of state. In MITgcm the model equations are written in terms of perturbations, therefore, a reference thermodynamic state needs to be specified. This is done using a reference potential temperature and/or salinity profile that is uniform in the horizontal. However, only the density profile is available from laboratory measurements (Figure 2.3.5). In the model the density is represented as (2.3.7) P = Po + p' where p is a constant density, and p' is the remainder. The variable that is dynamically important Q is p\ and is calculated using the linear equation of state as: (2.3.8) p' - Po (Ps (S) - ag (6)) 2.3. E X P E R I M E N T C O N F I G U R A T I O N 1.005 1.01 1.015 1.02 1.025 1.03 1.035 1.04 1.045 45 1.05 Density (1000 kg m" ) 3 FIGURE 2.3.5. and meters. Density profile used in the numerical simulation. Units in kg-m 3 where (3 is the haline contraction coefficient, ag is the thermal expansion coefficient, and S and 9 S are the reference salinity and temperature. Substituting (2.3.8) into (2.3.7): (2.3.9) P = Po + PoPsS - PoaeO thus (2.3.10) 6 = - (p- Po- po(3 S) s P ae 0 Equation (2.3.10) is used to obtain the reference temperature profile for the experiment runs. The profile used in the numerical simulations is presented in Figure 2.3.6. These temperatures are calculated using p — 1004.3 k g - m , (3 — 7.4x 10~ , ag = 2x 1 0 K , and a constant salinity 5 = -3 0 4 - 3 - 1 S 35. The thermal expansion coefficient, ag, is artificially increased by one order of magnitude from the normal value (2x 1 0 K ) in order to have a reasonable range in the perturbation temperature - 4 - 1 and to keep the numerical simulations stable. Since the same coefficient is used in the simulation code, the density in the model is the same as the density in the tank. 2.3.2.5. Float injection. To compare the model and the laboratory results, float particles were introduced in the model at the same position where particles are observed to start to move in the tank. The vertical position of the particles in the tank can only be estimated, by the color that they reflect (green or red), within 1 cm depth. However, the particles can move differently within the 1 cm increments due to the strong vertical shear that is observed in the tank (to be discussed 2.3. E X P E R I M E N T C O N F I G U R A T I O N 46 Potential Temperature fC) 2.3.6. Potential temperature profile used in the numerical simulation. Units in °C and meters. FIGURE in Chapter 3). For this reason, in the numerical model particles were planted every millimeter in depth for each horizontal position. Thus, if a particle in the tank started to move from the position (xl,yl), and from the green light reflected the depth was determined to be between 1.2 cm and 2.2 cm, in the numerical simulation there will be eleven particles planted at the position (xl,yl) in 1 mm intervals from 1.2 cm to 2.2 cm of depth. In this way, for 32 particles that are traced in the tank, in the numerical model there are 352 particles in total. Finally, to compare the laboratory and the numerical model results, the particles that best matched the final position of the particles in the tank were the ones chosen. The advection of floats in MITgcm was originally conceived to simulate floats that drift at an established depth and then move to the surface at a defined time interval, making profiles of salinity, temperature, pressure, and velocity of the flow. This means that a float had to be placed at a vertical model level where it would be advected horizontally by the currents. The time stepping of the float advection is done using a second order Runge-Kutta scheme (Press et al., 1992), whereby velocities and positions are interpolated between the grid points using bilinear interpolation in two dimensions. However, as the goal of this research is to model upwelling conditions, the floats or particles should be allowed to move in three dimensions. To achieve this, the original code was changed to use 3-D linear interpolation, and to use any initial depth in meters instead of a fixed level number. 2.3. EXPERIMENT CONFIGURATION TABLE 47 2. Physical and numerical parameters for the coarse resolution experiment. Description Vertical eddy viscosity coefficient for mixing Ay of momentum. In the laboratory case it is the molecular viscosity of water. Horizontal eddy viscosity coefficient for mixA ing of momentum. In the laboratory case it is the molecular viscosity of water. 1.5 x 10- m-s" Laplacian diffusion coefficient for mixing of K {T) heat horizontally. 1.5 x 10- m-s^ Laplacian diffusion coefficient for mixing of Kz(T) heat vertically. K {S) 1.5 x l(r m-sLaplacian diffusion coefficient for mixing of salt horizontally. K (S) 1.5 x 10- m-s-^ Laplacian diffusion coefficient for mixing of salt vertically. 0.52 sCoriolis parameter f 0 No /3-effect in this simulation 2xlO" KThermal expansion coefficient cue 7.4 x IO" Salinity contraction coefficient Ps hfacMin 0.2 States that the minimum grid thickness of a cell is 0.2Az Time step scheme Staggered Staggers the variables in time Advection scheme 3-DST with flux limiter Sets the advection scheme of the active tracers to a third order direct space-time with a flux limiting algorithm. EOS Linear Sets the equation of state to linear. 0 Neglects the non-hydrostatic terms. The opposite is true if it is set to 1. Number of time steps 1400 Number of time steps of the run, it is equivalent to 35 s. Ai 0.025 s Time step. Coordinate system Curvilinear Sets the coordinate system to a general curvilinear system instead of a Cartesian, or Polar spherical. Az 0.0028125 m Spacing between vertical levels. Parameter Value l(T m-s^ b h 9 2 h y y 2 h y Z 1 3 i 4 Other modifications to the code were necessary to allow the use of curvilinear coordinates instead of the default Cartesian coordinates. 2.3.2.6. Parameters. There are a large number of parameters, both physical and numerical, that are active for every run. However, only the parameters that are most significant for this experimental configuration, and are not default parameters are listed in Table 2. 2.3. E X P E R I M E N T CONFIGURATION 48 2.3.2.7. Numerical stability. The numerical stability of each term in the model can be analyzed, and although this does not guarantee the global stability of the model, it is a good guide for choosing the numerical parameters or for tracing errors. In his PhD thesis, Adcroft (1995) discusses the origin of the stability criterion relevant to models such as MITgcm. The following analysis is based on his work. To ensure that all wave numbers for the diffusion/dissipation terms decay without spurious oscillations, the following criterion should be satisfied ,„n, * £ < x where K represents the appropriate Laplacian diffusion coefficient or the eddy viscosity coefficient. For more than one dimension, each dimension contributes a term to the left hand side of the inequality. In this way, using the values of the parameters described above, the minimum horizontal and vertical spacing possible in the model, Arc = 0.16 cm, Az = 5.635 x 1 0 for the current run Ah = Kh(T) = Kh(S), and K (T) = K (S), V V - 2 cm , and noting that the stability criterion are 4n At = 0.0411 Ax h 2 An At = 1.9 x 10" Az AA At 0.316 Az v 5 2 v 2 It is seen that these are all below the stability limits. The C F L number (Section 2.2.1.2) for an extreme maximum horizontal velocity u = 3 c m s \u\At Ax - 1 is 0.48 This C F L number is very high by normal standards, but as it was explained in Section 2.2.4.2, when using the 3DST-flux limiting advection scheme the accuracy increases with the C F L number. It is seen in the results that this number does not introduce any instability in the model. For a non-hydrostatic run of the model, it is necessary to address the stability of the model for internal gravity waves. As an approximation to find the internal wave phase speed (c — a/k) to be used for the stability analysis, I use software by Dr. R. Pawlowicz, University of British Columbia, which solves the vertical mode equation for non-rotating linear internal waves. Due to the nature 49 2.3. E X P E R I M E N T C O N F I G U R A T I O N a) 1000 b) Buoyancy Frequency (s " ) 1 1500 2000 Wave number (m ) _1 F I G U R E 2.3.7. (a) Buoyancy frequency profile, (b) Dispersion relation for rotating internal waves in the laboratory tank. of the problem investigated here, rotation is likely to be important, so the software was modified to introduce the Coriolis frequency and the equation solved was: a w dz 2 (2.3.12) 2 N 2 a kw = 0 2 2 - / 2 where w is the vertical velocity, N is the buoyancy frequency, a is the internal wave frequency, k is the wave-number. The value of N (Figure 2.3.7a) is found using the initial density (Figure 2.3.5). To find a, several values between / and N are tested, choosing, for the stability analysis, the one that gives the highest phase speed. The dispersion relation is calculated (Figure 2.3.7b). When forcing a stratified fluid over topography, internal waves are more likely to be generated in locations where the slope of the bottom changes abruptly (this will be discussed with more detail in the next chapter). The fastest internal waves should occur at the depth of the abyssal plain (9 cm), where the continental slope starts (Ax = 0.0036 m), and the frequencies that give the highest phase speed should be those very close to / . For example, using a == 0.62 s , the maximum phase speed -1 obtained, corresponding to the first baroclinic mode, is 0.09 m - s . Then, the stability condition -1 for internal gravity waves is cAt Ax 0.64 therefore with the current time step the model should be stable.The shapes of the first three baroclinic modes are shown in Figure 2.3.8. Here I have assumed that internal waves are not generated near the inner wall of the tank. If this is the case the model should be unstable (the Courant 2.3. E X P E R I M E N T C O N F I G U R A T I O N 50 number would be 1.6 due to the smaller grid size). However I decided to take this risk in order to save computational time. In any case, if instabilities of this kind are detected, the time step just has to be reduced, and the problem would be solved. 2.3.2.8. High resolution experiment set up. Objective 2 in section 1.3 involved the study of the frictional effects in a submarine canyon. Perenne et al. (2001a) suggests that modeling studies with a resolution high enough to resolve the bottom boundary layer, would help to better determine the differences between numerical and laboratory models. In the current experiment the thickness of the bottom Ekman layer can be calculated as: = 0.196 cm with the initial resolution (Az = 0.28125 cm), this is approximately 70% of the cell thickness. Then, by doubling the vertical resolution, the boundary layer should be properly resolved. However, if just the vertical resolution is increased, there is a risk that the canyon shape would be altered. To avoid this, the resolution was doubled in every direction to 480 points in the azimuthal direction, 164 points in the radial direction, and 64 levels in the vertical direction. This gives a spacing at the shelf break of about 0.25 cm in the horizontal, and 0.14 cm in the vertical. Figures 2.3.9 and 2.3.10 depict the horizontal and vertical grids for the high resolution experiments. 2.3. E X P E R I M E N T C O N F I G U R A T I O N 51 0.34 0.32 0.3 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.1 0.15 0.2 0.25 0.3 FIGURE 2.3.9. Top view of the high resolution horizontal grid over the canyon. Distances are in meters. The bathymetric contours are are placed at 1 cm intervals. -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 -0.045 0.16 0.18 0.2 0.22 0.24 0.26 FIGURE 2.3.10. High resolution numerical grid of a vertical cross section of the submarine canyon in the laboratory tank. The bold line represents the real bathymetry. Scales in meters. A resolution this dense required the experiments to be run with multiple processors. This was done on the "Monster High-Performance Computing Super-Cluster" of the Geophysical Disaster Computational Fluid Dynamics Centre at the University of British Columbia. For the parallelization in MITgcm, the physical domain is decomposed into different regions. These regions are called tiles and each one of them is allocated to a logical processor. A processor can own more than one tile. For this experiment, the physical domain is decomposed into 2 0 tiles in the azimuthal direction, and one tile in the radial direction, using only one processor per tile. In 2.3. EXPERIMENT CONFIGURATION 52 this way, each tile has 24 x 164 x 64 points. Each tile shares three cells with the adjacent tile to ensure an accurate exchange of data structures (i.e. variable arrays) between tiles. For the parallelization of the float package some modifications were made because the model assumed by default that the domain is decomposed in both horizontal directions. However for this configuration the decomposition is only in the azimuthal direction and only in this direction the exchange of the floats between tiles is needed. All physical and numerical parameters are the same as for the lower resolution runs, with the exception of the time step, which was adjusted to At = 0.01 s (hence the number of time steps is increased to 3500), to obtain a similar numerical stability as in Section 2.3.2.7. CHAPTER 3 Results The results obtained with MITgcm are analyzed using figures that represent the horizontal flow vectors at different vertical levels, and vertical cross sections of model outputs (potential temperature, velocities, salinity, and at) across the mouth of the canyon, across the middle of the canyon, across the shelf upstream of the canyon, and along the canyon axis (Lines A , B , C, and D, respectively in Figure 3.0.1). A l l the figures are plotted at the end time of the run, 35 seconds after the initialization of the model. Particle tracks are used to compare the model outputs with the observations of the flow in the laboratory tank and with results obtained by Allen et al. (2003) from numerical simulations using the SCRUM model. LineD LineC - Una B I — \ 0 0? JO \> 0^ .7;;-— ^ Una A 0.03 0.04 o.05 o.06 0.07 0.08 ' • - - 20 40 60 80 100 120 140 160 180 200 220 240 3.0.1. Schematic view of the computational domain. The bold lines indicate the cross sections used to analyze model outputs. FIGURE 3.1. 3.1.1. V e l o c i t y field. Flow description The flow field is investigated initially using plots of the horizontal veloc- ity vectors (u and v) at different vertical levels. The azimuthal (u) and the radial components (v) of the horizontal velocity are not computed at exactly the same location due to the spatial staggering of the variables (Figure 2.1.1). Thus, for the vector plots, both variables are interpolated so that 53 3.1. F L O W DESCRIPTION 54 the middle point of every vector is located at the center of the cell, exactly on top of the active tracer points. The pattern of the resulting flow is in good agreement with observational studies (Hickey, 1997), and numerical studies (Allen, 1996; Klinck, 1996; Perenne et al. 2001b; Allen et al., 2003) of the circulation over submarine canyons. Near the surface the flow is not affected by the presence of the canyon (not shown). The influence of the canyon begins to be seen distinctively at around 1.2 cm depth, where the flow starts to turn inshore at the head of the canyon, and then deflects downstream (Figure 3.1.1). It is interesting to note that at the head of the canyon, inshore of the 2 cm isobath, the flow is almost completely radial. 0.34 0.32 0.3 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.1 0.15 0.2 0.25 0.3 3.1.1. Plot of the velocity vectors at level 5 (-1.27 cm). The solid lines represent the bathymetry in meters. The maximum velocity vector in this layer is 1.3 cm-s" . Note the deflection of the flow around the head of the canyon. Horizontal distances are in meters. FIGURE 1 The compression of water columns as they are upwelled onto the shelf causes negative vorticity downstream of the canyon and the return of the flow offshore. On the rest of the shelf, the flow is shoreward close to the bottom due to Ekman pumping. From about a depth of 1.5 cm downward a cyclonic eddy occurs near the head (Figure 3.1.2). The formation mechanism of this eddy, according to the theory presented in Section 1.2, is related to the increase in the stretching vorticity of the water columns as water falls from the shelf into the canyon. A t 2.39 cm, which is about the middle depth of the canyon (Figure 3.1.3), the eddy is stronger and quite elongated as it is basically following the bathymetric contours. At this level are 3.1. F L O W DESCRIPTION 55 observed the maximum cross shelf velocities (about 0.31 cm-s ) which occur near the upstream -1 and downstream walls. 0.1 0.15 0.25 0.2 0.3 F I G U R E 3.1.2. Plot of the velocity vectors at level 6 (-1.55 cm). The solid lines represent the bathymetry in meters. The maximum velocity vector in this layer is 1.36 cm-s" . A n asymmetric cyclonic eddy develops at the head of the canyon with strong flow on the downstream side. Horizontal distances are in meters. 1 0.1 0.15 0.25 0.2 0.3 F I G U R E 3.1.3. Plot of the velocity vectors at level 9 (-2.39 cm). The solid lines represent the bathymetry in meters. The maximum velocity vector in this layer is 1.44 cm-s outside the canyon, and 0.31 cm-s" at the canyon walls. At this depth the eddy basically follows the canyon walls. Horizontal distances are in meters. _1 1 Deeper in the canyon, the eddy is present (Figure 3.1.4), but becomes much weaker and more symmetric, until it eventually disappears at about a depth of 6.6 cm (not shown). 3.1. FLOW DESCRIPTION 56 However, a cyclonic eddy is not seen above the canyon rim as observations in the real ocean suggest (Hickey, 1997; Allen et al., 2001). This is not surprising since an eddy above the rim is not present in the laboratory tank either, due to scaling differences when compared to Astoria Canyon or Barkley Canyon. What is initially expected from the numerical model is to simulate the results in the tank, and eventually model more realistic oceanic conditions. 0.34 0.32 0.3 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.1 0.15 0.2 0.25 0.3 F I G U R E 3.1.4. Plot of the velocity vectors at level 12 (-3.23 cm). The solid lines represent the bathymetry in meters. The maximum velocity vector in this layer is 1.3 ciri'S" . The eddy is present but is more symmetric and the velocities are slower than higher in the water column. Horizontal distances are in meters. 1 3.1.2. Upstream conditions. A good way to test the solutions of the numerical model is to observe the conditions of the flow upstream of the canyon. In this region the bathymetry is simpler, and thus, the flow pattern should be simpler too. The dynamics in this region can be understood through Ekman layer theory (Ekman, 1905). When the water in the tank starts to rotate anticyclonically it causes a divergence in the bottom Ekman boundary layer. In this layer, the flow begins to rotate faster due to viscous stresses, creating a centrifugal force that is stronger than the pressure gradient, causing the radially outward Ekman transport. By continuity, this transport is fed with water from above the Ekman layer. The loss of water from the interior is compensated by a weak return flow toward the center of the tank. When the fluid is stratified, the Ekman transport causes dense fluid to be advected upslope into a region of lower density, thus inducing a gravitational downslope force that eventually can balance the rotational forces inhibiting further upslope transport. The timescale for the shutdown of the transport in the bottom boundary layer is T — S / [(1 + S ) S f cos 2 s s (MacCready and Rhines, 3.1. F L O W DESCRIPTION 57 1991), where S = ( i V s i n 7 / / c o s 7 ) is a slope Burger number, and 7 is the bottom angle. For 2 s the current simulations the shutdown timescale T is 4 s over the continental shelf (f—0.52 s , _1 N=3 s , 7 = 6 ° ) , and 6 . 4 x l G _1 3 over the continental slope (/=0.52 s , iV=1.8 s , 7=45°). This _1 s timescale reduce the fiux out of the boundary layer by a factor of ( T / t ) _1 1//2 (MacCready and Rhines, 1991), where t is the spin-up time. In the present case the time of the simulations is 35 s, therefore it is expected that the flux in the Ekman layer will be reduced around 70% over the shelf, and practically non-existent over the slope as a consequence of the stratification. For a nice analysis of the dynamics and evolution of homogeneous and stratified flow in a laboratory tank, with and without a submarine canyon, see Mirshak (2001). Vertical cross sections of the velocity components (u, v, w), potential temperature (#), and a t are used to analyze the regime upstream of the canyon. The contours are plotted across the shelf along Line C (Figure 3.0.1). Figure 3.1.5, shows the azimuthal velocity (u) with positive velocities coming out of the page. Offshore, in the deepest part of the tank, the velocity increases smoothly toward the shelf with almost no vertical shear until the near the bottom. The maximum velocity of 1.42 cm-s" occurs 1 at the shelf break. Over the shelf the vertical shear increases and the velocity starts to decrease toward the outer wall of the tank due to side and bottom friction. The wiggles in the velocity near the bottom are related to the fact that the values of a given variable in smaller cells (the ones in direct contact with the topography) are calculated closer to the bottom than the values calculated for full cells. This is can be observed in Figure 3.1.5, where the azimuthal velocity across the shelf is plotted in each cell without any interpolation. Figure 3.1.6, shows the radial velocity (v) across the shelf. Positive velocities are toward the outer wall of the tank. At the bottom over the abyssal plain and over the shelf, the water is pushed toward the outer wall due to Ekman dynamics, then, by continuity, the surface flow is directed toward the center of the tank. Figure 3.1.7, shows the vertical velocity across the shelf. There are more small scale features in this plot, but they coincide with features observed in the other velocity fields (especially in the v plot), as the vertical velocity is obtained diagnostically from the integration of the continuity equation. Ignoring the smaller features at about mid-depth above the shelf, the rest of the pattern shows the expected downward flow at the center that is supplying the water for the Ekman transport in the boundary layers. Also by continuity, water is upwelled over the slope and shelf. 3.1. F L O W DESCRIPTION 15 0.2 0.25 0.3 0.35 0.4 0.45 58 0.5 iniiiiii II iiHiiiiiiliiiiiNiiiiiijii .•is'is.JMss . .llllllitii.. Hllillli ttUi • i \-0.09' b) 15 f 02 0 25 03 0 35 0.4 0 45 1.4215 1.3504 1.2793 1.2082 1.1372 1.0661 0.995 0.924 0.8529 0.7818 0.7107 0.6397 0.5686 0.4975 04264 0.3554 0.2643 0.2132 0.1421 0.0711 05 3.1.5. Vertical cross section of the along shore velocity across the shelf (Line C). (a) Contour plot of u. Note the wiggles over the shelf, (b) Plot using each value of u in its corresponding cell. The thickness of each cell is represented by the relevant non-dimensional "h-factor" times the vertical spacing, that is, to plot a cross section of u, the thickness of each cell is h • Az; for v plots, the thickness is h • Az; and for w, and 6 plots the thickness is h • Az. It can be observed that the changes near the bottom coincide with the changes of vertical level. Positive velocities coming out of the page. Velocity is in cm-s . Depth and horizontal distance are in meters. FIGURE w a c -1 In the model, (3.1.1) (2.3.8) takes the form: p' = p 0 (ft (S(z) - S ) - ag (e(z) - 6 )) ref ref where S(z) is the salinity calculated by the model, S f re — 35 is the reference salinity, 9(z) is the potential temperature calculated by the model, and 6 f is the reference temperature obtained from re (2.3.10) and plotted in Figure 2.3.6. As the salinity is constant, the first term inside the brackets 3.1. FLOW DESCRIPTION 3.1.6. Vertical cross section of the radial velocity across the shelf (Line C). Positive (yellow and red) velocities are to the right. Velocity is in cm-s . Depth and horizontal distances are in meters. FIGURE -1 in (3.1.1) is zero, and the equation of state becomes: p' = -p Ctg 0 (6(z) - 0 ) ref then, the total density is (3.1.2) P = P' + Pref 3.1. FLOW DESCRIPTION where p j re 60 is the initial density profile from the tank (Figure 2.3.5). For convenience, in Figure 3.1.8 the quantity that is presented is ,n , ON (P ~ 3.1.3 a = t 1 0 0 0 k S • m " ) 3 : kg • m 13 It can be observed that the inclination of the isopycnals is toward the outer wall of the tank, which corresponds with the expected upwelling. The steeper isopycnals are located at the point where the continental slope meets the abyssal plain, and at the shelf break, that is, where the most pronounced changes in the bathymetry occur. FIGURE 3.1.8. Vertical cross section of the density across the shelf (Line C). Units are a . Depth and horizontal distances are in meters. t 3.1.3. Flow at the canyon mouth (Line A ) . Figure 3.1.9 presents a vertical cross section of the azimuthal velocity (u) at the mouth of the canyon (Line A). The view is toward the head of the canyon. The predominant feature of this plot is the strong vertical shear between the flow inside the canyon and the shelf. In the interior of the canyon, the negative u velocities represent an interior cyclonic eddy. The flow on the downstream side is stronger near the canyon rim in the negative direction. The flow near the bottom is almost zero. This plot is very similar to the plot of Allen et al. (2003) using the model SCRUM (Figure 3.1.10). The radial velocity plot at the mouth of the canyon (Figure 3.1.11) shows some negative (offshore) velocities at the upstream rim of the canyon, which are not present in the SCRUM results (Figure 3.1.12). Inside the canyon, the negative and positive radial velocities near the upstream and downstream walls, correspond to the interior eddy and are very similar in both SCRUM and 3.1. F L O W DESCRIPTION 61 3.1.9. Vertical cross section of the azimuthal velocity at the canyon mouth (Line A). Positive velocities are to the right. The view is toward the head of the canyon. Velocities are calculated at the middle of each cell face. Velocity is in cm-s . Depth and horizontal distances are in meters. Note the strong vertical shear. FIGURE -1 F I G U R E 3.1.10. Vertical cross section of the along shore velocity at the canyon mouth from Allen et al. (2003) using the model S C R U M . Positive velocities are to the right. The view is toward the head of the canyon. The vertical scale is 5 cm from top to bottom. Velocity units are c m s . Reproduced from Allen et al. (2003) with permission of Dr. S.E. Allen. A . - 1 MITgcm. The positive velocity layer at the bottom upstream and downstream of the canyon is the Ekman transport mentioned above. Note that at the upstream corner of the canyon the Ekman transport works against the cyclonic eddy, and can inhibit to some extent the formation of an eddy above the canyon rim, while on the downstream side over the shelf, both the Ekman transport 3.1. F L O W DESCRIPTION 62 and the flow from the cyclonic eddy are directed inshore, which adds to the already high positive velocities in this region. F I G U R E 3.1.11. Vertical cross section of the radial velocity at the mouth of the canyon (Line A). Positive velocities are into the page. The view is toward the head of the canyon. Velocity is in cm-s . Depth and horizontal distances are in meters. -1 -0.M FIGURE 3.1.12. A cross section of the radial velocity at the mouth of the canyon from Allen et al. ( 2 0 0 3 ) using the SCRUM model. Positive velocities are to the right. The view is is toward the head of the canyon. The vertical scale is 5 cm from top to bottom. Velocity units are cm-s . Reproduced from Allen et al. ( 2 0 0 3 ) with permission of Dr. S.E. Allen. A . -1 Figure 3.1.13 shows the vertical cross section of vertical velocity at the mouth of the canyon. There are many small scale features inside and outside of the canyon. However, the results from MITgcm show clear downward (negative) velocities at the upstream rim of the canyon, corresponding to the water falling from the shelf, and upward velocities in the downstream rim, corresponding to 3.1. F L O W DESCRIPTION 63 the upwelling of water onto the shelf. This is the opposite of what happens in S C R U M (Figure 3.1.14), where upward and downward velocities are contrary to the current understanding of canyon dynamics. As explained in Section 1.2.2.1, the reason for these inconsistencies is related to the errors in the vertical advection of density in a terrain-following model (Allen et al., 2003). 0.2 0.18 0.22 0.24 F I G U R E 3.1.13. Vertical cross section of the vertical velocity at the mouth of the canyon (Line A ) . Positive velocities are upward. The view is toward the head of the canyon. Velocity units in cm-s . Depth and horizontal distances are in meters. -1 X l O " ' F I G U R E 3.1.14. Vertical cross section of the vertical velocity at at the mouth of the canyon from Allen et al. (2003). Positive velocities are upward. The view is toward the head of the canyon. The vertical scale is 5 cm from top to bottom. Velocity units in c m s . Reproduced from Allen et al. (2003) with permission of Dr. S.E. Allen. A . - 1 In Figure 3.1.15, er is contoured at the mouth of the canyon (Line A ) . Although this figure is t similar to the one in Allen et al.(2003) (Figure 3.1.16), in the results from MITgcm, the isopycnals 3.1. F L O W DESCRIPTION 64 over the upstream rim of the canyon (top of yellow band) tilt slightly downward suggesting local stretching of the water columns above. At rim depth, the isopycnals are compressed (yellow and orange bands). Downstream the stretching is stronger as isopycnals tilt upward over the rim (yellow and green bands) and downward at rim depth (orange band). Downstream over the shelf the isopycnals are compressed near the bottom. Between about 2.5 and 3.5 cm depth inside the canyon, the isopycnals are stretched when compared to the initial values (not shown), with stronger stretching at the upstream wall. Between 3.5 and the bottom the isopycnals are compressed compared to the initial values. 0.16 0.16 0.2 0.22 0.24 0.26 Vertical cross section of the density at the mouth of the canyon (Line A ) . Units are at- The isopycnals over the canyon (yellow and green bands) bend downward but not as strongly as field observations suggest (Hickey, 1997). F I G U R E 3.1.15. 3.1.4. Flow at the middle of the canyon (Line B). In Figure 3.1.17 the contours of the azimuthal velocity show basically the same features as those at the mouth of the canyon (Figure 3.1.9): strong vertical shear at rim depth, and negative velocities inside the canyon. The differences are that the magnitudes of the negative velocities are reduced. For example, while the maximum negative azimuthal velocity at the mouth of the canyon is - 0 . 1 8 cm-s , in the middle of the canyon -1 it is -0.088 cm-s . Also, the positive velocities penetrate further down below the rim at the mouth -1 of the canyon. The radial velocity across the middle of the canyon is depicted in Figure 3.1.18. Offshore velocities are not present at the upstream rim of the canyon, and the interior of the canyon is practically isolated by a layer of positive radial velocity. The interior eddy is a predominant feature of the flow in this area, but is not as strong as at the mouth of the canyon. 3.1. F L O W DESCRIPTION 65 F I G U R E 3.1.16. A cross section of the density at the mouth of the canyon (Line A) from Allen et al. (2003). The vertical scale is 5 cm from top to bottom. Units are at- Very small isopycnal deflections are seen compared to field observations (Hickey, 1997). In particular, over the canyon near the upstream rim (yellow band) isopycnals tilt upward to the right in a region observed to have strongly descending isopycnals. Reproduced from Allen et al. (2003) with permission of Dr. S.E. Allen. A . 1.1891 1.1252 1 0613 0.9975 0.9336 0.8697 08059 0.742 0.6781 0.6143 0 5504 0.4865 0 4227 0.3588 0.2949 0.2311 0.1672 0.1033 0.0395 -0.0244 -0.0883 F I G U R E 3.1.17. Vertical cross section of the azimuthal velocity at approximately the middle of the canyon (Line B). Positive velocities are to the right. The view is toward the head of the canyon. Velocity units are in cm-s . Depth and horizontal distances are in meters. -1 The vertical velocity across the middle of the canyon has positive values at rim depth (Figure 3.1.19). Negative velocities at the upstream rim do not penetrate into the canyon, suggesting that there is no water from the shelf that is falling into the canyon. The rest of the features are similar to those at the mouth of the canyon (Figure 3.1.13) where the velocity is downward at mid-depth inside the canyon, and upward at the bottom. 3.1. F L O W DESCRIPTION 66 FIGURE 3 . 1 . 1 8 . Vertical cross section of the radial velocity at approximately the middle of the canyon (Line B). Positive velocities are into the page. The view is toward the head of the canyon. Velocity units are in cm-s . Depth and horizontal distances are in meters. -1 0.18 0.19 0.2 0.21 0 23 0.22 0 24 0.25 0.26 FIGURE 3 . 1 . 1 9 . Vertical cross section of the vertical velocity at approximately the middle of the canyon (Line B). Positive velocities are upward. The view is toward the head of the canyon. Velocity is in cm-s . Depth and horizontal distances are in meters. -1 In Figure 3.1.20 are shown the contours of a across the middle of the canyon. Although it is very t similar to Figure 3.1.15, there are differences. The stretching of the isopycnals at the downstream rim compared with the initial values (not shown) is greater here than at the mouth of the canyon, however, at the upstream rim the curvature of the isopycnals is not as noticeable. Inside the canyon 3.1. F L O W D E S C R I P T I O N 67 the stretching of the isopycnals at mid-depth and the compression at the bottom compared with the initial values are also more marked than at the mouth of the canyon. 0.21 0.22 0.23 0.24 0.25 FIGURE 3.1.20. Vertical cross section of the density at approximately the middle of the canyon (Line B ) . The view is toward the head of the canyon. Units are in a . Depth and horizontal distances are in meters. t 3.1.5. Flow along canyon axis (Line D). The interior of the canyon has small negative azimuthal velocities (Figure 3.1.21). The maximum positive velocity (1.47 c m - s ) occurs at the - 1 position were the shelf break outside of the canyon would be if there was no canyon. 0.15 0.2 0.25 0.3 FIGURE 3.1.21. Vertical cross section of the azimuthal velocity along the canyon axis (Line D ) . The bold dashed line represents the position of the shelf where there is no canyon. Positive velocities are out of the page. Velocity is i n c m - s . Depth and horizontal distances are in meters. - 1 3.1. FLOW DESCRIPTION 68 The radial velocity along the canyon axis (Figure 3.1.22), shows more distinct features. The positive velocities are stronger than upstream (compare with Figure 3.1.6). Also, in the region outside of the canyon (outside of the dashed line), the positive velocities occupy a much broader area and stretch almost continuously from the abyssal plain to the coast (right hand side of the plot), reaching a maximum velocity of 0.2 cm-s at the head of the canyon. Inside the canyon -1 most of the flow is off-shore with positive velocities found only at the bottom, at the bottom at mid-canyon and just below the rim between the shelf break and the head of the canyon. The surface flow is off-shore, just as in the upstream plots. 3.1.22. Vertical cross section of the radial velocity along the canyon axis (Line D). The dashed line represents the position of the shelf where there is no canyon. Positive velocities are to the right. Velocity is in cm-s . Depth and horizontal distances are in meters. FIGURE -1 The vertical velocity contours (Figure 3.1.23), show many small scale features; however as in most of the plots the correlation of the vertical velocity with the radial velocity is easy to spot. Positive velocities are located above the canyon rim, inside the canyon near the bottom at midcanyon, and outside the canyon from mid-depth toward the bottom, which roughly coincides with the positive radial velocities discussed in the previous paragraph. Negative velocities are strongest just below the rim near the head of the canyon, and at the position where the shelf break would be if there was no canyon. The isopycnal distribution along the canyon axis is depicted in Figure 3.1.24. Noticeable is the strong tilt in the isopycnals toward the head of the canyon at rim depth, which suggests enhancement of upwelling in this area. Below the rim at approximately 2.5 cm depth, the isopycnals bend 3.1. FLOW DESCRIPTION 69 3.1.23. Vertical cross section of the vertical velocity along the canyon axis (Line D). The bold dashed line represents the position of the shelf where there is no canyon. Positive velocities are upward. Velocity is in cm-s . Depth and horizontal distances are in meters. FIGURE -1 downward as it would be expected from the negative vertical velocities in this area. Below the shelf break depth, the isopycnals tilt upward where the slope would be, and once inside the canyon they flatten. The distribution outside the canyon is very similar to the distribution upstream shown previously (Figure 3.1.8). 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 FIGURE 3.1.24. Vertical cross section of the density along the canyon axis (Line D). The dashed line represents the position of the shelf where there is no canyon. Units are in a - There is a strong lifting of the isopycnals toward the head of the canyon due to enhanced upwelling. Depth and horizontal distances are in meters. t 3.2. COMPARISON W I T H L A B O R A T O R Y RESULTS 70 3.2. Comparison with laboratory results As described in Section 2.3.2.5, to compare the model and the laboratory results, float particles were introduced into the model at the same horizontal position where particles are observed to start to move in the tank. The vertical position of the particles in the tank can only be estimated within 1 cm depth, by the color that they reflect (i.e. green from 1.2 to 2.2 cm depth or red 2.2 to 3.2 cm depth). It is expected that the particles move differently within the 1 cm increments due to the significant vertical shear (Figure 3.1.9). For this reason, in the numerical model particles were planted every millimeter in depth for each horizontal position. Then, to the 32 particles injected in the tank, correspond 352 particles injected in the numerical model. The tracers are planted 30 s after the start of the simulation and the positions were captured every 0.5 s during 5 s. From the 11 numerical tracers corresponding to each tracer in the laboratory tank, the one that best matched the final position of the real tracer was chosen for this comparison. Figure 3.2.1 shows the tracks from the laboratory tank (empty red squares), from MITgcm (filled blue squares), and from S C R U M (empty black squares). The positions are plotted every 0.5 s. Near the surface (Figure 3.2.1A) over the upstream rim of the canyon, the laboratory tracks show a stronger off-shore velocity component than either S C R U M or MITgcm. The numerical models show a systematic deviation compared to the laboratory model. Near the head of the canyon MITgcm does a better job than S C R U M simulating the tank, and although it is not appreciable from the plot, the rest of the surface tracks on the upstream side are closer from MITgcm than from S C R U M . Downstream of the canyon, both models underestimate the along shore velocity component. Between 1.4 and 1.63 cm depth (Figure 3.2.1B) and between 1.8 and 2.2 cm (Figure 3.2.1C) both models reproduce the laboratory conditions quite well, especially at the downstream side of the canyon, where the upwelling is clear. The numerical tracks at the upstream rim of the canyon show again the systematic difference with the laboratory tracks; however, the results from MITgcm are much closer than the results from SCRUM. This agrees with the representation of the radial velocity across the mouth of the canyon discussed above (Figures 3.1.11, and 3.1.12). In the deep flow, the interior cyclonic eddy is present in both the laboratory and the numerical models, but is shorter and narrower in the models (Figure 3.2.ID). At the mouth of the canyon, the flow seems to be slower in the laboratory than in both models. Again, the downstream tracks are are the ones that match better. 3.2. COMPARISON W I T H L A B O R A T O R Y RESULTS 71 0.1 0.15 0.2 0.25 0.3 0.1 0.15 0.2 0.25 0.3 0.1 0,15 0.2 0.25 0.3 0.1 0.15 0.2 0.25 0.3 F I G U R E 3.2.1. Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm (filled blue squares), and S C R U M (unfilled black squares). (A) tracks between 1.2 and 1.3 cm depth, (B) tracks between 1.4 and 1.63 cm depth, (C) tracks between 1.8 and 2.2 cm depth, (D) tracks between 2.2 and 3.2 cm depth. Lines are bathymetric contours. Distances and bathymetric depths are in m. A quantitative comparison was made using the same method as Allen et al. (2003). A normalized error is calculated using the expression: J(Xlaboratory — Xmodel) - (Ylaboratory Error = — LabFluctuationV elocity * time 2 Ymodel) 2 where the laboratory fluctuation velocity is the error in the velocity between runs, estimated as 0.2 cms - 1 (Allen et al., 2003), and the time is the time lapse since a particle started moving. The results from this error calculation are shown in Table 3. The particle with the biggest error is the one that is near the downstream rim at the mouth of the canyon in Figure 3.2.ID. This particle is located very close to a point where the flow that goes inside the canyon separates from the flow that moves along the slope (Figure 3.1.3, and Figure 3.1.4); thus a very small difference in the horizontal position of the float would make the float change 3.2. COMPARISON WITH L A B O R A T O R Y RESULTS 72 T A B L E 3. Initial and final depths of floats in meters (numerical) and errors between MITgcm and laboratory simulations. Float number Initial Depth Final Depth Aver Error End Error 27 1 8 11 12 18 13 23 25 7 26 17 30 2 28 33 21 22 19 24 29 6 20 4 3 15 14 5 31 9 10 16 -0.0320 -0.0320 -0.0320 -0.0310 -0.0270 -0.0240 -0.0240 -0.0220 -0.0220 -0.0220 -0.0220 -0.0205 -0.0200 -0.0200 -0.0190 -0.0186 -0.0180 -0.0180 -0.0160 -0.0160 -0.0160 -0.0154 -0.0152 -0.0150 -0.0150 -0.0140 -0.0120 -0.0120 -0.0120 -0.0120 -0.0120 -0.0120 -0.0310 -0.0320 -0.0320 -0.0309 -0.0272 -0.0241 -0.0240 -0.0217 -0.0220 -0.0220 -0.0208 -0.0194 -0.0200 -0.0208 -0.0193 -0.0177 -0.0167 -0.0174 -0.0157 -0.0151 -0.0162 -0.0154 -0.0147 -0.0152 -0.0151 -0.0133 -0.0116 -0.0121 -0.0120 -0.0120 -0.0120 -0.0120 0.6357 1.2420 2.7112 0.5772 0.8564 1.1120 1.0903 0.3430 5.1769 0.6153 0.5873 1.1902 0.4012 0.3957 0.4471 0.2460 0.2839 0.6791 0.3456 0.2934 0.6663 0.4183 0.4062 0.7507 0.1894 0.4144 0.9043 0.4428 2.6461 0.8290 0.6264 0.3152 0.0851 1.0394 2.1379 0.5064 0.7505 0.9294 0.9309 0.6572 6.0221 0.7456 0.4283 0.0399 0.1993 0.3270 0.4430 0.0499 0.1686 0.8583 0.1486 0.4544 0.1196 0.0624 0.0241 0.5979 0.1768 0.2890 0.7656 0.4287 2.4005 0.6049 0.5382 0.2044 paths. For this reason this particle was dismissed from the calculation of the average error of the final position. This error is 0.55 in MITgcm, and 0.65 in S C R U M . In this run, MITgcm has been configured as closely as possible to the configuration of S C R U M used by Allen et al. (2003); that is, both models use the same values for the physical parameters, both use the hydrostatic approximation, and have a similar vertical resolution inside the canyon. However, the horizontal scales are different. Allen et al. (2003) increased the horizontal scales in the numerical model by ten-fold, because the slopes are too steep for S C R U M to model accurately the 3.2. COMPARISON WITH LABORATORY RESULTS 73 horizontal pressure gradients. To make sure that this scaling did not introduced additional errors, a test run using MITgcm was done in which the horizontal scales were increased by ten-fold, and the vertical Laplacian diffusion coefficients and vertical eddy viscosities were reduced by one order of magnitude, just like in Allen et al. (2003). The results of this modified run (not shown here) introduced no significant differences. For this reason, it can be assumed that the main difference between MITgcm and SCRUM is the representation of the vertical coordinate. So far, the improvements seen from the use of MITgcm to model the flow in the laboratory tank are: a more realistic representation of the vertical velocity around the canyon and the presence of negative radial velocities at the upstream rim; both of which contribute to reduce the error in the final positions of the tracer particles. These improvements can be regarded as a consequence of the use of a z-coordinate model over a terrain-following coordinate model, and demonstrate that the errors in the vertical advection of density in SCRUM, observed by Allen et al. (2003), are, at least, part of the problem that leads to the differences between laboratory and numerical results. However, although errors are diminished, there are still differences between the laboratory and the model results. Summarizing, the main differences include consistent underestimation of negative radial velocities near the upstream wall of the canyon, and the wider and broader eddy observed in the laboratory. The main concern for this work are the differences at the upstream rim, as the consistency of these differences indicate that a basic feature of the flow is being missed by the models. The differences in the interior of the canyon are treated as secondary problems because in this region the flow is changing within very short horizontal and vertical scales, and any small error in the position of a particle would result in a very different pattern. In addition to errors in the vertical advection of density due to the use of terrain-following models (Allen et al., 2003), Perenne et al. (2001a) suggest that differences between laboratory and model results near the rim of a canyon could be due to a models misrepresentation of the bottom boundary layer. This is consistent with the errors found in MITgcm and discussed in Section 2.3.2.3. To analyze the influence of the bottom drag in the model results, various tests using different values for the bottom drag are described next. 3.2.1. Sensitivity to changes in the bottom drag. For the numerical model tracks to agree with the laboratory model at the upstream side of the canyon, the numerical flow must have high azimuthal velocities (at least 1 cm-s ) and high off-shore velocities. In the run with the original -1 drag parameterization, this is not achieved. The azimuthal velocity of the particles is not high 3.2. COMPARISON WITH L A B O R A T O R Y RESULTS 74 enough at the depth where the off-shore velocities are observed at the upstream side (Figure 3.1.9, and Figure 3.1.11 at about 1.5 cm), and the off-shore velocities are not high enough at shallower depths. As an artificial means to increase the velocities over the shelf, for this test, the bottom drag was arbitrarily reduced to one fourth of the original drag, that is, (2.3.5) becomes: and similarly with (2.3.6). In this new run, due to the reduction in the bottom drag, the match between the numerical tracks and the laboratory tracks at the upstream rim of the canyon is improved (Figure 3.2.2, note that the unfilled black squares correspond now to the full drag run). The average error in the final positions of the tracers located in this area goes from 0.48 to 0.35 (around 30% improvement). The numerical particles at the upstream side of the canyon that best match the laboratory particles for this run are at greater depth than those in the full drag run as a result of the increase in the azimuthal velocity. However, inside the canyon the numerical model tracks are not closer to the laboratory tracks and in many cases are farther apart. The tracks of those particles away from the shelf are basically unchanged, as the flow is not in contact with topography. The average error for the final position of the tracers is 0.59 in this run. To better understand the changes in the flow regime induced by the reduction of the bottom drag, the difference between the radial velocity is compared (Figure 3.2.3). Figure 3.2.3 (a) shows contours of the radial velocity at the mouth of the canyon (Line A) for the run with reduced drag. The main features of the plot are very similar to Figure 3.1.11. However, there are important differences at the upstream side of the canyon. In the run with full drag, there are positive velocities in most of this area, especially in the layer in contact with the bottom (which is attributed to the Ekman transport in the bottom boundary layer), while in the run with the bottom drag reduced the velocities in this area are in the opposite direction. The differences between the two runs can be easily observed in Figure 3.2.3 (b), where the contours show the result of subtracting the the reduced drag from the original run. The highest anomalies are located at the upstream side of the canyon over the shelf, where the velocities in the reduced drag run are more negative, and in the downstream interior of the canyon, where the velocities are much more positive in the reduced drag run. 3.2. COMPARISON W I T H L A B O R A T O R Y RESULTS 75 FIGURE 3.2.2. Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm reduced drag run (filled blue squares), and MITgcm full drag run (unfilled black squares), (a) tracks between 1.2 and 1.5 cm depth, (b) tracks between 1.6 and 1.7 cm depth. Lines are bathymetric contours. Distances and bathymetric depths are in m. These results can be predicted using classic Ekman layer theory. In the boundary layer, the time derivative in the momentum equations (2.1.1), (2.1.2), and (2.1.3) can be neglected by assuming that the boundary layer flow remains in a quasi-steady state with the interior flow, which evolves on a longer time scale. The balance of the horizontal momentum equation is between the viscosity and the Coriolis terms, that is: (3.2.1) fkxu = d . du —A — v 3.2. C O M P A R I S O N W I T H L A B O R A T O R Y R E S U L T S 76 FIGURE 3.2.3. (a) Vertical cross section of the radial velocity at the mouth of the canyon (Line A ) with the drag reduced minus the original value. Positive velocities are into the page, (b) radial velocity with drag reduced minus original radial velocity. The view is toward the head of the canyon. Velocity is i n c m - s . Depth and horizontal distances are i n meters. - 1 where k is a unit vertical vector, A is the vertical eddy diffusivity, which i n our case is the molecular v viscosity of water ( 1 0 - 6 m - s ) , and u — (u - U ) i + vj is a departure, induced by the presence of - 2 g a rigid bottom and friction, from the interior geostrophic flow pdy 9 Using the fact that the stress exerted by the rigid surface on the fluid is f — A du/dz v , then (3.2.1) becomes: fkxu = Assuming the stress f is zero outside the boundary layer, then, integrating vertically across the layer we get fkxf E =f where TE is the transport of the boundary layer or the E k m a n transport. Thus, the E k m a n transport can be expressed as: (3.2.2) f E = ^ * Now, when the bottom drag is reduced, the E k m a n transport i n the bottom boundary layer is reduced as well. Over the shelf at the upstream side of the canyon this transport is directed inshore. Due to the reduction of the transport i n the bottom boundary layer, the off-shore flow generated 3.2. COMPARISON WITH L A B O R A T O R Y RESULTS 77 FIGURE 3.2.4. Vertical cross section of the density at the mouth of the canyon (Line A) for the reduced drag run. Units are at- The isopycnals at rim level inside the canyon (yellow and orange bands) bend downward but much more similar to field observations (Hickey, 1997) than for the full drag run. by the vorticity inside the canyon is more dominant. This explains why the bottom Ekman layer is no longer visible at the upstream side. Downstream, as the bottom drag is reduced, the up-canyon flow that comes from the rupture of the geostrophic balance inside the canyon finds less resistance from the bottom drag, and hence, the positive radial velocities are increased as well. The density distribution across the mouth of the canyon (Figure 3.2.4) shows greater similarities to the observations across the mouth of Astoria Canyon (Figure 1.2.1) than the run using the full drag, especially for the isopycnals close to the rim of the canyon (yellow and orange bands). This alone suggest that the frictional effects are too strong in the numerical model, and could be overwhelming other effects inside the canyon such as stretching of the isopycnals. Other runs with the bottom drag further reduced show basically the same features, exaggerated as the configuration approaches that of free slip at the bottom. As an example of the effects that the free slip condition has on the resultant flow, a plot with some distinctive tracers is presented in Figure 3.2.5. Note that the matching of the final positions at the upstream side is almost perfect, while on the downstream side has worsened considerably. Perenne et al. (2001b), conducted tests with different values of bottom drag coefficients and interior mixing (A ), obtaining similar results as the ones shown here. v They "calibrate" their model for subsequent tests using the values of bottom drag and interior mixing that offer the best agreement between laboratory and numerical results. Although in the present work this could also 3.2. COMPARISON W I T H L A B O R A T O R Y RESULTS 78 0.34, 0.32 0.3 0.28 0.26 0.24 0.22 0.2 0.18 0.16 0.1 0.15 0.25 0.2 0.3 FIGURE 3.2.5. Particle tracks of floats in laboratory experiment (unfilled red squares), MITgcm free slip run (filled blue squares), and MITgcm full drag run (unfilled black squares). The tracks are selected to show the most characteristic features of the flow with a free slip bottom boundary condition. Lines are bathymetric contours. Distances and bathymetric depths are in m. be achieved, this approach is not used, because modifying the values for interior mixing or bottom drag in this way is not realistic, and it can conceal other important errors that affect the agreement of laboratory and model results. It is also important to recall that having a carefully tuned model for the laboratory is not the goal of our research, as what works well for the laboratory may not be applicable for the real ocean. For this reason, the simulation should be as realistic as possible within the limits of a finite-difference numerical model. Although these tests of reducing the bottom drag do not provide an explanation for the differences between laboratory and numerical results, they are successful in showing that changes in the bottom drag could seriously affect the pattern of the resultant flow in and around a submarine canyon. In the real ocean, turbulence can be regarded as a form of drag. Many studies suggest that in submarine canyons, turbulence mixing is enhanced within the bottom boundary layer (e.g. Hotchkiss and Wunsch, 1982; Lueck and Osborn, 1985; Kunze et al., 2002), mainly due to focusing of internal waves. The highest internal-wave energy levels observed are related to internal tides (Gardner, 1989; Petruncio et al., 1998), therefore one could imagine that the flow in a submarine canyon could change significantly between spring and neap tides, due to the changes in the thickness of the bottom boundary layer (i.e. changes in the bottom drag). Kunze et al. (2002) also point out that with eddy diffusivities on the order of 1 0 0 x l 0 m s , water can communicate diffusively _ 4 2 _ 1 over 100 m in the vertical in only ten days (compared to a year or more in the open ocean). This 3.3. HIGH RESOLUTION RESULTS 79 time scale is comparable to the frequency of upwelling events. However not many studies to date have addressed the influence that the changes in the thickness of the turbulent bottom boundary layer would have on the flow pattern near a submarine canyon during an upwelling event. Now that the bottom drag is recognized as an important factor in determining the pattern of the flow in a submarine canyon, it is important to understand whether the numerical discretization of the bottom boundary layer could be introducing additional errors in the numerical model. In the following section, the results from a run with increased resolution are presented in order to investigate this question. 3.3. H i g h resolution results As discussed in section 2.3.2.8, the bottom Ekman layer (~0.196 cm) occupies about 70% of the bottom cells (~0.28 cm) in the low resolution runs discussed above. For this new run the spatial resolution is increased by a factor two in every direction, with a resultant vertical cell thickness of ~0.14 cm. With this resolution, the flow in most of the Ekman layer is explicitly resolved, and any numerical discretization is only needed for a very thin portion of the depth (a maximum of 0.07 cm from the bottom), constraining any errors in the parameterization to this layer. In Figure 3.3.1 are summarized the results from this run upstream across the shelf (Line C). The main characteristics of the flow remain unchanged, and, as it could be expected, the solutions are much smoother than in the coarser resolution run. The differences are subtle yet observable, for example in Figure 3.3.1(a), the wiggles observed in the coarser resolution run have practically disappeared, because the distance between cell centers (where the velocity is calculated) has been reduced. In Figure 3.3.1(b), the radial velocity contours show smaller velocities than in the coarse resolution run, also it can be seen that the thickness of the bottom boundary layer over the shelf has been reduced as a result of the increased vertical resolution. This is consistent with the vertical velocities depicted in Figure 3.3.1(c). The displacement of the isopycnals also shows more upwelling in this case (compare Figure 3.3.1(d) with 3.1.8) . The flow along the canyon axis (not shown) behaves in a similar manner. Figure 3.3.2 depicts the flow across the canyon mouth (Line A), and Figure 3.3.3 does the same across the middle of the canyon (Line B). Both show very similar features as the flow in the coarser resolution run. In Figure 3.3.2(a), the azimuthal velocities are basically the same, and exhibit the same features that those in the coarser resolution run (Figure 3.1.9). Also in Figure 3.3.2(b), 3.3. HIGH RESOLUTION RESULTS 80 FIGURE 3.3.1. Results from the high resolution run (480x164x64 points) upstream from the canyon along Line C. (a) Azimuthal velocity, (b) Radial velocity, (c) Vertical velocity, (d) Density. Velocity units are in m - s , and density is in a . Depth and horizontal distances are in meters. -1 t the radial velocity contours show some off-shore velocities near the upstream rim of the canyon with roughly the same magnitude as those in the coarse resolution run (Figure 3.1.11). As it was explained above, these velocities are not strong enough nor shallow enough to modify the tracks of the particles (not shown). What is needed is either that the azimuthal velocity becomes stronger near the bottom (about 0.9 c m s , which is highly unlikely), or that negative radial velocities of - 1 0.1 cm-s -1 reach higher up in the water column (up to about 1.2 cm depth). The vertical velocity in Figure 3.3.2(c), shows that negative velocities are much more confined to the upstream wall of the canyon. In addition, the two patches of positive vertical velocity observed in this region in the coarse run (Figure 3.1.13) are not present in the higher resolution run. Maybe the most interesting feature from this increased resolution run is observed in the contours of density (Figures 3.3.2(d), and 3.3.3(d)). The tilt up and the tilt down in the isopycnals toward the downstream wall of the canyon show remarkable similarities with Hickey's observations in Astoria 3.3. HIGH RESOLUTION RESULTS 81 canyon (compare with upper panel in Figure 1.2.1). This tilt was already clear in Figure 3.1.20 across the middle of the canyon, because the isopycnal displacement in this region is stronger due to enhanced upwelling (upwelling increases toward the head of the canyon). However, with the low resolution run this feature was not observed as clearly at the mouth of the canyon (Figure 3.1.15). FIGURE 3.3.2. Results from the high resolution run (480x164x64 points) at the mouth of the canyon (Line A), (a) Azimuthal velocity, (b) Radial velocity, (c) Vertical velocity, (d) Density. Velocity units are in m - s , and density is in a . Depth and horizontal distances are in meters. -1 t Summarizing, although the flow is resolved with more accuracy when the spatial resolution is increased, and the evolution of the density field is improved, there are no significant improvements upstream of the canyon, which is the area of biggest concern. This clearly indicates that if there are truncation errors in the parameterization of the bottom boundary layer, they are not influential for the general circulation pattern upstream of the submarine canyon being modeled. 3.4. NON-HYDROSTATIC R U N 82 FIGURE 3.3.3. Results from the high resolution run (480x164x64 points) at the middle of the canyon (Line B). (a) Azimuthal velocity, (b) Radial velocity, (c) Vertical velocity, (d) Density. Velocity units are in m - s , and density is in a . Depth and horizontal distances are in meters. -1 t 3.4. Non-hydrostatic r u n Perenne et al. (2001a) argue that another possible source of errors, that could lead to the differences between laboratory and numerical results, is non-hydrostatic effects. However, in a later study (Perenne et al., 2001b), the good match obtained between (hydrostatic) numerical and laboratory results for upwelling and downwelling experiments in a rotating tank with a submarine canyon leads them to suggest that if non-hydrostatic effects do occur, they do not significantly modify the flow. However, it should be noted that the canyon modeled by Perenne et al. (2001a,b) is less steep than the canyon used in this study and, as discussed in the previous section, they calibrated their interior mixing and bottom drags to improve the agreement with laboratory observations. Therefore their conclusion in this aspect may not be applicable for our canyon. Allen et al. (2003) argue that the strong stratification in the laboratory tank should be enough to prevent non-hydrostatic flows, and evaluate the strength of the non-hydrostatic terms in the 3.4. NON-HYDROSTATIC R U N 83 momentum equation, to obtain an up-canyon, non-hydrostatic flow of 0.04 c m - s . This error is -1 one order of magnitude too small compared to an observed difference of 0.2 c m - s -1 between the numerical results and the laboratory flow. These limits the possibility of non-hydrostatic effects affecting the observations of the flow in the laboratory. However, this analysis is only an estimate of the magnitude of the non-hydrostatic effects. To further investigate what kind of non-hydrostatic effects are present in the canyon, what their spatial structure and temporal evolution is, and how much they affect the flow pattern, MITgcm was set up using the non-hydrostatic capability with the coarser resolution (243x82x32 points). The non-hydrostatic run is considerably more expensive computationally, because a three dimensional elliptic equation must be inverted to solve for the non-hydrostatic pressure perturbation. The method used in MITgcm is called The Pressure Method (Harlow and Welch, 1965; Williams, 1969). The complete discussion of the algorithm used for this inversion and for the non-hydrostatic formulation on MITgcm is found in Adcroft (1995), and in the MITgcm web site (http://mitgcm.org/). The resultant flow from this run (not shown) is not qualitatively different from the hydrostatic run. Therefore the non-hydrostatic effects do not significantly modify the flow, as predicted by Perenne et al. (2001b), and Allen et al. (2003). Nevertheless, this result is interesting, as studies of a submarine canyon (or flow over rough topography) using a non-hydrostatic model are very limited. These studies include Koscher (2001), who in his undergraduate thesis, explores nonhydrostatic effects in a submarine canyon using the Mesoscale Compressible Community model (MC2, Tanguay et al., 1990), and Legg (2004) who uses MITgcm to investigate the effects of internal tides generated at the shelf break over a corrugated continental slope. For this reason and to evaluate the importance of non-hydrostatic effects in the numerical canyon, the results from the non-hydrostatic run are subtracted from the results of the hydrostatic run. This difference (that represents the non-hydrostatic effects) is shown in Figure 3.4.1. The non-hydrostatic effects are only a very small part of the solution representing about 10% of the horizontal velocities, and a maximum of 15% of the vertical velocity. The azimuthal flow in the non-hydrostatic run is stronger than in the hydrostatic flow in most areas (Figure 3.4.1a). However, off-shore velocities are weaker at the upstream side of the canyon using the non-hydrostatic equations (Figure 3.4.1b). In the hydrostatic run, the vertical velocity is obtained diagnostically from the continuity equation, and in the non-hydrostatic run, a prognostic equation is solved instead (2.1.3). Thus, the fact that the highest relative differences are found in the vertical velocity field is 3.4. NON-HYDROSTATIC R U N 84 not a surprise. However it is very interesting to observe the sequence of positive-negative anomalies that are shown in Figure 3.4.1c. These features suggest the presence of internal waves over the canyon, which seem to be almost completely dominated by a second vertical mode. The signature of this wave-like anomaly can be traced to the density difference in Figure 3.4.Id. d)-°° FIGURE 3.4.1. Differences between hydrostatic and non-hydrostatic runs at the mouth of the canyon (Line A), (a) Azimuthal velocity difference, (b) Radial velocity difference, (c) Vertical velocity difference, (d) at difference. Velocities are in cm-s; density is in at. Depth and horizontal distances are in meters. It is almost certain that when a stratified flow is forced over abrupt topography, internal waves will be generated. As mentioned above when discussing the importance of turbulence induced drag in submarine canyons, these are known places for the generation of internal waves. In a continuously stratified medium, the energy associated with a plane internal wave travels at an angle to the horizontal. This angle is determined by the wave frequency, the inertial frequency, and the buoyancy frequency. The slope of the energy flux or group velocity is called the wave characteristic. Internal wave energy incident upon bathymetry that is steeper than the wave characteristic (i.e., 3.4. NON-HYDROSTATIC R U N 85 supercritical bathymetry) will reflect back toward deeper water with some dissipation due to nearbottom mixing. Waves incident on sub-critical bathymetry will be focused into shallower water, with the wavenumber and particle velocity amplitudes increasing in inverse proportion to water depth (Wunsch, 1968). When the slope of the bathymetry is critical (topographic slope equals wave characteristic slope), there is a singularity in the linear solution, the energy density goes to infinity and the group velocity vanishes. Internal wave generation is strongest where the bottom slope is critical for the given wave frequency. In this case, even very weak currents can generate strong internal waves. Hotchkiss and Wunsch (1982) show that for internal waves propagating in a canyon along the axis toward the head, the energy density increases due to the depth change and the converging walls. Thus, submarine canyons can be viewed as conduits that focus and amplify internal wave energy. In an effort to try to understand the spatial, temporal, and horizontal propagation of these waves, in Figure 3.4.2 are contoured the differences at different times between the vertical velocity from the hydrostatic and non-hydrostatic runs at a depth of 1.27 cm. It can be seen that 10 s after the initialization of the model (Figure 3.4.2a), the disturbances are very weak and seem to be radiating from the canyon in every direction, but concentrating on the head. At 20 s (Figure 3.4.2b), as the flow becomes stronger, the disturbances grow and are advected by the mean flow toward the downstream side of the canyon. The peak amplitude of these disturbances occurs at 30 s (Figure 3.4.2c), coinciding with the time of the maximum forcing of the flow. It can be seen that they have propagated further downstream of the canyon. The increase in the amplitude at the downstream corner of the canyon is expected as the waves are propagating from deep to shallow water. At the end of the run (Figure 3.4.2d), the amplitude has decreased as the external forcing has already been removed. From the fact that the crests and troughs of the waves seem to be approximately at the same position relative to the canyon, it is expected that their phase speeds should be on the same order as the speed of the flow at 1.2 cm depth (u~ 1 cm-s ), but moving in the opposite direction. -1 The group speed should be weaker than the flow speed, and this is why these waves get advected downstream. To obtain a better approximation of the phase speed of these waves, (2.3.12) was solved using the final stratification at the upstream corner of the canyon, as in these region the vertical structure of the internal waves looks simpler. The vertical modes were normalized by multiplying them by the maximum vertical velocity divided by the maximum mode value, and then plotted on top of 86 3.4. NON-HYDROSTATIC R U N 0.019 0.0168 0.0147 00126 0.0104 0.0083 0.0062 0.004 0 0019 -0.0002 -0.0024 -0.0045 -0.0066 -0.00B8 -0.0109 -0.013 0.0231 0.0199 0.0167 0.0134 0.0102 0.007 00038 0.0006 -0.0027 -0.0059 -0.0091 -0.0123 -0.0155 -0.0187 -0.022 -0.0252 FIGURE 3.4.2. Top view of the differences in the vertical velocity between hydrostatic and non-hydrostatic runs at 1.2 cm depth, (a) t=10 s. (b) t=20 s. (c) t=30 s. (d) t=35 s. Velocities are in c m s . Depth and horizontal distances are in meters. - 1 the vertical velocity profile. The best matches were obtained for the second baroclinic mode using frequencies between 0.62 and 1.42 s , which correspond to phase speeds between 0.87 and 0.77 - 1 cm-s , and wavelengths of 8 and 3 cm, respectively. This agrees with the observed wavelength of -1 approximately 4 cm. The mode comparison is plotted in Figure 3.4.3a, where the black dashed line is the second baroclinic mode calculated by solving (2.3.12), and the solid blue line is the difference between hydrostatic and non-hydrostatic vertical velocities from MITgcm. Here we can see the very good agreement, especially in the upper part of the water column above the zero crossing. The same procedure was used for the middle part of the mouth of the canyon (Figure 3.4.3b), but the results are not as good. This could be expected given that the vertical modes are calculated using linear theory for a flat bottom, while in the interior of the canyon the topography is far from flat and non-linearities can play an important role (Allen, 1996). In this plot the red dash-dotted line is the sixth baroclinic vertical mode, the black dashed line is the second baroclinic mode calculated 3.4. NON-HYDROSTATIC R U N 87 assuming that the depth was the same as the shelf-break (as if the canyon was removed), and the solid blue line is the non-hydrostatic vertical velocity from MITgcm. b) FIGURE 3.4.3. Comparison between calculated baroclinic vertical mode shapes and vertical velocity profiles from the numerical model at (a) Upstream rim of the canyon, and (b) in the middle of the canyon. Solid blue lines represent the vertical velocity difference between hydrostatic and non-hydrostatic runs in MITgcm, dashed black lines are the second vertical mode, and the red dash-dotted line in (b) represents the 6th vertical mode. Depths are in meters, and velocities in m - s . -1 As it can be seen, the generation, propagation, and dissipation of internal waves in a submarine canyon is not a simple matter, and understanding the complete spatial and temporal structure of these features, although very interesting, is beyond the scope of this work. However, MITgcm proved to be a useful tool to investigate non-hydrostatic phenomena in canyons, and it could be used for this purpose in future research. Furthermore, the results of this test show clearly that nonhydrostatic effects in the laboratory submarine canyon are not strong enough to affect the general pattern of the flow, and are definitely not responsible for the discrepancies between the numerical and laboratory observations. Up to this point, have examined most of the possible causes of error suggested by previous studies as the reasons for the observed differences between numerical and laboratory results. These include errors in the vertical advection of density due to the use of terrain-following numerical models, errors in the parameterization of the bottom boundary layer, and non-hydrostatic effects. However, even though these errors do exist, they not seem to be strong enough to account for the observed differences. In the course of this research, the choice of the advection scheme used for the evolution of the tracer field has proven to introduce changes in the resultant flow. In the 3.5. COMPARISON OF A D V E C T I O N SCHEMES 88 next section, a comparison of the different advection schemes available for the model is conducted focusing on the differences and the impact on the final results. 3.5. Comparison of advection schemes In section 2.2 the possibility of dispersion errors and false extrema in the solution of the advection equation for active tracers and how they depend on the numerical advection scheme used in the discretization of (2.2.6) was described. To illustrate the impact of such errors for the simulations conducted in this work, consider the initial temperature in the surface layer of the model. According to the temperature distribution used to initialize the model (Figure 2.3.6) this temperature is 12.95° C. Now, because in the model there are no additional sources of heat, this temperature should be the maximum observed at any time anywhere in the model simulations. It can be seen from Figure 3.5.1, that depending on the choice of advection scheme, false extrema are observed to appear in this layer. In this figure is plotted the potential temperature at the surface level at the end of the run (35 s). The temperatures are taken across the head of the canyon (inshore of Line B ) . It can be seen that only the advection schemes using a flux limiting algorithm have a smooth solution, and are below the physical limit of 12.95° C. Even though for all the runs using linear advection schemes (2nd order centered, 3rd order up-wind, and 4th order centered) the time step is reduced from 0.025 to 0.01 s in order to ensure that the C F L condition is not violated, these advection schemes develop false extrema and spurious oscillations, especially in the area marked by the canyon width (bold blue lines in Figure 3.5.1). The schemes that use centered differences (2nd and 4th order) show oscillations propagating from the canyon toward the boundaries. Even order methods tend to have dispersion errors such that short wavelength waves propagate upstream against the flow, as observed here. Additionally, Haidvogel and Beckmann (1999) demonstrate that for centered in time (leap frog) centered in space difference schemes, like the ones used here, the finite difference solution always lags the true solution, and that the finite difference solution has a finite phase error no matter what value of At is used. To reduce the amplitude of these oscillations these schemes can be used with a finite amount of diffusion. However, although tests using increased values of the Laplacian diffusion coefficients reduce the appearance of false extrema, and diminish the amplitude of the oscillations, they do not make them disappear completely. Another possibility is introducing bi-harmonic diffusion, but as 89 3.5. COMPARISON OF A D V E C T I O N SCHEMES 13.25j 12.95 12.9 - * - 2nd order centered •••*- 3rd order up-wind ••->-- 4th order centered 2nd order centered with flux limiter — 3-DST ••<->-- 3-DST withfluxlimiter 12.85 12.8 1 50 100 Azimuthal index 150 200 250 FIGURE 3.5.1. Temperature in the surface layer of the tank approximately across the head of the canyon (inshore from Line B) for different advection schemes. The bold blue lines mark the edges of the canyon. The bold black line is the maximum temperature physically possible. Note that the only schemes that are below this limit are the ones that use a flux limiter algorithm. the flow in the laboratory is not turbulent this approach is not used. The principal reason of concern with these schemes, beyond the appearance of false extrema, is that this oscillations are amplified within the canyon limits, casting doubts on any solution obtained for this area using these schemes. Third order up-wind biased schemes (discussed in Section 2.2.3.2) were originally developed to reduce numerical diffusion of low order up-wind schemes and spurious damping and phase errors associated with central difference schemes (Pietrzak, 1998). These improvements are clear in Figure 3.5.1. Small oscillations are still present, with the maximum amplitude of the oscillations located in within the canyon limits. If these are dispersion errors they can pollute the solution of the flow in the area that is the focus of this investigation, however, it is uncertain whether these are errors or naturally occurring oscillations in the tracer field. It is also interesting to find that the solution obtained by using the 3rd order direct space time scheme (3-DST) is virtually the same as that obtained with the 3rd order up-wind biased scheme. Thus, the only advantage derived from the use of 3-DST is the possibility of running the model with a higher Courant number. Although Figure 3.5.1 is useful to observe the presence of spurious oscillations and false extrema, it says little about the impact of these errors in the dynamical response of the flow inside, and around the canyon. The other visualizations of the results used previously (cross sections across Line A and B) are not very useful either, as most of the solutions can be considered to be "close enough" 3.5. COMPARISON OF A D V E C T I O N SCHEMES 90 to the qualitative understanding of the canyon dynamics in these areas. Subtracting one field from the other (like it was done when analyzing the non-hydrostatic results) would imply that one of the schemes is more accurate than the other, and even though the schemes that use flux limiters appear to be better, that cannot be stated yet. Comparison of the particle tracks is not useful either, as the final positions (not shown here) are all too close to allow us to draw any conclusion. The alternative approach used is to investigate the depth of the deepest upwelling observed in the model for each one of the advection schemes. I define the depth of the deepest upwelling as the greatest vertical excursion of a water parcel that is found on the shelf; this means that the water parcel has to be out of the canyon and onto the shelf to be considered upwelled. To identify the area of the greatest vertical excursion, an anomaly in the potential temperature is calculated as the difference between the final and initial temperatures. Recalling that for this work the temperature is an inverted approximation of the density, this anomaly gives a measure of the displacement of the isopycnals as well. Consistent with the conceptual picture described in Section 1.2, the area of the deepest upwelling is located near the head of the canyon at 1.27 cm depth. The results of this comparison are shown in Figure 3.5.2. A negative anomaly indicates that the final temperature is colder than the initial temperature, and thus the water parcel was originally in a deeper layer of the tank. It can be seen that all the schemes show upwelling downstream of the canyon head. As an average, the schemes without flux limiters (Figures 3.5.2a, b, c, and e), show negative anomalies on the order of -3.5°C; this corresponds to water parcels originally located at 2.6 cm depth, and a vertical excursions of 1.3 cm. For comparison, consider that in Barkley canyon, Allen et al. (2001) reports vertical displacements of the isopycnals of 21 to 35 m. These are calculated relative to the isopycnals at the same position on the shelf away from the canyon. Calculating the vertical scale for processes inside the canyon, DH = fL/N (Section 1.2) gives a vertical scale for Barkley Canyon and the laboratory canyon of 138 m and 1.9 cm, respectively. Using these scales, the displacement of the isopycnals in the model is roughly equivalent to a vertical motion of 94.4 m in a canyon such as Barkley Canyon. Considering that a perfect match between observed and modeled vertical excursions cannot be expected due to differences in the scaling of the flow (the Rossby numbers RL = U/fL , and R = ii/ fR are two 0 and 8 times higher in the tank, respectively), this upwelling depth is not unrealistic. Allen and Hickey (in preparation) say that for small R the scale for the upwelling depth is DHV^O • RL, this Q linear scale provides an upper limit to the upwelling in the laboratory of 1.3 cm which is exactly the 3.5. COMPARISON OF A D V E C T I O N SCHEMES 91 maximum upwelling observed in the model. However, the difference between the temperature of the water upwelled by canyon effects and the surrounding water on the shelf (that is moved vertically by coastal upwelling), is not too big and is rather diffuse (i.e. downstream of the canyon head in Figure 3.5.2f). Take for example the temperature anomaly observed when using the 3rd order direct space-time advection scheme, where the water upwelled by normal coastal upwelling induces a temperature anomaly of about -1.96°C on the shelf upstream of the canyon (yellow contour in Figure 3.5.2e), while the highest anomaly induced by canyon effects is -3.27°C (aqua contour downstream of the canyon in Figure 3.5.2e). This gives a maximum difference of 1.31°C between coastal and canyon effects. The advection schemes using a flux limiter (Figures 3.5.2d, and f) show a difference of about 1.96°C between the water upwelled due to canyon effects and the water upwelled by coastal processes in the area downstream of the canyon head. Although the temperature difference is only 0.65 °C greater than that obtained using the 3rd order direct space-time advection scheme, the upwelling due to canyon effects is very easy to identify as it is clearly localized downstream of the canyon head. Using these schemes the deepest upwelled water has a temperature anomaly of -2°C corresponding to an original depth of 2 cm and representing a vertical excursion of 0.75 cm. This is approximately equivalent to 54.5 m in a canyon like Barkley Canyon. This is generally, more consistent with the current understanding of the dynamics of upwelling conditions in submarine canyons. Despite of the qualitative agreement with field observations and the general good match with the tracks of the laboratory particles obtained using advection schemes with flux limiters, there are still systematic differences between the modeled and the observed flow. In the following I propose a possible explanation for these inconsistencies, based on the mechanics of flux limiting algorithms. A non-dissipative, finite-difference advection scheme is unachievable (Shchepetkin and McWilliams, 1998) even without any physically induced dissipation, because the use of a fixed grid induces false oscillations associated with dispersive errors (as it is shown above). To control these unphysical oscillations, unphysical dissipative mechanisms are used. One of the most used of these mechanisms are the non-linear, flux limiting algorithms. Flux limiting algorithms act as locally imposed dissipation that suppresses dispersive overshoots in regions with strong gradients in the advected property. This demands that the algorithm should locally lower its accuracy to first-order upstream (or up-wind), however, global accuracy is in principle unaffected. It is noted that this kind of dissipation can be as strong as necessary, even at the 3.5. COMPARISON OF A D V E C T I O N SCHEMES 92 FIGURE 3.5.2. Temperature anomaly at 1.2 cm depth (this is the depth of the maximum temperature anomaly on the shelf) calculated using different advection schemes, (a) centered second order with Adams-Bashforth. (b) upwind biased third order with Adams-Bashforth. (c) centered fourth order scheme with AdamsBashforth. (d) second order with flux limiter. (e) third order direct space-time method, (f) third order direct space-time method with flux limiter. Temperature anomalies are in ° C . Horizontal distances are in meters. expense of losing accuracy in the solution. The flux limiting algorithm finally used in this study was the Sweby flux limiter (Sweby, 1984) (Figure 3.5.2f) which depends on the local ratio of the advected quantity (for example (2.2.22), and (2.2.23)). 3.5. COMPARISON OF A D V E C T I O N SCHEMES 93 Now, in most of the studies discussed in the background section of this thesis (e.g. Hickey, 1997; Allen et al., 2001; Allen et al., 2003), the explanation proposed for the observed positive vorticity at the upstream side of the canyon in the laboratory and in field observations, relies on the stretching of the isopycnals as water from the shelf falls into the canyon. In this process very strong gradients of temperature (and thus density, and pressure) are generated. In fact, the relation between the distribution of the isopycnals and pressure gradients inside submarine canyons is fundamental to the explanation of the dynamics of the observed upwelling that is related to these topographic features (Allen et al., 2001). Flux limiting algorithms can constrain these temperature gradients in order to avoid over or under shoots of the advected field, and in doing so, they could constrain the stretching of isopycnals necessary for the development of the stretching vorticity on the upstream side of the canyon. A l l numerical models use some form of advection scheme to simulate the evolution of thermodynamic properties. Therefore, when using linear advection schemes or non-linear advection schemes without a flux limiting algorithm, numerical models used to simulate the flow in a submarine canyon are subject to either overshooting, phase or damping errors that pollute the solution in the areas of greatest interest. The alternative of using advection schemes with gradient-limiting or flux-limiting algorithms gives a better solution, but the smoothing effects of flux limiting algorithms can restrict the natural evolution of the flow in the areas where a steep gradient in the solution is really needed. This could explain (at least to some extent) why the eddy above the canyon rim observed in both Astoria Canyon and Barkley Canyon (Hickey, 1997; Allen et al., 2001), or the off-shore velocities at the upstream rim of the laboratory canyon (Allen et al., 2003; this thesis), are not predicted by any published model. To analyze the validity of this argument we have to investigate the effects of the flux limiting algorithm in the canyon in greater detail. A n initial comparison along Line A just above rim depth (1.8 cm) (Figure 3.5.3) indicates that for the low resolution run the differences between the 3-DST scheme, and the 3-DST with flux limiter are not as marked as they should be if the flux limiting algorithm really imposed an important limitation to the stretching of the isopycnals, that is, if the argument proposed above is true. On the other hand, it can be seen from Figure 3.5.3, that the effects of the flux limiting algorithm (dashed blue line) are greater for the high resolution run, smoothing much more the temperature variations at the mouth of the canyon when compared with the same scheme without a flux limiter (magenta line). This could explain why the high resolution 94 3.5. COMPARISON OF A D V E C T I O N SCHEMES 3-DST flux limiting Low res - - 3-DST flux limiting High res — 3-DST Low res 3-psT High res 60 Azimuthal index 80 100 FIGURE 3.5.3. Temperature along Line A above the rim depth (1.8 cm) for the 3-DST advection schemes with and without using a flux limiting algorithm for the low and high resolution runs. The bold blue lines mark the edges of the canyon. results do not introduce major improvements to the results. It is important to make clear that the intention of this argument is not to suggest that the 3-DST without flux limiter does a better job simulating the flow in the submarine canyon. The analysis of Figure 3.5.2 shows that the dynamics of the upwelling is captured better using the flux limiter. Furthermore, the average error for the tracers tracks for the 3-DST without flux limiter (not shown), is only different from the average error for the 3-DST with flux limiter in the third decimal place, and upstream of the canyon the error is less for the latter. The main conclusion that can be drawn from Figure 3.5.3 is that the flux limiter is effectively reducing the steepness of the temperature gradients inside the canyon, and by doing so it can be limiting the stretching of the isopycnals in this region, at least in the high resolution case. To quantify the effects of the flux limiting algorithm at rim depth, the Sweby limiter (2.2.29) used in high resolution run with the 3-DST flux limited scheme, is calculated for the vertical flux, as this flux has the greatest influence on the advection of density. As the vertical index k increases downward, the gradients in the vertical direction are calculated as: 0fc-3/2 ~ flfc-1/2 k 0fc-l/2 r = Afc+1/2 0fc-l/2 _ ~ _ 0fc+l/2 fyc+3/2 0fc+l/2 and the vertical flux is calculated as: (3.5.1) F = - (W + l \w\) (e k+1/2 + V> (0fc-l/2 - Qk+l/2)) + \ (W - M ) (0fc-l/2 - i> ( ) r+ (0fc-l/2 - Ok+l/i) 3.5. COMPARISON OF ADVECTION SCHEMES 95 The flux limiter is plotted in Figure 3.5.4 for the temperature gradient downstream, ip(r~), and upstream ip(r ). + In general, the limiter acts in this way: • When ip(r~) — %jj(r ) = 1, + the limiter reduces the accuracy of the scheme to first order upstream: F = \ ( fc+l/2 + k-l/2) w fl + \ \M {Ok-l/2 e - or F = w6 + k • When ip(r~) = 0, ^\w\6 9 k the limiter reduces the accuracy of the scheme to: F = ^ (w + \w\) 9 k + 1 / 2 + \{w- \w\) (d _ k 1/2 - rf> (r+) (9 _ k 1/2 - 9 )) k+1/2 or F = (w — \w\) 6 + ^ (10 k |io|) V (r ) S 6 + k similarly for ip(r~) = 0: F = (w + \w\) 0 + ^(wk \w\) $ (r~) 5 9 k • When ib(r~) , or ip(r ) are equal to d\ + d r the 3DST advection scheme is recovered. + 2 To understand better the effects of the limiter in the simulations ip(r ), ip(r~), and the flux using the + 3DST scheme with and without flux limiter are evaluated. The calculations are performed using the high resolution results for five different points (Figure 3.5.4) located in the layers immediately above and below the shelf break depth, and the results are summarized in Table 4. The vertical velocities and temperatures used in the calculations correspond to the 3DST flux limited run. Additionally in Table 4, the differential heat fluxes are calculated at the top and bottom of the cell. These are obtained by subtracting the temperature at the middle of the cell from all the temperatures involved in the calculation of the fluxes. • Point 1: There is more flux into the cell than out of the cell for both schemes. The vertical velocity is positive therefore only ip(r~) is active, and as it is equal to d\ + d r then at 2 the top of the cell the flux is the same for both schemes. The velocity at the bottom of the cell is also positive, therefore the cell is cooling. The heat flux at the bottom of the cell is reduced for the 3-DST flux limited scheme (3-DSTfl). As the heat flux coming from below is reduced by the limiter, the cell is cooled more than using the non-limited 96 3.5. COMPARISON OF A D V E C T I O N SCHEMES T A B L E 4. Effects of the flux limiting algorithm for five points (Figure 3.5.4) located in the layers immediately above and below the shelf break depth. Depths are in m, velocities in m - s , and fluxes in ° C - m - s . -1 3 Point Depth w 1 0.019 0.0204 0.019 0.0204 0.019 0.0204 0.019 0.0204 0.019 0.0204 2.5xl0" l.OxlO" -8.5xl0 3-OxlO" -7.5xl0-6.6xl0~ 1.2xl0~ -1.3xl0~ -1.7X10 1.45xl0~ 2 3 4 5 yj(r ) ip(r ) + 0.4 0 5 4 5 0.44 4 5 0 y 10 y 0.57 0.5 4 -5 y 10 0.47 1 b 5.3xl0~ 1.93xl0~ -1.85xl05.88xl0~ -1.3xl0~ -8.9xl0 2.24xl0-2.2xl0 -3.4xl0" 2.67xl0- 3.4X10- IU 9 0 5 5.3xlO~ 1.88xl0-1.85xl05.77xl0-1.3xl0~ -9.0xl02.24xl0~ -2.2xl0~ -3.4xl02.68xl0~ 9 ilJ 10 AFlux 3-DST AFlux 3-DSTfl Flux 3-DST Flux 3-DSTfl iU 0.52 _f) -1 9 y 10 9 - 1 0 y - 9 iU 10 3.4xl0 11 -1.3xl0-1.18xlO-3.17xl0-2.04xlO" 6.41xl01.29xlO1.25xl0- 10 -8.6X10 iU n iU n iU 10 -3.21X10- 3.04xl0~ 11 12 _ i l - 1 1 -1.18xl0-2.11xl0 -2.04xlO~ 7.83xl01.29xlO1.25 x l O iU _11 iU n iU - 1 0 -3.21X10" 2.02xl0" 11 12 scheme. This is confirmed by the differential heat flux, which at the bottom of the cell is negative and stronger for the 3-DSTfl scheme, meaning that more cold water is brought into the cell using this scheme. It is important to note that this result does not contradict the results presented in Figure 3.5.2 where the schemes without flux limiter show stronger coastal upwelling than the schemes that use a limiter. These upwelling calculations show stronger upwelling in the layer that is in contact with the topography, and the heat flux calculations are done one layer above the topography. In fact, the same calculations for the layer below Point 1 (z= 0.0218 m depth), show no flow at the bottom of the cell (the flux at the top is balanced by lateral flux), and a differential heat flux at the top of the cell that is zero for the limited scheme and 4.3 x l O - 1 1 °C-m -s (greater heat flux out of the 3 cell) for the non-limited scheme, meaning that in this layer the cell is cooled more using the non-limited scheme, which is consistent with the results presented in Figure 3.5.2. • Point 2: As the vertical velocity is negative in the upstream corner of the canyon mouth only ip(r ) is active. The negative flux (carrying warmer water) is equal for both schemes. At + the bottom of the cell the flux is also into the cell, only ^ ( r ) is active. If the temperature - gradient is constant, the stronger flux from the top indicates that the cell is warming. The heat flux is reduced for the limited scheme, this indicates that the cell is warms less when using the limiter. The stronger negative differential heat flux at the bottom of the cell shows that more cold water is brought into the cell using the limiter. This confirms that indeed the cell is warmed less using the limited advection scheme. 3.5. COMPARISON OF A D V E C T I O N SCHEMES 97 • Point 3: The vertical velocity is negative at the top and bottom of the cell indicating that the cell is warming. The negative flux at the top of the cell is the same for both schemes. At the bottom of the cell, the flux limiter is enhancing the outward heat flux, therefore the cell should warm less using the flux limited scheme. The differential heat flux at the bottom indicates that less cold water is being taken out of the cell using the flux limiter making the cell become less warm with the 3-DSTfl scheme. • Point 4: Both schemes show the same results. The flux limiter is not active. • Point 5: It was expected that the effects of the limiter would be the same as at Point 2, and judging from the velocity and heat flux this is true. However, the differential heat flux at the bottom is opposite to that at Point 2. The reason for this is that there is a temperature inversion in this area of the domain, and the temperature at the lower level is warmer than the temperature at the upper level. As a result of this, more warm water is brought into the cell using the limited scheme. The model was run without any convective adjustment so whenever an instability occurs, the stratification is not reorganized. Similar inversions were observed in other regions of the domain close to topography, but the biggest temperature differences were observed in the downwelling area downstream of the canyon. There are reasons to believe that this is not a problem with the numerical stability of the model, but a physical phenomena instead, as parcels of fluid in the bottom boundary layer have been observed to be warmer than parcels of fluid immediately above, during relaxation from coastal upwelling over a sloping bottom (Moum et al., 2004). Although this is a very interesting result, its investigation is out of the scope of this thesis and thus it was not pursued any further. Similar tests were repeated in several places throughout the domain giving always similar results. This suggests that the flux limiting algorithm is actually reducing downwelling and thus the stretching of the isopycnals. At the upstream side of the canyon, the reduction of the downwelling (and thus stretching vorticity) could explain at least to some degree the consistent mismatch between laboratory and numerical observations. This analysis suggests that the choice of an advection scheme with a flux limiting algorithm affects the final evolution of the thermodynamic variables in different ways. In general, the flux limiters allow the model to obtain a more reasonable explanation from the dynamical point of view, as is reflected in the analysis of the depth of the deepest upwelling. However, the flux limiting 3.6. VORTICITY FIELD 98 FIGURE 3.5.4. Contour plot of the Sweby limiter for the vertical flux at 2.04 cm depth, ( a ) ^ ( r ) . (b)ip{r ). - + algorithm seems to be constraining the development of stretching vorticity at the upstream side of the canyon, precisely where the thermodynamic gradients are higher. In this manner, flux limiting algorithms can partially be accountable for the differences at the upstream side of the canyon, and hence for the discrepancies between the laboratory and the numerical model results. 3.6. Vorticity field The results from the comparison of the different advection schemes indicate that the schemes using a flux limiting algorithm reduce the stretching of the isopycnals at the upstream rim of the canyon, and therefore this stretching is not enough to generate the strong cyclonic vorticity in this area required for the matching of the particle tracks. To understand how the model is simulating the stretching of the isopycnals, I compare an estimate of the stretching vorticity against the relative vorticity in the canyon. Using the definition given by Hickey (1997) initially described in (1.2.1), the stretching vorticity (SV) is: where the initial (h) and final (h ) thickness of a layer are calculated by comparing the initial and c final temperature distributions. The relative vorticity (RV = dv/dx — du/dy) can be discretized in MITgcm as: (3.6.1) RV = - J - (8iAy v c where A% is the area of the vorticity cell (Figure 2.1.2). SjAx u) c 3.6. VORTICITY FIELD 99 These two quantities are contoured at the canyon mouth in Figure 3.6.1, for the low and high resolution runs. We do not expect these quantities to be the same as in the model run frictional effects are included, but we should expect, that for the same reason, the estimates of stretching vorticity should be higher than the relative vorticity. However, as it can be seen from Figure 3.6.1(a), for the low resolution run, the relative vorticity has higher positive values than the estimated stretching vorticity (Figure 3.6.1(b)), and shows unexpectedly high values of negative vorticity inside the canyon near the walls due to the side and bottom friction that is slowing down the interior eddy. All over the shelf upstream and downstream of the canyon, the relative vorticity shows positive values due to the strong friction on the sloping bottom. The stretching vorticity for the low resolution run shows a very different pattern from the relative vorticity. Here isopycnals are stretched inside the canyon and compressed everywhere at rim depth due to the upwelling of deeper denser water. This is opposite of what the previous discussions suggested about the stretching of isopycnals at the upstream side due to water falling from the shelf into the canyon. These results suggest two possibilities: one, that the model is completely misinterpreting the stretching and compression of the isopycnals, or two, that the frictional effects are strongly overwhelming any other physical processes in the area of study. The latter explanation seems more feasible as frictional effects have already shown a very strong influence on the simulation results. The same calculations using high resolution results (Figure 3.6.1(c), and (d)), show that the stretching vorticity is now higher than for the low resolution run and in at the upstream wall it is even higher than the relative vorticity, also at the upstream rim the isopycnals all seem to be more stretched or less compressed compared to the low resolution results, while the relative vorticity field is practically the same. The reason for this improvement is that using a higher resolution the advection of density is calculated t with more accuracy by the advection schemes, and thus the stretching and compressing of isopycnals is better represented. This also explains why the density contours shown in Figure 3.3.2(d) are much more similar to the observations of Hickey (1997). Now, in addition to the effects of the flux limiting algorithm, if the stretching vorticity does not explain the observed vorticity in the canyon, then the remaining question is: what can explain the vorticity pattern observed in the model? Perenne et al. (2001b) argues that the main mechanism that generates vorticity in their study is flow separation from the upstream rim of the canyon. They get to this conclusion based on the fact that the vortex in their canyon forms first at the upstream side and then is advected into the canyon interior. The same situation occurs in our simulations 3.6. VORTICITY FIELD 100 FIGURE 3.6.1. Vorticity across canyon mouth, a) Relative vorticity low resolution, b) Stretching vorticity low resolution, c) Relative vorticity high resolution, d) Stretching vorticity high resolution. Note that in the high resolution run the positive stretching vorticity is higher than the relative vorticity. Vorticity units are s , horizontal distances and depth in meters. -1 (Figure 3.6.2). They also recall the criterion according to which an along-slope flow should not separate by inertia, if the inertial radius u/f is smaller than the isobath radius of curvature. In our case this can explain the flow detachment from the upstream side of the canyon. Our internal radius of curvature at the shelf break isobath is 1.4 cm. After 10 s from the beginning of the simulation (Figure 3.6.2(a)), at the upstream corner u — 0.65 cm-s , so u/f -1 = 1.25 cm, and therefore the flow is following the upstream isobaths. After 20 s (Figure 3.6.2(b)), u = 1.13 c m - s , u/f = 2.19 -1 cm, the flow is already detached from the upstream corner, and an eddy starts to form. At times 30 and 35 s, the flow separation is evident and the eddy has been advected downstream and covers practically all the length of the canyon. Although the formation of this eddy at rim depth can be then accounted for by flow separation, it must be noted that if an eddy was formed by stretching of the isopycnals it would be advected in exactly the same way into the interior of the canyon. 3.6. VORTICITY FIELD FIGURE 3.6.2. d) t=35 s. 101 Velocity vectors at 2.1 cm depth, a) t=10 s. b) t=20 s. c) t=30 s. It has to be concluded then that in the numerical simulations conducted in this study the main mechanism for the generation of vorticity at the upstream rim of the canyon is also flow separation. This may not be the case in the real ocean or in the laboratory tank, as the observations suggest that an eddy (or off-shore velocities in the tank) develops above the canyon rim depth (Hickey, 1997; Allen et al., 2001; Allen et al., 2003). This eddy should be generated by stretching of the isopycnals, and not by flow detachment from the shelf, as flow separation occurs mainly where the flow is in contact with the topography. Therefore, the differences between the laboratory and the numerical results at the upstream rim of the canyon, can be accounted for by a combination of the use of a flux limiting algorithm for the advection of the active tracers, and the overwhelming frictional effects in the layers in contact with the topography. This finally leads to the poor representation of the stretching of the isopycnals in the numerical model. CHAPTER 4 Summary and Conclusions The dynamic response of a submarine canyon to an impulse-like upwelling-favorable body force was studied using the z-coordinate numerical model MITgcm. The emphasis was directed to the quantitative comparison of the numerical results with laboratory data obtained by Allen et al. (2003); however, qualitative comparisons with field observations were made when it was considered appropriate. The body force was exerted over 30 s and the total duration of the simulation was 35 s, which, when scaled, is equivalent to the duration of an upwelling event in the real ocean. Our results show that near the surface, the flow is basically unaffected by the presence of the canyon. Above the canyon rim near the head, the flow is deflected shoreward, crossing the downstream rim near the head of the canyon. Once on the shelf, there is small return of the flow off-shore probably due to the anticyclonic vorticity from compression of isopycnals in this area, produced by the upwelling of water from within the canyon. Inside the canyon, a cyclonic eddy develops, spanning the whole depth of the canyon. This eddy is stronger, elongated and asymmetrical just below canyon rim depth, and becomes weaker and more symmetric as the depth increases. Upwelling is observed throughout the domain at shelf break depth, but is stronger within the canyon. These numerical results show general qualitative agreement with field studies such as those described in Hickey (1997), and Allen et al. (2001), and results from modeling studies such as Allen (1996), Klinck (1996), Perenne et al. (2001b), and Allen et al. (2003). The comparison of the model with laboratory results shows an improvement in the agreement compared to the previous results obtained using a terrain following coordinate model (Allen et al., 2003). The error in the final positions of the tracers is reduced from 0.65 using the S C R U M model to 0.55 using MITgcm. As both models use basically the same configuration, it is concluded that the improvements observed are due to the use of a z-coordinate model as opposed to a terrain-following model. Other improvements related to the use of MITgcm over S C R U M , include: 102 4. S U M M A R Y A N D CONCLUSIONS 103 • A more realistic representation of the vertical velocities near the canyon, showing downward velocities on the upstream side and upward velocities on the downstream side of the canyon, which is consistent with the theory of the dynamics of a submarine canyon. • The presence of negative radial velocities at the upstream rim of the canyon. Although reduced, consistent differences still exist between the laboratory and numerical results at the upstream rim of the canyon. For this reason other possible sources of errors were explored. These include: errors in the bottom boundary layer parameterization, errors due to non-hydrostatic effects, and errors related to the choice of advection scheme for the density. The bottom drag proved to be of great influence on the results of the simulation. Runs using an arbitrarily reduced drag were performed and showed that the agreement between the model and laboratory results can be significantly increased in this way. It is noted that even though it is possible to calibrate the model by testing different values for the bottom drag and interior mixing, this was not done in the present work, as manipulating the model to obtain good results for the tank does not certify that the model will do well for problems in the real ocean. On the contrary, the importance of comparing a numerical and a laboratory model lies in identifying the differences between them and investigating possible reasons for those differences. Motivated by the important influence that the changes in the bottom drag had on the simulations, and by the results of Perenne et al. (2001a, b), I investigated the possibility of errors in the bottom boundary layer parameterization by running the model with increased resolution using parallel processors. This test was aimed at reducing the size of the bottom cells so that the bottom Ekman layer would be explicitly resolved. However, the results of this test did not give any significant improvement in the simulation, especially in the agreement to the particle tracks. For this reason it is concluded that if these kind of errors do exist they are not strong enough to account for the observed differences between laboratory and numerical results. However, the results from this run showed noticeable similarities with field observations by Hickey (1997), especially in the isopycnal distribution at the mouth, and at the middle of the canyon (Figures 3.3.2d, and 3.3.3d). Next a run using the non-hydrostatic capability of MITgcm was configured in order to investigate the characteristics of the non-hydrostatic effects for the flow in the canyon. The results from this test did not improve the agreement with laboratory observations. Thus, although non-hydrostatic effects are present in the flow, they do not significantly modify the general pattern of the flow, and thus cannot be responsible for the observed differences with laboratory measurements. Nevertheless, 4. S U M M A R Y A N D CONCLUSIONS 104 some interesting non-hydrostatic features are observed. A packet of standing internal waves was identified. These waves are generated at the upstream rim of the canyon and get advected downstream by the mean flow. Estimates, using the linear vertical mode equation for internal waves, showed good agreement with the observed waves for a second baroclinic vertical mode, with corresponding phase speeds on the order of 0.8 cm-s , and wavelengths of approximately 4 cm (about -1 half the canyon width at the mouth). The choice of advection scheme for the active tracers in the model (salinity and temperature) also proved to have an influence on the final results of the simulations. Test runs were performed using the following advection schemes: second order centered (2-CO), third order upwind (3-UW), fourth order centered (4-CO), second order with flux limiter (2-COfl), third order direct space time (3-DST), and third order direct space time with flux limiter (3-DSTfl). In these tests, all the schemes without a flux limiting algorithm developed false extrema in the temperature field. The even order schemes (2-CO and 4-CO) also showed dispersion errors. A l l these errors had an important impact on determining the depth of the deepest upwelling observed in the model. Using the schemes without flux limiters the maximum vertical excursion of the isopycnals was calculated as 1.3 cm, which is equivalent to 94.4 m in a canyon such as Barkley Canyon. Although considering the different scaling used for the simulation in the tank, this depth is acceptable, the difference between the water upwelled due to the presence of the canyon and the surrounding water on the shelf is not as marked as observations suggest it is. Using the advection schemes with flux limiting algorithms the results are improved. The maximum upwelling was calculated as 0.75 cm, which is equivalent to about 54 m in the real ocean, and the difference between the water upwelled by the canyon and the water upwelled by coastal processes is as marked as the observations indicate. Flux limiting algorithms are designed to prevent under or overshoots in the regions where the gradients of the advected quantity are steep. This raised concerns about the advection of density near the canyon, especially at the upstream rim, since this is a region of steep gradients. If a naturally occurring gradient in this area is unnecessarily reduced, the stretching of isopycnals (which is considered one of the main reasons for the observed vorticity in canyons) would be limited. This limitation could constrain the development of an eddy above the rim of the canyon or even the limit the increase in negative radial velocities at the upstream rim that is needed for the agreement between numerical and laboratory results. To explore this possibility the limit functions for the 3DSTfl scheme were analyzed. Plots of these functions showed that the areas with stronger limitation 4. S U M M A R Y A N D CONCLUSIONS 105 were indeed located near the rim of the canyon. The results of this test suggest that the action of these effects limits the warming of the cells at the upstream rim (and thus constraining the stretching of the local isopycnals). For this reason it is concluded that the choice of advection scheme does have an important influence in the physical response of the flow in the canyon and some of the effects of the flux limiting algorithms can be related to the observed discrepancies between laboratory and numerical results. Comparison between estimated stretching vorticity and observed relative vorticity indicate that the stretching and compressing of the isopycnals is not accurately represented in the numerical model. This may be due to the overwhelming effects of friction on the numerical results. This effects should be constrained to the boundary layers in the model but they are effectively seem to be reaching far beyond that. By analyzing the time evolution of the canyon eddy it is concluded that, in the numerical simulations, the main mechanism for the generation of vorticity at the upstream rim of the canyon is flow separation. However, even though this is what is happening in the numerical model, this eddy formation mechanism it is not believed to be the primary mechanism in the tank. In the tank the stretching of the isopycnals should be the reason for the observed off-shore velocities above the upstream rim of the canyon, as flow separation occurs in the layers in contact with the topography and not above that. In general, observations in the ocean (Hickey, 1997; Allen et al., 2001), and in the laboratory (Allen et al., 2003) indicate that the stretching and compressing of isopycnals are one of the main mechanisms for the generation of vorticity in submarine canyons. However the inability of numerical models to accurately represent this physical process, the adverse effects of the use of flux limiting algorithms for the advection of the thermodynamic variables, and the strong influence of frictional effects, lead to discrepancies. This is especially true at the upstream area above the rim of the canyon, where the stretching of the isopycnals is crucial for the development of cyclonic vorticity. More work is needed to understand the physical and numerical processes that result in this shortcoming of the numerical models, and not to limit oneself to arbitrarily manipulating parameters (i.e. the bottom friction coefficient) to obtain a satisfactory result. Bibliography [1] Adcroft, A. J., C. N. Hill, and J. Marshall, 1997: Representation of topography by shaved cells in a height coordinate ocean model. Mon. Wea. Rev. 125: 2293-2315. [2] Adcroft, A., 1995: Numerical algorithms for use in a dynamical model of the ocean, PhD thesis. Imperial College, London. [3] Allen, S. E., 1996 : Topographically generated subinertialflowswithin a finite length canyon. J. Phys. Oceanogr. 26 (8): 1608-1632. [4] Allen, S. E., 2000 : On subinertial flow in submarine canyons : Effect of geometry. J. Geophys. Res. 105: 1285 1297. [5] Allen, S. E., C. Vindeirinho, R. E. Thomson, M. G. G. Foreman, and D. L. Mackas, 2001: Physical and biological processes over a submarine canyon during an upwelling event. Can. J. Fish. Aquat. Sci. 58: 671-684. [6] Allen, S. E . , M. S. Dinniman, J. M. Klinck, D. D. Gorby, A. J. Hewett, and B. M. Hickey, 2003: On vertical advection truncation errors in terrain-following numerical models : Comparison to a laboratory model for upwelling over submarine canyons. J. Geophys. Res. 108 (Cl): 3003. [7] Bitz, C.M. and W.H. Lipscombe, 1999: An energy-conserving thermodynamic model of sea ice. J. Geophys. Res. 104: 15,669 - 15,677. [8] Chen, X., and S. E . Allen, 1996: Influence of canyons on shelf currents - a theoretical study. J. Geophys. Res. 101: 18043-18059. [9] Cherubin L., X. Carton, J. Paillet, Y. Morel & A. Serpette, 2000: Instability of the Mediterranean water undercurrents southwest of Portugal: Effects of baroclinicity and of topography. Oceanologica Acta. 23(5): 551573. [10] Cherubin L. M., N. Serra and I. Ambar, 2003: Low-Frequency variability of the Mediterranean undercurrent downstream of Portimao canyon. J. Geophys. Res. 108(C3): 3058, doi:10.1029/2001JC001229. [11] Cox, M. D., 1987: Isopycnal diffusion in a z-coordinate ocean model. Ocean modeling. 74: 1-5. [12] Danabasoglu, G., and J. C. McWilliams, 1995: Sensitivity of the global ocean circulation to parameterizations of mesoscale tracer transports. J. Climate. 8: 2967-2987. [13] Douglas, C. C , G. Haase, and S. Reitzinger, 2003: Special solution strategies inside a spectral ocean element model. Math. Models Meth. Applied Sci. 13: 1-14. [14] Ekman, V. W. 1905: On the influence of the earth's rotation on ocean currents. Arkiv. Matem., Ast. Fysik. 2:1-52. 106 BIBLIOGRAPHY 107 [15] Follows, M . J., T. Ito, and J. Marotzke, 2002: The wind-driven, subtropical gyres and the solubility pump of CO-2, Global Biogeochem. Cycles, 16(4), 1113. [16] Freeland, H. L., and K. L. Denman, 1982: A topographically-induced upwelling center off southern Vancouver Island. J. Mar. Res. 40: 1069-1093. [17] Gardner, W. D., 1989: Periodic resuspension in Baltimore Canyon by focussing of internal waves. J. Geophys. Res. 94: 18185-18194. [18] Gent, P. R., J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr. 20 (1): 150-155. [19] Haidvogel, D. B., J. L. Wilkin, and R. Young, 1991: A semi-spectral primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear horizontal coordinates. J. Comput. Phys. 94: 151-185. [20] Haidvogel, D. B., and A. Beckmann, 1999: Numerical ocean circulation modeling. Imperial College Press, London. [21] Hall G., and J.M. Watt, 1976: Modern numerical methods for ordinary differential equations, Claredon Press, Oxford. [22] Harlow, F. H., and J.E. Welch. 1965: Time dependent viscous flow. Phys. Fluids. 8: 2182-2193. [23] Hickey, B. M., 1995: Coastal submarine canyons. Topographic effects in the ocean: Aha Huliko'a Hawaiian winter workshop, Honolulu, HI, University of Hawaii at Manoa. 95-110. [24] Hickey, B. M . , 1997: The response of a steep-sided, narrow canyon to time-variable wind forcing. J. Phys. Oceanogr. 27 (5): 697-726. [25] Hotchkiss, F. S., and C. Wunsch, 1982: Internal waves in Hudson Canyon with possible geological implications. Deep-Sea Res. 29: 415-442. [26] Hundsdorfer W., and R. A. Trompert, 1994: Method of lines and direct discretization: a comparison for linear advection. Appl. Numer. Math. 13: 469-490. [27] Koscher, C. M. 2001: Adaptation of the MC2 model to simulate steep-sided submarine canyon flow and determination of non-hydrostatic effects, Honours thesis. University of British Columbia. [28] Kinsella, B. D., A. E. Hay, and W. W. Denner, 1987: Wind and topographic effects on the Labrador current at Carson canyon. J. Geophys. Res. 92: 10853-10869. [29] Klinck, J.M., 1989, Geostrophic adjustment over submarine canyons, J. Geophys. Res. 94: 6133-6144. [30] Klinck, J. M., 1996: Circulation near submarine canyons: A modeling study. J. Geophys. Res. 101: 1211-1223. [31] Kunze, E., L. K. Rosenfeld, G. S. Carter, and M. C. Gregg, 2002: Internal waves in Monterey Submarine Canyon. J. Phys. Oceanogr. 32 (6): 1890-1913. [32] Legg, S. 2004: Internal tides generated on a corrugated continental slope. Part I: Cross-slope barotropic forcing. J. Phys. Oceanogr. 34 (1): 156-173. [33] Large, W. G., J . C. McWilliams, and S. C. Doney, 1994: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Rev. Geophys. 32: 363- 403. BIBLIOGRAPHY 108 [34] Lueck, R. G., and T. R. Osborn, 1985: Turbulence measurements in a submarine canyon. Cont. Shelf Res., 4: 681-698. [35] MacCready, P., and P. B. Rhines, 1991: Buoyant inhibition of Ekman transport on a slope and its effect on stratified spin-up. J. Fluid Mech., 223: 631-661. [36] Marshall, J., A. Adcroft, C. Hill, L. Perelman, and C. Heisey, 1997(a): A finite-volume, incompressible navier stokes model for studies of the ocean on parallel computers. J. Geophys. Res. 102(C3): 5753-5766. [37] Marshall, J., C. Hill, L. Perelman, and A. Adcroft, 1997(b): Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. J. Geophys. Res. 102(C3): 5733-5752. [38] Maso, M., La Violette, PE and Tintore, J.,1990: Coastal flow modification by submarine canyons along the NE Spanish Coast. Sci. Mar. 54 (4): 343-348. [39] McKinley, G. A., M . J. Follows, and J. Marshall, 2000: Interannual variability of the air-sea flux of oxygen in the North Atlantic. Geophys. Res. Let. 27: 2933- 2936. [40] Mirshak, R., 2001: Spin-up over steep topography and the effects of a submarine canyon. MSc thesis. University of British Columbia. [41] Mourn, J. N., A. Perlin, J. M. Klymak, M. D. Levine, T. Boyd, P. M . Kosro, 2004: Convectively driven mixing in the bottom boundary layer. J. Phys. Oceanogr., 34 (10): 2189-2202. [42] Perenne, N., D. B. Haidvogel, and D. L. Boyer, 2001(a): Laboratory - numerical model comparisons of flow over a coastal canyon. J. Atmos. Oceanic Technol. 18: 235-255. [43] Perenne, N., J. W. Lavelle, D. C. Smith, D. L. Boyer, 2001(b): Impulsively started flow in a submarine canyon: comparison of results from laboratory and numerical models. J. Atmos. Oceanic Technol. 18: 1698-1718. [44] Petruncio, E. T., L. K. Rosenfeld, and J. D. Paduan, 1998: Observations of the internal tide in Monterey Canyon. J. Phys. Oceanogr. 28 (10): 1873-1903. [45] Pietrzak, J. 1998: The use of TVD limiters for forward-in-time upstream-biased advection schemes in ocean modeling. Mon. Wea. Rev. 126 (3): 812-830. [46] Press, W. H., S. A. Teulosky, W. T. Vetterling, and B. P. Flannery, 1992: Numerical recipes in Fortran 77: The art of scientific computing. Cambridge University Press. 963 p. [47] Redi, M . H., 1982: Oceanic isopycnal mixing by coordinate rotation. J. Phys. Oceanogr. 12 (10): 1154-1158. [48] Roe, P. L., 1985: Some contributions to the modeling of discontinuous flows. Lect. Notes Appl. Math. 22:163-193. [49] Shchepetkin, A. F., and J. C. McWilliams, 1998: Quasi-monotone advection schemes based on explicit locally adaptative dissipation. Mon. Wea. Rev. 126: 1541-1580. [50] Shepard, F.P., N . F. Marshall, P. A. McLoughlin, 1974: Currents in submarine canyons. Deep-Sea Res. 21: 691-706. [51] Song, Y . H., and Haidvogel, D. B., 1994: A semi-implicit ocean circulation model using a generalized topographyfollowing coordinate system. J. Comput. Phys. 115: 228-244. [52] Sweby, P. K . 1984: High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. Num Anal. 21: 995-1011.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical simulation of flow in a laboratory tank using...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical simulation of flow in a laboratory tank using a z-coordinate numerical model Jaramillo, Sergio 2005
pdf
Page Metadata
Item Metadata
Title | Numerical simulation of flow in a laboratory tank using a z-coordinate numerical model |
Creator |
Jaramillo, Sergio |
Date Issued | 2005 |
Description | A z-coordinate numerical model is employed to study the dynamic response of an upwelling favorable flow in a rotating laboratory tank with a submarine canyon. A vertically uniform body force is exerted on the fluid inside the tank over a time of 30 s, equivalent to the duration of an upwelling event in the real ocean. The results show a flow pattern that is in qualitative agreement with previous numerical, theoretical and field studies. The agreement between laboratory and numerical results is an improvement on previous results obtained using a terrain-following numerical model. However consistent differences at the upstream rim of the canyon exist. A number of processes are investigated to understand the reason for these differences. Results show that the model is extremely sensitive to changes in the bottom friction and interior mixing coefficients. Although the model can be "tuned" to obtain better agreement, this approach is not used because it could mask other important problems. Instead, results from an increased resolution run using distributed processors show that the observed differences are not related to errors in the bottom boundary layer parameterization. Non-hydrostatic effects are estimated, showing a package of standing internal waves that is generated at the upstream side of the canyon. These internal waves are dominated by a second baroclinic vertical mode that is well predicted by the linear theory. However the non-hydrostatic effects are not strong enough to influence the main characteristics of the flow in the canyon. Different advection schemes for the active tracers are compared, showing that non-linear schemes with a flux limiting algorithm produce the best results from a physical point of view. Flux limiters are designed to regulate the flow in areas with steep gradients of the advected quantity; the influence of these algorithms is investigated and it is shown that they constrain the stretching of isopycnals at the upstream rim and thus they account at least in part for the observed errors in this area. Comparisons between stretching and relative vorticity show that the model results are overwhelmed by frictional effects and that the stretching and compressing of isopycnals is poorly represented. Through these modeling steps, it is concluded that the vorticity observed in the numerical model is more related to flow separation than to stretching vorticity effects. This also contributes to the observed differences between the numerical and laboratory results. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052361 |
URI | http://hdl.handle.net/2429/16521 |
Degree |
Master of Science - MSc |
Program |
Oceanography |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2005-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2005-0228.pdf [ 18.83MB ]
- Metadata
- JSON: 831-1.0052361.json
- JSON-LD: 831-1.0052361-ld.json
- RDF/XML (Pretty): 831-1.0052361-rdf.xml
- RDF/JSON: 831-1.0052361-rdf.json
- Turtle: 831-1.0052361-turtle.txt
- N-Triples: 831-1.0052361-rdf-ntriples.txt
- Original Record: 831-1.0052361-source.json
- Full Text
- 831-1.0052361-fulltext.txt
- Citation
- 831-1.0052361.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052361/manifest