UBC Faculty Research and Publications

Solution convergence of flow over steep topography in a numerical model of canyon upwelling. Dawe, Jordan T.; Allen, Susan E. May 31, 2010

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

Item Metadata


52383-Allen_AGU_2009JC005597.pdf [ 907.13kB ]
JSON: 52383-1.0041906.json
JSON-LD: 52383-1.0041906-ld.json
RDF/XML (Pretty): 52383-1.0041906-rdf.xml
RDF/JSON: 52383-1.0041906-rdf.json
Turtle: 52383-1.0041906-turtle.txt
N-Triples: 52383-1.0041906-rdf-ntriples.txt
Original Record: 52383-1.0041906-source.json
Full Text

Full Text

Click Here  JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 115, C05008, doi:10.1029/2009JC005597, 2010  for  Full Article  Solution convergence of flow over steep topography in a numerical model of canyon upwelling Jordan T. Dawe1 and Susan E. Allen1 Received 23 June 2009; revised 10 December 2009; accepted 28 December 2009; published 11 May 2010.  [1] The convergence as resolution is increased is examined for a numerical model simulation of upwelling over a continental shelf canyon. A series of idealized one‐ and two‐dimensional models are used to plot the dependence of bottom Ekman layer velocity structure and transport on model resolution. Using these results as a guide, a bottom Ekman layer‐resolving three‐dimensional numerical model is constructed. This model is used to simulate a laboratory model of canyon upwelling. Alongshore velocity, volume of fluid upwelled through the canyon onto the shelf, and change in depth of the canyon water are examined at three grid resolutions. Numerical model particle trajectories are compared with a laboratory model of canyon upwelling to validate the model. The results suggest that current computational cluster power is sufficient to accurately simulate all aspects of upwelling in a steep canyon with the exception of flow separation from the canyon wall. Citation: Dawe, J. T., and S. E. Allen (2010), Solution convergence of flow over steep topography in a numerical model of canyon upwelling, J. Geophys. Res., 115, C05008, doi:10.1029/2009JC005597.  1. Introduction [2] The interaction of flows with steep topography is important to a wide range of ocean processes, including continental shelf fluid exchange [Allen, 2004], tidal energy conversion [Laurent et al., 2003], and deep ocean mixing [Kunze and Smith, 2004]. Canyons on the edge of the continental shelf are a good test case for the study of these kinds of flow. First, canyons enhance shelf exchanges during upwelling or downwelling events, which has impacts on shelf biology and deep ocean water properties [Allen et al., 2001; Allen, 2004; Bosley et al., 2004; Ladd et al., 2005; Williams et al., 2006], making the problem of interest. Second, the importance of canyons in shelf exchange processes means that there have been several field studies done on a variety of coastal canyon systems which give a reasonable picture of the dynamics of the system [Hickey, 1997; Sobarzo et al., 2001; Wang et al., 2003; Palanques et al., 2005]. Third, coastal canyon systems are easily represented in a laboratory tank model [Pérenne et al., 2001a; Allen et al., 2003; Boyer et al., 2006], giving us an idealized theoretical test bed to compare with model simulations. [3] The elements needed to represent a canyon system are a continental shelf cut by a canyon, density stratification, and an alongshore geostrophic current. Coastal upwelling or downwelling events in the ocean are triggered when the wind and Coriolis forces conspire to push fluid toward or away from the land. This sets up a barotropic pressure gra1 Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada.  Copyright 2010 by the American Geophysical Union. 0148‐0227/10/2009JC005597  dient that drives a geostrophic current that runs parallel to the coastline. Friction then creates a bottom Ekman layer that transports fluid onshore or offshore, relaxing the barotropic pressure gradient [Thorpe, 1987]. However, fluid that is transported up or down the steep continental shelf slope advects the vertical density gradient and creates a baroclinic pressure gradient which eventually arrests the bottom Ekman flow [MacCready and Rhines, 1993]. [4] Canyons serve to enhance these events in two ways: first, they act like a ramp and advect fluid less vertically than the shelf, delaying the shutdown of the boundary layer; second, the walls of the canyon prevent alongshore flow within the canyon, breaking geostrophic balance and exposing the interior of the canyon to the barotropic pressure gradient [Freeland and Denman, 1982; Klinck, 1996]. Canyons can further be divided into “long” canyons, in which the breaking of geostrophy within the canyon is largely accomplished by convergence and divergence of the topographic contours, and “short” canyons, in which it is the result of advection [Allen, 2000; Waterhouse et al., 2009]. [5] Laboratory models that include all these effects can be created using a cylindrical rotating tank containing a cylindrically symmetrical topography representing the continental shelf and canyon and filled with a vertically stratified fluid [Pérenne et al., 1997]. Alongshore currents can be created by varying the tank rotation rate [Pérenne et al., 2001a]. Several studies have been performed using this type of setup and have been compared with numerical models [Boyer et al., 2000; Pérenne et al., 2001b; Allen et al., 2003; Boyer et al., 2004; Jaramillo, 2005; Haidvogel, 2005]. The results, in general, show good agreement between tank and numerical models for canyons with shallow slopes and poor agreement for canyons with steep slopes.  C05008  1 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  Figure 1. Model topography viewed (top) in cross section along the canyon axis and (bottom) from above. The tank is 1 m in diameter and filled to depths of 9.1 cm over the abyss, 2.1 cm at the shelf break, and 3 mm at the outside wall. The base of the shelf slope has a radius of 20.5 cm, and the shelf break has a radius of 27.5 cm. The central cylinder added to the numerical model to ensure numerical stability is indicated by dashed lines and is 10 cm in radius. The dotted line indicates the height of the canyon walls. [6] These continuing problems in accurately representing the flow field in these models are the motivation for the following work. Here we concentrate on a z level model, the Massachusetts Institute of Technology general circulation model (MITgcm) [Marshall et al., 1997], using it to simulate a laboratory tank experiment of an upwelling event over a continental shelf canyon. By running the model in a variety of configurations we attempt to estimate the grid spacing necessary to resolve the numerical solution and identify which dynamical features are hardest to resolve. Then we use data from the laboratory experiment to examine how well our numerical solution represents reality. Our results suggest that the accuracy of numerical models of flow over steep topography may be limited by the impact of small‐scale nonlinear processes such as flow separation.  2. Two Models [7] In this study we discuss two models: a laboratory model of an upwelling continental shelf canyon system and a numerical model simulation of the laboratory model. The laboratory model data we refer to are the result of previous work by Allen et al. [2003], Mirshak and Allen [2005], Bowie [2006], Waterhouse [2006], and Waterhouse et al. [2009].  C05008  2.1. Laboratory Model [8] The laboratory model consists of a 1 m diameter cylindrical Plexiglas tank placed upon a rotating table. The tank contains a radially varying topography intended to simulate a continental shelf (Figure 1). A 22° section of this shelf is designed to be removable to accommodate a variety of canyon inserts. [9] A vertical density gradient is created in the tank by varying the salinity of the water used to fill the tank using the two‐bucket method of Oster [1965] [see also Hill, 2002]. This method uses two buckets, one filled with fresh water and the other filled with salt water, and the two are connected by a small tube with a valve that is initially closed. Fluid from the freshwater bucket is slowly drained into the laboratory tank while the valve between the buckets is opened, allowing salt water into the freshwater bucket where a mechanical stirrer ensures that the bucket is well mixed; thus, the fluid entering the tank increases in salinity linearly with time. The initial salinity of the saltwater tank was chosen so that the buoyancy frequency N at the model shelf break would be 2.2 s−1. [10] The tank is filled while rotating with Coriolis parameter f = 0.4 s−1 and once filled is left to spin up for 3 h to allow the fluid to equilibrate into solid body rotation. Upwelling events are simulated in the tank by varying the rotation rate of the table to create an alongshore flow. During an experiment the tank’s rotation rate is increased from f = 0.4 s−1 to f = 0.52 s−1 with a constant angular acceleration over a period of 27.3 s, a duration chosen to simulate a day‐long upwelling event. Data were collected immediately after the rotational acceleration ended. 2.2. Numerical Model [11] We use the MITgcm, a primitive equation z coordinate model, to construct a numerical model of the laboratory tank. Since the flow in the tank is laminar, we set the model viscosity and diffusion coefficients to the molecular values for water: n = 10−6 m 2 s −1 for viscosity and  = 1.5 × 10−9 m2 s−1 for salinity diffusion. Temperature is set to be constant, and the model is initialized with a vertical salinity gradient similar to that which would be produced via the Oster [1965] two‐bucket method, modified by 3 h of vertical diffusion (Figure 2). A linear equation of state is used for density calculations. [12] The model is configured to run in nonhydrostatic mode in order to capture internal wave dynamics. However, comparison of the nonhydrostatic runs presented here with equivalent hydrostatic runs shows essentially identical results in all model fields, save for the vertical velocity. Maxima in the hydrostatic vertical velocity field tend to lie directly above or below minima, while the nonhydrostatic field has a more complex spatial pattern. However, the vertical velocity differences occur in regions where the vertical velocity is small and thus have little effect on the simulation. Nevertheless, we use the nonhydrostatic dynamics to ensure that our simulation is as realistic as possible. [13] The model is run in cylindrical coordinates to take advantage of the cylindrical symmetry of the tank. Measurements of the laboratory tank are used to create the numerical model topography, the accuracy of which is  2 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  also be incorrect, and this will affect the strength of the pressure gradients which drive upwelling on the shelf and in the canyon. Furthermore, the Ekman transport also affects the advection of density gradients up the continental shelf slope. This in turn will alter the rate at which the Ekman transport is shut down by the baroclinic pressure gradient created as dense water advects up the slope [MacCready and Rhines, 1993]. In sections 3.1 and 3.2 we use a set of simple one‐ and two‐dimensional models to estimate the model resolution needed to resolve bottom Ekman transports. 3.1. One‐Dimensional Model [17] Consider an infinite f plane with a steady geostrophic flow U in the x direction, with a no‐slip condition at the bottom. The equations governing the evolution of this system are  Figure 2. Model initial density profile. enhanced through the use of MITgcm’s partial grid cell representation of topography [Adcroft et al., 1997]. The tank topography is modified in the numerical model by the addition of a 10 cm radius cylinder in the center of the tank, in order to avoid Courant‐Friedrichs‐Lewy criterion violations where the cylindrical grid spacing is small. Friction with the wall of this cylinder is disabled to minimize its influence on the model flow. A variety of different grid spacings were used in this study to examine the convergence of the numerical simulation. [14] To simulate the change in rotation rate in the tank experiments that drives the tank flow, an azimuthal body force is applied to the tank for 27.3 s. The Coriolis parameter is maintained at f = 0.52 s−1 throughout the model run; analytical calculations indicate that this constant Coriolis parameter results in negligible errors.  3. Resolving the Ekman Layer [15] If one ignores nonlinear effects, such as flow separation within the canyon, the smallest scale that must be resolved in the tank is the viscous bottom Ekman layer thickness. This thickness is set by the Ekman depth, sffiffiffiffiffi 2 ¼ ; f  ð1Þ  which equals 1.96 mm for our values of n and f. Thus, one might expect that a vertical resolution of the order of 2 mm would be necessary to resolve the velocities in the Ekman layer. [16] However, the velocity structure of the Ekman layer is of secondary importance compared to the amount of fluid transported. Ekman transport is the primary mechanism by which fluid moves from areas of high pressure to low pressure, reducing surface height gradients and spinning down the geostrophic flow. If the Ekman transport is incorrect, the spin‐down rate of the geostrophic flow will be  @u 1 @P @2u ¼ fv À þ 2 @t 0 @x @z  ð2Þ  @v 1 @P @2v ¼ Àfu À þ 2; @t  @y @z  ð3Þ  where u and v are the velocities in the x and y directions, respectively; P is the pressure; and r is the density of the fluid. The initial conditions for this system are then u = U, v = 0, and P = P0 − rfUy, where P0 is a constant dependent on the depth of the fluid. Bottom friction will quickly result in the creation of an Ekman velocity spiral, reducing the bottom u velocity and driving a cross‐flow v velocity within the frictional boundary layer. [18] Discretization of these equations in the z direction and integration forward in time with a fourth‐order Runge‐ Kutta scheme allow us to examine the effect of vertical resolution on the representation of this Ekman layer. A 2 cm deep model was configured with the velocity at the top of the domain fixed to u = U and the viscosity, density, and f parameters set to appropriate values for the laboratory tank. The model was then run with a variety of vertical resolutions, from Dz = 2 cm to Dz = 0.02 mm, for 300 s of model time (about 12 inertial periods). This was sufficient to let the solution converge to a steady state where velocities were changing <0.01% s−1. [19] Comparison of the theoretical solution of the bottom Ekman boundary layer [Kundu, 1990] with a low‐resolution (Dz = 2 mm) run shows that while the velocity structure of the run is reasonable, the cross‐geostrophic flow Ekman transport is overestimated (Figure 3a). Calculating the model’s Ekman transport as a function of resolution (Figure 3b) indicates that in order for the numerical transport to be within 1% of the theoretical value a vertical resolution of ∼0.2d is required; for the transport to be within 2%, ∼0.3d is required. 3.2. Two‐Dimensional Model [20] The experimental tank topography is far more complicated than an infinite flat plane, with several changes in slope from the central abyss to the shallow continental shelf. The continental shelf’s 4.5° slope allows for more complicated dynamics such as upslope advection of the vertical density gradient and flow interactions with the edges of the  3 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  transport over the flat abyssal plane relative to the higher‐ resolution model. However, over the sloped continental shelf the situation is reversed, with the low‐resolution model having less Ekman transport than the high‐resolution model. One obvious reason for this is the staircase‐like representation of the smooth continental shelf in the lower‐resolution model, which obstructs flow within the Ekman layer. The effects of these sudden depth transitions are mitigated in MITgcm via the inclusion of a partial grid cell representation of bottom topography, but this is not enough to completely alleviate the problem. [24] A clearer view of how the Ekman transport is affected by resolution and topographic slope can be found by examining how the net Ekman transport varies over the radius of the tank (Figure 5). To calculate the cross‐geostrophic Ekman transport, we assume that all flow directed radially outward within 1 cm of the bottom is in the Ekman layer and sum all positive radial transports over the entire model run. We do this for five of the two‐dimensional model runs: one with low vertical and low horizontal resolution (6.25 mm horizontal × 2 mm vertical), one with high vertical and high horizontal  Figure 3. One‐dimensional Ekman layer model. (a) Steady state solution of one‐dimensional Ekman layer model cross‐ mean flow velocity in mm s−1 with Dz = 2 mm (bars) compared to the analytical solution (line). Transport is represented by the area of the bars or the area inside the curve. (b) Model cross‐mean flow transport at a variety of vertical resolutions, normalized by the theoretical transport. Vertical resolution is given in units of the Ekman depth. shelf and the tank. To simulate these effects, we turn to the MITgcm. [21] For now, we ignore the effect of the canyon and simply look at the Ekman layer in a tank with an unbroken shelf. This allows us to run the MITgcm in a two‐dimensional mode under the assumption that the flow is cylindrically symmetric. This simplification lets us perform many model runs and compare the Ekman layer representation as the model resolution changes. We make a set of runs with a radial resolution varying from 1.25 to 6.25 mm and a vertical resolution varying from 0.4 to 2.0 mm, with an additional ultrahigh‐resolution run done at 0.625 mm horizontal and 0.2 mm vertical resolution in order to judge the convergence of the lower‐resolution models. [22] First, we integrate the radial velocity at the end of the ultrahigh‐resolution run (t = 30 s) from the bottom of the tank to the surface to generate the r‐z stream function. The r‐z stream function shows two distinct overturning cells, one in the abyss and one on the shelf, separated by the steep continental slope (Figure 4a). No bottom layer flow is apparent between these cells as the steepness of the 45° slope quickly results in baroclinic pressure gradients that shut down the Ekman transport. [23] Comparison of the 6.25 mm horizontal × 2.0 mm vertical resolution model stream function with the ultrahigh‐ resolution stream function (Figure 4b) shows several differences. In agreement with the one‐dimensional model results, the 2 mm vertical resolution model overestimates the Ekman  Figure 4. Two‐dimensional radial‐depth stream functions for models run at (a) 0.5 mm horizontal × 0.2 mm vertical and (b) 6 mm × 3 mm grid spacings.  4 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  Figure 5. Bottom Ekman layer transport (m2 s−1, given in units of flux per centimeter of azimuthal distance) as a function of radius for a variety of model resolutions. resolution (1.25 mm horizontal × 0.4 mm vertical), two mixed resolution runs (6.25 mm × 0.4 mm and 1.25 mm × 2 mm), and the ultrahigh‐resolution run (0.625 mm × 0.2 mm). [25] The resulting transports show two maxima, one in the abyss and one over the shelf, and a minima over the steep continental slope. In the abyss the transport shows little sensitivity to horizontal resolution but shows the same overestimation of transport at low vertical resolution that the one‐dimensional model did. At 0.4 mm vertical resolution (0.2d) the transport is almost identical to the ultrahigh‐ resolution run, and the model representation appears to have converged. Over the slope all but the very lowest resolution run appear to have converged, indicating that steep topographic slopes need little resolution to be well resolved, despite the presence of a sizable transport. Even in the low‐ resolution run the transports only differ from the ultrahigh‐ resolution run near the shelf break and appear to be more related to the shelf break Ekman transports which pull water up from the abyss. [26] The shelf transports show the strongest sensitivity to resolution, with the low‐resolution runs underestimating the transport by ∼15%. Increasing either the horizontal or the vertical resolution improves the transports roughly equally. However, the 2 mm vertical resolution models both display a small radial oscillation in the transport strength due to tiny vertical recirculation cells that occur where the flow meets a step in the ill‐resolved topography. This suggests that the vertical resolution is still slightly more important to achieving a realistic representation of the transport. The high vertical and high horizontal resolution run shows good overall agreement with the ultrahigh‐resolution transport, although the transport is still slightly underestimated. [27] Clearly, the shelf Ekman layer is the most difficult part of the flow to resolve, and so we focus on this to help pick a suitable resolution for a full three‐dimensional model. We run the two‐dimensional model at a variety of horizontal and vertical resolutions and measure the total Ekman transport over the shelf at r = 0.3625 m, a point chosen to be near  the middle of the shelf and to be at the minimum of the low‐ resolution transport oscillations, in order to maximize the transport errors. Comparing these values as a percentage of the ultrahigh‐resolution (0.625 mm × 0.2 mm) transport allows us to plot out a two‐dimensional space representing how shelf transport varies with resolution (Figure 6). From this it is clear that both horizontal and vertical resolution are equally important in resolving the shelf Ekman layer flow.  4. Representation of Upwelling in a Continental Shelf Canyon [28] Having examined the grid spacing needed to resolve two‐dimensional radial Ekman transport over a model continental shelf, we move to a full three‐dimensional simulation  Figure 6. Net midshelf bottom Ekman layer transport as a function of vertical and horizontal grid resolution. Transports are normalized by the model transport at a resolution of 0.625 mm horizontal × 0.2 mm vertical.  5 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  driving a baroclinic pressure gradient that creates the offshore jet; increased Ekman transport thus results in a stronger offshore jet. [30] Between the 5 × 2.8 mm model (Figure 8a) and the 3 × 1 mm run (Figure 8b), the strength of the offshore jet increases by 2%, while the velocity over the shelf decreases by 3%–10%. Between the 3 × 1 mm (Figure 8b) and 1.5 × 0.5 mm (Figure 8c) runs, the changes are much smaller, with a <1% change in the jet velocity and a reduction of 1%–4% over the shelf. Furthermore, the velocity structure over the shelf in the 5 × 2.8 mm model appears to be significantly influenced by the low‐resolution step changes in the shelf topography. This effect is much smaller at 3 × 1 mm reso-  Figure 7. Stretched grid used in the 3‐D MITgcm model runs showing the concentration of resolution over the canyon region. The grid shown here only has a resolution of 12.5 mm, for clarity. The grid cells were evenly subdivided to produce the higher resolutions used in the model runs. of a continental shelf canyon upwelling system. To minimize our computational requirements, we employ a stretched grid in the azimuthal direction in order to concentrate model resolution over the canyon (Figure 7). Grid resolution is unevenly distributed, with more grid cells downstream of the canyon to better represent the canyon’s wake. The spacing of the stretched grid was chosen so that grid cells would be approximately square over the canyon. Three main runs are performed with grid sizes of 64 azimuthal × 81 radial × 32 depth (5 mm horizontal grid spacing × 2.8 mm vertical grid spacing), 128 × 134 × 91 (3 mm × 1 mm), and 256 × 268 × 182 (1.5 mm × 0.5 mm). Referring to Figure 6, these resolutions resulted in Ekman transports of ∼85%, ∼92.5%, and ∼97% of the ultrahigh‐resolution run. We choose a variety of physical quantities that are meaningful in the simulation of canyon upwelling and examine their representation in each model run for evidence of convergence of the solution. [29] First, we look at the alongshore velocity on the opposite side of the domain from the canyon. Far away from the area of flow disturbed by the canyon, the alongshore velocity structure is fairly similar at all three model resolutions with a radially increasing velocity in the abyss, a velocity maximum in a small jet just offshore of the shelf break, and a broad region of roughly constant velocity over the shelf (Figure 8). Viscous boundary layers are apparent close to the topography. As the model resolution increases, the maximum speed within the jet slightly increases while the velocity over the shelf decreases. These changes are consistent with the two‐dimensional model results that show Ekman transport increasing with resolution. More upslope Ekman transport over the shelf results in greater reduction of surface pressure gradients, which in turn reduces the alongshore barotropic geostrophic velocity. At the same time, the downslope return flow above the Ekman layer results in a depression of isopycnals at the edge of the shelf break,  Figure 8. Azimuthal velocity in cm s−1 on the opposite side of the tank from the canyon at (a) 5 mm horizontal × 2.8 mm vertical, (b) 3 mm × 1 mm, and (c) 1.5 mm × 0.5 mm grid resolutions.  6 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  Figure 9. Vertical displacement in millimeters of fluid between the start and end of each model run along the canyon axis at (a) 5 mm horizontal × 2.8 mm vertical, (b) 3 mm × 1 mm, and (c) 1.5 mm × 0.5 mm grid resolutions. lution and is not apparent at all in the 1.5 × 0.5 mm run. Overall, it appears that the velocity structure over the abyss and continental slope regions is well resolved in all the models and is likely within 5% of being resolved over the shelf, consistent with the results from the two‐dimensional model. This is unsurprising as far away from the canyon the low‐resolution laminar flow should behave much like the two‐dimensional model results. [31] Next, we turn to measures of the upwelling within the canyon and examine the depth change of water upwelled within the canyon over the 30 s of the model run. Assuming that diffusion of density is negligible implies that density changes over the model run are solely due to advection of the initial model density structure. Since the initial density is  C05008  stably stratified and varies only with depth, every density is associated with a single initial depth. By subtracting the depth of each point by the initial depth associated with the density at that point, we get the depth change of the water at that model point. [32] All models show a general pattern of strong upwelling within the canyon with a maximum at the canyon head, smaller upwelling in the Ekman layer over the shelf, and a slight downwelling elsewhere (Figure 9). In the 5 × 2.8 mm model (Figure 9a), the maximum depth change at the canyon head is ∼9 mm, while in the higher‐resolution runs (Figures 9b and 9c) the maximum upwelling is ∼11 mm. In addition, the higher‐resolution runs show a larger region of downwelling over the shelf than the 5 × 2.8 mm model. In contrast to these differences, the 3 × 1 mm and 1.5 × 0.5 mm resolution model upwelling fields appear nearly identical, suggesting that the canyon upwelling solution is fully converged. [33] For a second measure of the canyon upwelling, we look at the total volume of water upwelled through the top of the canyon over the model run. We calculate this by summing the vertical and radial volume fluxes through the top of the canyon, defined as the surface of the topography if there was no canyon. The fluxes are summed only for radii >27.9 cm (the shelf break radius) so as to only include upwelling out the top of the canyon. Fluxes are sampled once per second and are assumed to be constant over that second. We perform this calculation on the three main model runs and additionally include the results from a run performed at 2 × 0.67 mm resolution (Table 1). The total fluid volume upwelled increases by ∼19% between the 5 × 2.8 mm and 3 × 1 mm runs but in contrast is nearly identical between the three highest‐resolution runs, with a variation of <1% in the volume of fluid upwelled. Again, this suggests that the 1.5 × 0.5 mm model’s canyon upwelling solution is fully converged. [34] The structure of the vertical velocity across the canyon shows the structure of standing internal waves generated by the flow over the canyon and serves as a test of the convergence of the MITgcm’s nonhydrostatic dynamics (Figure 10). A series of areas of alternating upward and downward motion are apparent as the flow moves over the canyon. The 5 × 2.8 mm run tends to overestimate the magnitude of these velocities (Figure 10a) relative to the higher‐ resolution runs. The magnitudes and overall patterns of the vertical velocity in the 3 × 1 mm and 1.5 × 0.5 mm runs are very similar, but many small‐scale differences are still apparent, especially where the flow enters the canyon on the upstream side. Nevertheless, the overall changes between the 3 × 1 mm and 1.5 × 0.5 mm runs are small and suggest that the 3 × 1 mm model solution has nearly converged. [35] Next, we examine the density field at the end of the model runs. The model density is interesting for a few reaTable 1. Net Fluid Volume Upwelled Onto the Shelf From Within the Canyon During Each MITgcm Model Run Resolution (mm)  Upwelling (cm3)  5 × 2.8 3×1 2 × 0.67 1.5 × 0.5  16.63 20.43 20.57 20.48  7 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  Figure 10. Vertical velocity in mm s−1 across the canyon mouth at a radius of 30 cm at the end of the model run at (a) 5 mm horizontal × 2.8 mm vertical, (b) 3 mm × 1 mm, and (c) 1.5 mm × 0.5 mm grid resolutions. Figures 10a–10c are plotted looking toward the head of the canyon at 30 cm radius; mean alongshore flow is from left to right. sons. First, it is the result of the integrated three‐dimensional density advection over the model run and thus is sensitive to the details of the time evolution of the system. Second, it lets us compare the density structure from this z coordinate model with the density from a previous s coordinate model examined by Allen et al. [2003]. Third, Hickey [1997] measured temperature structure over Astoria Canyon, and these measurements allow us to compare the model density structure to ocean data. [36] During upwelling events in Astoria Canyon, Hickey [1997] found that the density field within the canyon formed a wedge‐shaped structure (Figure 11a), with fluid  C05008  layers near the top of the canyon thinning on the upstream side and growing thicker on the downstream side as the flow pushed fluid downstream. A similar effect occurs in the model simulations (Figures 11b–11e), where we consider the “wedge” to be the fluid bounded by the 21 s and 24 s surfaces. [37] Models of s coordinates are generally considered more suitable for simulating flows over topographic features than z coordinate models. Allen et al. [2003] performed a study similar to the present one using the s coordinate Rutgers University model (SCRUM) [Song and Haidvogel, 1994] and found systematic problems with SCRUM’s representation of the laboratory model. To compare our current results with this previous research, we performed a simulation with SCRUM using our current model’s topography and forcing. All other aspects of the model were configured as given by Allen et al. [2003]. The SCRUM model run was configured with 16 s levels, which results in a vertical resolution of ∼4 mm within the canyon, making the SCRUM resolution slightly coarser than the resolution of the 5 × 2.8 mm MITgcm simulation. [38] In both the SCRUM simulation (Figure 11b) and the 5 × 2.8 mm MITgcm simulation (Figure 11c) the density wedge is not well defined, with neither model doing a particularly good job representing the density wedge. This suggests that terrain coordinate models are not necessarily better than z coordinate models in representing this type of flow over steep topography. As we increase the model resolution (Figures 11d and 11e), the thickness difference between the upstream and downstream sides of the wedge consistently increases, and a well‐defined wedge structure emerges. The change in the downstream wedge thickness between the 3 × 1 mm and 1.5 × 0.5 mm runs is small, but the upstream thickness changes significantly between the runs. This thickness, however, is likely limited by the vertical resolution of the two models. [39] Finally, we compare the representation of potential vorticity (PV) in the model as resolution changes. PV conservation represents a strong constraint on geostrophic flows and furthermore depends upon gradients in both velocity and density. For these reasons, PV should be more sensitive to changes in the flow structure than the velocity or density fields themselves. Because the model topography is much steeper than topography in the ocean, the horizontal vorticity in the model is of comparable magnitude to the vertical vorticity, and thus, we must use the full Ertel potential vorticity in our comparison: PV ¼  ð2W þ r  uÞ Á r:   ð4Þ  [40] To examine the representation of PV in the model, we calculate the difference between the initial and final PV in the model and then subtract from this the PV change due to the models’ azimuthal body force. This results in the PV changes due to advection and molecular diffusion (Figure 12). A three‐layer structure is apparent in the background PV changes, with increased PV in the canyon, decreased PV just over the shelf, and increased PV near the surface. The initial PV structure of the model is horizontally uniform and negative, consisting of near‐zero PV near the surface and bottom  8 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  Figure 11. Density structure (a) across the canyon mouth of Astoria Canyon during an upwelling event (based on Figure 3 from Hickey [1997]) and simulated by the (b) SCRUM model run and the MITgcm model runs at (c) 5 mm horizontal × 2.8 mm vertical, (d) 3 mm × 1 mm, and (e) 1.5 mm × 0.5 mm grid resolutions. The 26.6 st and 26.8 st surfaces in Figure 11a and the 21 s and 24 s surfaces in Figures 11b– 11e are highlighted to emphasize the “density wedge” structure that forms in an upwelling canyon. Figures 11a–11e are plotted looking toward the head of the canyon at 30 cm radius; mean alongshore flow is from left to right. (due to the diffusion of the surface and bottom density gradients during the experiment setup) with a PV minimum near the shelf break. The PV changes in the canyon and near the surface can be explained as vertical advection of this background gradient as higher PV is advected upward in the canyon and downward near the surface. The middle band of decreased PV is the result of friction with the shelf topography. [41] Significant small‐scale changes exist between the PV representations at all model resolutions. These changes  mostly represent frictional effects and so can be thought of as a map of the influence of the boundary layer on the flow. At 5 × 2.8 mm resolution the PV change has little small‐ scale structure, though there are some isolated hot spots of PV increase due to viscous effects near the walls of the canyon. At 3 × 1 mm resolution the PV change shows more viscous forcing, with regions of increased PV developing on the upstream and downstream flanks of the upper canyon. Additionally, a dipole structure of PV decrease lying over PV increase develops at the upstream edge of the canyon.  9 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  converged at 3 × 1 mm resolution; it may have converged by 1.5 × 0.5 mm resolution, but it is impossible to tell this from the model runs we have performed. Thus, small‐scale details of the simulation of advection and viscosity still may not be adequately represented.  5. Comparison With Lab Experiment  Figure 12. Ertel potential vorticity change between the start and end of each model run at (a) 5 mm horizontal × 2.8 mm vertical, (b) 3 mm × 1 mm, and (c) 1.5 mm × 0.5 mm grid resolutions. The effect of the model forcing on the potential vorticity has been removed. Figures 12a–12c are plotted looking toward the head of the canyon at 30 cm radius; mean alongshore flow is from left to right. This dipole is the PV signature of the density wedge structure and results from viscosity acting on vorticity forced by density surface squashing and stretching in the density wedge. At 1.5 × 0.5 mm resolution the overall pattern remains the same, but regions with fine PV structure, such as near the canyon walls, show more definition and stronger gradients. The PV dipole at the upstream edge of the canyon is larger in magnitude and spatial extent, nearly reaching the downstream canyon wall. Clearly, the model’s simulation of PV has not  [42] Over sections 3 and 4 of this paper we have examined the convergence of the model solution; however, we have not examined the accuracy of the solution that the model is converging to. To do this, we compare the horizontal flow field of the model with the flow field generated by a laboratory tank experiment performed by Bowie [2006]. [43] Bowie utilized neutral density wax particles to track the flow field of an upwelling flow in the laboratory tank. The particles were illuminated by horizontal sheets of colored light that allowed estimation of the pellets’ vertical locations: blue illuminated particles were between 14.5 and 17 mm depth, green particles were 17–18.6 mm, yellow particles were 18.6–19.5 mm, white particles were 19.5– 22.8 mm, and red particles were 22.8–25.5 mm. A video camera mounted above the tank was then used to record the trajectories of particles moving in and around the canyon. The density gradient was slightly weaker in Bowie’s laboratory experiments compared to our previous model runs (N = 2.0 s−1 at the shelf break), but all other physical parameters were the same. We performed a model run at 1.5 × 0.5 mm resolution (the highest resolution in section 4) that duplicated Bowie’s setup and recorded the model velocity field at 0.5 s increments. Numerical particles were then placed in the saved model velocity fields, and their trajectories were tracked from 28 to 30 s after the start of the model run, using a fourth‐order Runge‐Kutta algorithm. [44] Comparing the particle tracks from the laboratory experiment with the model experiment (Figure 13) shows that the model does a reasonable job of reproducing the laboratory tank’s flow field. In the abyssal area away from the shelf, particle tracks are smooth circular arcs indicating solid body rotation. Flow entering the canyon from the shelf turns left and flows up the canyon. Near the head of the canyon, the flow is directly up the canyon axis onto the shelf, with stronger flow on the downstream side. Some flow from the canyon onto the downstream shelf is also apparent. [45] However, there are two features of the model particle tracks that do not agree with the laboratory experiment. First, the model tracks show many more midlevel particles than are apparent in the laboratory tracks. This is due to oversaturation of the camera CCD during the laboratory experiment, resulting in many particles being misidentified as white [Bowie, 2006]. Second, the laboratory experiment shows particles on the upstream flank of the canyon mouth moving offshore into the abyss across isobaths before turning back onshore in a wide arc and flowing up the canyon. This offshore flow is a robust feature of the laboratory experiments and also appears in field data from Astoria Canyon [Hickey, 1997]. In contrast, the model particles tend to flow directly off the upstream flank of the shelf and then make a hard turn when they reach the downstream canyon wall. Thus, while the model solution is identical to the laboratory experiment in  10 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  C05008  Figure 13. Particle tracks from 28 to 30 s from (a, b) laboratory experiments and (c, d) MITgcm model output. Figures 13a and 13c show particles below the shelf break depth, and Figures 13b and 13d show particles above the shelf break depth (1.7 cm). Contours represent model topography, with the thicker contour indicating model shelf break depth. Particles are shown every 0.3 s and become slightly larger with time. most respects, it fails to simulate the flow on the upstream canyon flank realistically.  6. Discussion [46] We believe the failure of the model to simulate the flow on the upstream flank of the canyon is due to inadequate representation of flow separation. Temperature data from Astoria Canyon [Hickey, 1997] suggest that water drops slightly into the canyon as it leaves the shelf (Figure 11a, shown by the depth of the 26.6 st isopycnal). This drop results in stretching of the water column as it leaves the shelf, which would tend to induce a counterclockwise, positive relative vorticity. Positive relative vorticity would create an offshore flow near the canyon wall and would increase the onshore flow within the canyon. In contrast, the model flow hardly drops at all as it crosses the canyon (Figure 11e, the 21 s isopycnal). This suggests that the model flow is separating from the upstream flank of the canyon too early.  [47] We can estimate the amount of stretching needed to produce the observed flow by estimating the vorticity of the flow from Figure 13a. Particles move ∼5 mm offshore over the 2 s displayed in Figure 13 as they leave the upstream flank of the canyon and move a similar distance back onshore as they near the canyon axis; thus, we take the radial velocity change in this flow to be 0.25 − (−0.25) = 0.5 cm s−1. The distance from the upstream canyon edge to the canyon axis is ∼3 cm. Thus, the vorticity z of the flow is ¼  @v ¼ 0:17 sÀ1 : @x  ð5Þ  [48] Using the potential vorticity of a homogeneous layer of thickness h, (z + f)/h, we can estimate the amount of vortex stretching needed to generate this much vorticity. Before leaving the upstream flank of the canyon the flow has a vorticity of 0, so  11 of 13  f þf ¼ ; hi hf  ð6Þ  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  where hi is the initial thickness of the layer and hf is the thickness of the layer after stretching. Then we write hf  ¼ ¼ 1:32: hi f þ 1  ð7Þ  Examining Figure 11a, a 30% stretching is consistent with the stretching observed in Astoria Canyon. [49] One significant dynamical difference between our laboratory canyon and an ocean canyon is turbulence. True coastal canyon systems are far from laminar, with large turbulent flows that result in eddy diffusivities that vary in space and time and are of comparable magnitude for both viscosity and material properties of the fluid. Laboratory vertical and horizontal viscosities, once scaled to the size of the real ocean, are the same order of magnitude as measured ocean values, but salinity diffusivity will be 3 orders of magnitude higher. Bottom boundary layer depths are ∼20 m in the coastal ocean and thus are the same relative size as the Ekman layer in the laboratory tank, but the bottom turbulence will rise and fall in intensity with the currents and the tides, resulting in a time‐varying bottom layer depth. Despite this deficiency, our lab results compare favorably with field observations, suggesting that these turbulent effects are small.  7. Conclusions [50] This paper examined the convergence of a numerical model simulation of an idealized laboratory model of an upwelling continental shelf canyon system with increasing resolution. A series of increasingly complex models were used to examine the conditions necessary to resolve flow features. A simple 1‐D model of bottom Ekman layer transport indicates that achieving a transport within 1% of the true transport value requires a vertical resolution of ∼0.2d, where d is the Ekman layer thickness. [51] A cylindrically symmetric, 2‐D, z level model shows that both horizontal and vertical resolution are important to resolution of the Ekman layer over sloped topography, while only vertical resolution is important over flat topography. This indicates that two concerns are important: achieving sufficient vertical resolution to resolve the velocity structure of the Ekman layer and sufficient horizontal resolution to achieve a smooth representation of the shelf slope, without abrupt steps in the topography. The partial cell parametrization present in the MITgcm helps the representation of topography but does not allow for improvements in the Ekman layer velocity structure and thus does not resolve problems with simulating Ekman layer transports. [52] Three full 3‐D MITgcm runs were performed with horizontal and vertical resolutions of 5 × 2.8 mm, 3 × 1 mm, and 1.5 × 0.5 mm. Significant changes in the model velocity, total canyon upwelling, and density advection occur between the 5 × 2.8 mm and 3 × 1 mm resolution runs, but there is little change in these fields between the 3 × 1 mm and 1.5 × 0.5 mm runs. Representation of potential vorticity, on the other hand, shows significant differences between the 3 × 1 mm and 1.5 × 0.5 mm runs and does not appear to have converged. Comparison of the 1.5 × 0.5 mm resolution model with an identically configured laboratory experiment shows that we are able to simulate most features of the laboratory tank flow field. However, the model un-  C05008  derestimates the amount of offshore flow on the upstream flank of the canyon. This problem is likely due to the flow separating from the canyon wall too early. For problems that are not dependent on the details of the flow field on the upstream flank of the canyon, such as advection over the shelf or total upwelling through the canyon, modern computational clusters appear to have sufficient power to adequately represent an upwelling continental shelf canyon. [53] Acknowledgments. Funding for this research was provided to the second author via a NSERC Discovery grant. We would like to thank Mike Dinniman and John Klinck from Old Dominion University for providing SCRUM model output and the two anonymous peer reviewers whose comments significantly improved aspects of this manuscript. This research has been enabled by the use of WestGrid computing resources, which are funded in part by the Canada Foundation for Innovation, Alberta Innovation and Science, BC Advanced Education, and the participating research institutions. WestGrid equipment is provided by IBM, Hewlett Packard, and SGI. Figures 1–13 were generated using the matplotlib library in the Python programming language.  References Adcroft, A., C. Hill, and J. Marshall (1997), Representation of topography by shaved cells in a height coordinate ocean model, Mon. Weather Rev., 125, 2293–2315. Allen, S. E. (2000), On subinertial flow in submarine canyons: Effect of geometry, J. Geophys. Res., 105, 1285–1297. Allen, S. E. (2004), Restrictions on deep flow across the shelf‐break and the role of submarine canyons in facilitating such flow, Surv. Geophys., 25, 221–247. 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, doi:10.1139/cjfas-58-4-671. 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(C1), 3003, doi:10.1029/2001JC000978. Bosley, K. L., J. W. Lavelle, R. D. Brodeur, W. W. Wakefield, R. L. Emmett, E. T. Baker, and K. M. Rehmke (2004), Biological and physical processes in and around Astoria submarine canyon, Oregon, USA, J. Mar. Syst., 50, 21–37. Bowie, A. W. (2006), Upwelling in short submarine canyons: A physical study, M.S. thesis, Univ. of B. C., Vancouver, B. C., Canada. Boyer, D. L., X. Zhang, and N. Pérenne (2000), Laboratory observations of rotating, stratified flow in the vicinity of a submarine canyon, Dyn. Atmos. Oceans, 31, 47–72. Boyer, D. L., D. B. Haidvogel, and N. Pérenne (2004), Laboratory‐numerical model comparisons of canyon flows: A parameter study, J. Phys. Oceanogr., 34, 1588–1609. Boyer, D. L., J. Sommeria, A. S. Mitrovic, V. K. C. Pakala, S. A. Smirnov, and D. Etling (2006), The effects of boundary turbulence on canyon flows forced by periodic along‐shelf currents, J. Phys. Oceanogr., 26, 813–826. Freeland, H., and K. Denman (1982), A topographically controlled upwelling center off Vancouver Island, J. Mar. Res., 40, 1069–1093. Haidvogel, D. B. (2005), Cross‐shelf exchange driven by oscillatory barotropic currents at an idealized coastal canyon, J. Phys. Oceanogr., 35, 1055–1067. Hickey, B. M. (1997), The response of a steep‐sided, narrow canyon to time‐variable wind forcing, J. Phys. Oceanogr., 27, 697–726. Hill, D. F. (2002), General density gradients in general domains: The “two‐ tank” method revisited, Exp. Fluids, 32, 434–440. Jaramillo, S. (2005), Numerical simulation of flow in a laboratory tank using a z‐coordinate numerical model, M.S. thesis, Univ. of B. C., Vancouver, B. C., Canada. Klinck, J. M. (1996), Circulation near submarine canyons: A modeling study, J. Geophys. Res., 101, 1211–1223. Kundu, P. K. (1990), Fluid Mechanics, 638 pp., Academic, San Diego, Calif. Kunze, E., and S. G. L. Smith (2004), The role of small‐scale topography in turbulent mixing of the global ocean, Oceanography, 17(1), 55–64.  12 of 13  C05008  DAWE AND ALLEN: SOLUTION CONVERGENCE OF FLOW OVER CANYON  Ladd, C., P. Stabeno, and E. D. Cokelet (2005), A note on cross‐shelf exchange in the northern Gulf of Alaska, Deep Sea Res., Part II, 52, 667–679. Laurent, L. S., S. Stringer, C. Garrett, and D. Perrault‐Joncas (2003), The generation of internal tides at abrupt topography, Deep Sea Res., Part I, 50, 987–1003. MacCready, P., and P. B. Rhines (1993), Slippery bottom boundary layers on a slope, J. Phys. Oceanogr., 23, 5–22. Marshall, J., A. Adcroft, C. Hill, L. Perelman, and C. Heisey (1997), A finite‐volume, incompressible Navier Stokes model for studies of the ocean on parallel computers, J. Geophys. Res., 102, 5753–5766. Mirshak, R., and S. E. Allen (2005), Spin‐up and the effect of a submarine canyon: Applications to upwelling in Astoria Canyon, J. Geophys. Res., 110, C02013, doi:10.1029/2004JC002578. Oster, G. (1965), Density gradients, Sci. Am., 213(2), 70–76. Palanques, A., et al. (2005), General patterns of circulation, sediment fluxes and ecology of the Palamós (La Fonera) submarine canyon, northwestern Mediterranean, Prog. Oceanogr., 66, 89–119. Pérenne, N., J. Verron, D. Renouard, D. L. Boyer, and X. Zhang (1997), Rectified barotropic flow over a submarine canyon, J. Phys. Oceanogr., 27, 1868–1893. Pérenne, N., D. B. Haidvogel, and D. L. Boyer (2001a), Laboratory‐numerical model comparisons of flow over a coastal canyon, J. Atmos. Oceanic Technol., 18, 235–255. Pérenne, N., J. W. Lavelle, D. C. Smith IV, and D. L. Boyer (2001b), Impulsively started flow in a submarine canyon: Comparison of results from laboratory and numerical models, J. Atmos. Oceanic Technol., 18, 1698–1718.  C05008  Sobarzo, M., M. Figueroa, and L. Djurfeldt (2001), Upwelling of subsurface water into the rim of the Biobío submarine canyon as a response to surface winds, Cont. Shelf Res., 21, 279–299. Song, Y., and D. B. Haidvogel (1994), A semi‐implicit ocean circulation model using a generalized topography‐following coordinate system, J. Comput. Phys., 115, 228–244. Thorpe, S. A. (1987), Current and temperature variability on the continental slope, Philos. Trans. R. Soc. London, Ser. A, 323(1574), 471–517. Wang, D.‐P., L.‐Y. Oey, T. Ezer, and P. Hamilton (2003), Near‐surface currents in DeSoto Canyon (1997–99): Comparison of current meters, satellite observation, and model simulation, J. Phys. Oceanogr., 33, 313–326. Waterhouse, A. F. (2006), A physical study of upwelling flow dynamics in long canyons, M.S. thesis, Univ. of B. C., Vancouver, B. C., Canada. Waterhouse, A. F., S. E. Allen, and A. W. Bowie (2009), Upwelling flow dynamics in long canyons at low Rossby number, J. Geophys. Res., 114, C05004, doi:10.1029/2008JC004956. Williams, W. J., E. C. Carmack, K. Shimada, H. Melling, K. Aagaard, R. W. Macdonald, and R. G. Ingram (2006), Joint effects of wind and ice motion in forcing upwelling in Mackenzie Trough, Beaufort Sea, Cont. Shelf Res., 26, 2352–2366. S. E. Allen and J. T. Dawe, Department of Earth and Ocean Sciences, University of British Columbia, 6339 Stores Rd., Vancouver, BC V6T 1Z4, Canada. (sallen@eos.ubc.ca; jdawe@eos.ubc.ca)  13 of 13  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items