JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 108, NO. C1, 3003, doi:10.1029/2001JC000978, 2003 On vertical advection truncation errors in terrain-following numerical models: Comparison to a laboratory model for upwelling over submarine canyons S. E. Allen,1 M. S. Dinniman,2 J. M. Klinck,2 D. D. Gorby1 A. J. Hewett,1 and B. M. Hickey3 Received 14 May 2001; revised 27 December 2001; accepted 10 April 2002; published 3 January 2003. [1] Submarine canyons which indent the continental shelf are frequently regions of steep (up to 45°), three-dimensional topography. Recent observations have delineated the flow over several submarine canyons during 2–4 day long upwelling episodes. Thus upwelling episodes over submarine canyons provide an excellent flow regime for evaluating numerical and physical models. Here we compare a physical and numerical model simulation of an upwelling event over a simplified submarine canyon. The numerical model being evaluated is a version of the S-Coordinate Rutgers University Model (SCRUM). Careful matching between the models is necessary for a stringent comparison. Results show a poor comparison for the homogeneous case due to nonhydrostatic effects in the laboratory model. Results for the stratified case are better but show a systematic difference between the numerical results and laboratory results. This difference is shown not to be due to nonhydrostatic effects. Rather, the difference is due to truncation errors in the calculation of the vertical advection of density in the numerical model. The calculation is inaccurate due to the terrain-following coordinates combined with a strong vertical gradient in density, vertical shear in the horizontal velocity and topography with strong INDEX TERMS: 4255 Oceanography: General: Numerical modeling; 3337 Meteorology and curvature. Atmospheric Dynamics: Numerical modeling and data assimilation; 4520 Oceanography: Physical: Eddies and mesoscale processes; KEYWORDS: model, sigma, error, topography, upwelling, canyon Citation: Allen, S. E., M. S. Dinniman, J. M. Klinck, D. D. Gorby, A. J. Hewett, and B. M. Hickey, 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, 2003. 1. Introduction [2] Both laboratory and numerical models have been used to study processes in coastal ocean circulation. Recent studies have compared test simulations from different numerical models with the disquieting result that the solutions differ not only quantitatively but, in some cases, qualitatively [Haidvogel and Beckmann, 1998]. Given the wide variety of processes and parameterizations in these models, the cause(s) of these differences have not been identified. [3] The traditional method for verifying numerical models is to compare to an analytical solution. However, only a limited number of analytic solutions exist for stratified flow over steep topography. An alternative evaluation method is to compare results from different models. However dis1 Oceanography, Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada. 2 Center for Coastal Physical Oceanography, Old Dominion University, Norfolk, Virginia, USA. 3 School of Oceanography, University of Washington, Seattle, Washington, USA. Copyright 2003 by the American Geophysical Union. 0148-0227/03/2001JC000978$09.00 agreements between the models frequently tend to raise questions instead of answering them. [4] We chose to compare results from numerical simulations to those from a physical (laboratory) model of both homogeneous and stratified fluids flowing over a submarine canyon on a continental shelf. This problem activates a number of physical processes (discussed below) making it a good test case [Haidvogel and Beckmann, 1998]. The circulation we consider starts impulsively; the adjustment to steady state is compared between the two models. [5] Laboratory models are thought to be better analogs for geophysical scale flow than numerical models since a wider range of length scales are active. Regional circulation is produced by forces acting on scales from 100 km to 1 mm (the turbulence dissipation scale), a range of 8 orders of magnitude. A laboratory model has scales from 1 m to 1 mm, a range of three orders of magnitude. Numerical models represent scales from the domain size to the grid spacing, the ratio of which is typically 100 (the number of grid points in one direction), a range of two orders of magnitude. One might argue that sub-grid-scale parameterizations represent scales smaller than the grid, but these formulations range in quality from approximately correct to merely plausible but unverified [Gent et al., 1995]. These problems with formulations 3-1 3-2 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 1. A schematic of flow over a canyon. The near surface and deep layer are shown explicitly. The deep shelf flow and the upwelling current together form the active layer. The shelf-break depth is approximately 200 m in the field and is 2.2 cm in the lab. The near surface flow extends down about 150 m from the surface in the field and 1.6 cm from the surface in the lab. The near surface flow-passes over the canyon with little deflection. In the active layer the flow is strongly constrained by the topography. Flow in the deep layer turns cyclonically. [after Allen et al., 2001] are added to issues of truncation errors in numerical schemes, domain distortion (e.g., bottom-following coordinates), and other parameterizations (e.g., equation of state, estimates of bottom or surface stress, etc.). Therefore, laboratory models should represent good solutions (assuming appropriate laboratory technique) which can be used as a stringent test for numerical models. [6] The extrapolation of laboratory analogs to a geophysical scale system involves measures of similarity. Similarity theory applied to fluid dynamics has yielded rich results over the centuries [Batchelor, 1967]. Flows with the same values of nondimensional ratios, which compare sizes of physical processes, will be identical other than having different scales. One need not match all nondimensional numbers in such a comparison, only the ‘‘important’’ ones. In many cases, ‘‘important’’ is easy to determine; in others it is open to interpretation. For advection over a canyon the important nondimensional numbers are a Rossby Number, a Burger Number, a Froude Number and a geometric number describing the shape of the canyon [Allen et al., 1999]. [7] First we present an overview of the flow regime under consideration. Next the two flow representations (laboratory and numerical) are presented and their results are compared for both homogeneous fluid and stratified fluids. The discussion includes an estimate of the effect of nonhydrostatic terms in the laboratory results and an estimate of the effect of errors in density advection on the numerical results. [8] Submarine canyons are common features of the continental shelf. A typical canyon cuts from the continental slope into the continental shelf about halfway to the coast. Data sets including off-axis information are available over Astoria Canyon [Hickey, 1997], Lydonia Canyon [Butman, 1983], Barkley Canyon [Allen et al., 2001] and Carson Canyon [Kinsella et al., 1987]. Submarine canyons are regions of enhanced upwelling during upwelling favourable conditions and thus are important in evaluating cross-shelf exchange in coastal regions. The problem used to compare the two models is a short duration upwelling event. Previous numerical studies have shown qualitative agreement with the observations [Klinck, 1996; Allen, 1996; She and Klinck, 2000]. [9] After acceleration of the shelf break current (in the field due to the wind or shelf waves) a flow pattern forms over the canyon with strong variation in the vertical (Figure 1). In particular: level 1, the near surface current is not affected or is only weakly affected by the canyon and travels straight over the canyon; level 2, flow just above the canyon rim moves over the canyon, descends into the canyon, turns toward the canyon head and crosses the downstream canyon rim near the head; level 3, slightly deeper flow over the continental shelf bends into the canyon near the mouth, traverses up the canyon and crosses the downstream canyon rim near the head; and level 4, deep flow within the canyon turns cyclonically. Flows numbered 2 and 3 above form the active layer (Figure 1). Observations [Hickey, 1997; Allen et al., 2001] show a trapped eddy over the ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Table 1. Scales for Astoria Canyon [Hickey, 1997] and the Laboratory and Numerical Models Considered Here Number Hs U f L W r N Wsb Ro Fr Bu Ge Astoria Lab Numerical 150 m 0.2 m sÀ1 1.05 Â 10À4 sÀ1 21.8 km 8.9 km 4.5 km 7.5 Â 10À3 sÀ1 15.7 km 0.42 0.18 0.49 0.57 2.2 cm 1.2 cm sÀ1 0.52 sÀ1 8 cm 2.4 cm 1.4 cm 2.2 sÀ1 6.9 cm 1.6 0.25 1.2 0.35 2.2 cm 1.2 cm sÀ1 0.052 sÀ1 80 cm 24 cm 14 cm 2.2 sÀ1 69 cm 1.6 0.25 1.2 0.35 canyon (in flow numbered 2 above). The presence or absence of the eddy is due to the amount of stretching near the canyon rim [Allen et al., 2001]. [10] The parameters governing the advective flow over a submarine canyon during an upwelling event are [Allen et al., 1999] (1) a Rossby number: Ro = U/fr where U is the incident velocity at the half length of the canyon and at the depth of the shelf break Hs, f is the Coriolis parameter and r is the radius of curvature of the shelf break isobath; (2) a Froude number: U/NHs where N is the buoyancy frequency at the depth of the shelf break; (3) a Burger number: (N Hs)/ fL. The length of the canyon L is measured, at shelf break depth, from the shelf break to the head (last isobath deflected significantly by the canyon); (4) a Geometric parameter Ge = W/Wsb where Wsb is the width of the canyon at the shelf break. The width of the canyon W is measured at half length. The geometric scales, strength of stratification and rotation rate for the laboratory model were chosen to give values of the nondimensional parameters similar to those for Astoria Canyon (Table 1). 3-3 world. The second major laboratory limitation is the necessity of exaggerating the vertical scale. The total vertical scale must be large enough to reduce viscous effects; yet the total horizontal scale is limited by the size of the apparatus. In some cases (here, in the homogeneous case), vertical exaggeration leads to nonhydrostatic effects that would not be expected in the real world. Although the eventual goal of our research is to model the real world, here where we are comparing two models, the numerical model is configured to use uniform molecular viscosity. [14] The laboratory experiments were performed on a 1-m diameter, computer-controlled rare earth motor-driven rotating table at the University of British Columbia. The outside tank is a cylindrical Plexiglas tank of radius 50 cm. It was fitted with a Plexiglas insert to give a generic coastal bathymetry consisting of a flat abyssal plane, a steep continental slope (slope 45° in laboratory) and a continental shelf (slope 5° in laboratory). The continental shelf and slope were fitted to the outside of the tank with the abyssal plane in the center. Without the canyon, the topography has cylindrical symmetry. A 22° slice was cut out of the generic bathymetry to allow insertion of a canyon shape (Figure 2). The canyon was made of layers of 1 mm thick sheets of polystyrene. The layers were joined using plastic cement and the canyon model was covered with automotive filler and sanded to produce a smooth finish and so join smoothly with the generic shelf topography. The canyon was painted flat black to allow flow visualization. The exact bathymetry 2. Model Details [11] The test problem is largely determined by the laboratory setup, in which the flow is started with a rapid change in rotation rate of the table. The numerical model (SCRUM Version 3.1, a terrain-following coordinate model, Song and Haidvogel [1994]) is configured as closely as possible to the laboratory setup with an initial, impulse-like body force starting the flow. The body force is applied over one rotation period to match the time over which the rotation rate of the table is changed. [12] Two sets of lab experiments are considered, one with homogeneous fluid, and the other with stratified fluid. The laboratory and numerical geometries using a homogeneous fluid are identical. However, in the stratified case the horizontal dimension was increased by a factor of 10 to reduce the bottom slope and allow the simulation to proceed. Appropriate adjustments were made to other parameters (detailed below) to preserve the relevant nondimensional ratios. 2.1. Laboratory Scale Model [13] Physical or laboratory models have their specific limitations. Generally, the scale of the model implies moderate Reynolds number (order 1000). Either the viscous effects are minimized or the molecular viscosity in the tank is assumed to model a uniform eddy viscosity in the real Figure 2. Schematic bathymetric map of the topography in the tank. Contours are 1 cm apart. The vertical scale is exaggerated compared to the horizontal scale by a factor of 10 (radius of tank is 50 cm). The ‘‘shelf’’ slope is 5° and the ‘‘slope’’ slope is 45°. The numerical geometry is similar except: (1) for the stratified simulations the horizontal scale is increased by a factor of 10 (so there is no vertical exaggeration in the numerical model) and (2) only 60% of the tank is simulated. 3-4 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 3. Density profile in the tank and numerical simulation. The initial density profile (solid) is calculated assuming no extra mixing in the filling process. The expected diffusion after 1 hour of spin-up is shown by the dashed line. The numerical simulation used the solid profile as its initial condition. was determined by digitizing photographs. Horizontal light sheets were used to illuminate the canyon surface. The line of illumination marked an isobath with depth dependent on the height of the horizontal light sheet. Each picture produced one bathymetric contour and the contours were taken 1 cm apart. This bathymetry was used in the numerical model. [15] For the homogeneous experiments, the tank was filled with fresh water and then accelerated to its starting rotation speed ( f = 0.4 sÀ1). The table was allowed to rotate until the fluid reached solid body rotation and was stationary with respect to the tank (about an hour). [16] For the stratified case the table was first accelerated to its starting rotation speed ( f = 0.4 sÀ1). Then the tank was filled by gravity feed through a tube coming down from the ceiling along the axis of rotation. The tank was filled from the bottom using the Osler method of stratification [Osler, 1965]. In a tank with straight sidewalls the Osler method results in a linear stratification. Because of the geometry of the shelf slope topography, here the stratification is stronger near the top (Figure 3). The table was allowed to rotate until the fluid reached solid body rotation and was stationary with respect to the tank (several hours). The tank stratification was replaced each day to reduce the formation of ‘‘mixed’’ layers at the top and bottom of the tank. [17] The flow was forced by accelerating the table from f = 0.4 sÀ1 to f = 0.52 sÀ1 linearly over 27.3 s. This period is about one revolution of the tank and is comparable to a wind event of one day in the real world. The change in rotation speed causes the bathymetry to accelerate relative to the water, which due to inertia continues to rotate at the slower rate. Relative to the tank, the water moves with shallow water to its left. This forcing results in flow at the shelf break of about 1.2 cm sÀ1. [18] The flow in the vicinity of the canyon was visualized using neutrally buoyant particles (pliolite) which were added upstream of the canyon during the acceleration phase of the experiment. The tank was lit by a slide projector with a bicolored slit. The projector was mounted above the tank attached to the table (i.e., it rotated with the table). The light bounced off a 45° mirror set in the abyssal plain near the axis of rotation. Two horizontal sheets of light of 1 cm thickness resulted. The green light illuminated depths of 12 to 22 mm and the red light illuminated depths of 22 to 32 mm. A Hi-8 video camera mounted above the tank (rotating with the table) recorded the movement of the particles. [19] The flow was analyzed for a 5 s interval, 2.7 s after the acceleration of the tank was complete. Using frame capture software, frames from the video were captured onto computer at 0.5 s intervals. Particles were traced manually from frame to frame. The experiment was repeated 5 – 15 times. Composite photos of all the frames for the homogeneous case and the stratified case were generated and digitized versions of these composite figures are used for model comparison. [20] Two major sources of error occur in the laboratory. First, motion in the initial flow; second, vertical motion of the pliolite with respect to the fluid. Errors are likely to be different between different runs of the model. To estimate error, two tracks very close in space but from different runs were compared (see section 4, Figure 7b.) The location difference after 2.5 s is 0.5 cm, with a resultant velocity error of 0.2 cm sÀ1. Other minor sources of error include the depth of the light sheets (estimated as plus or minus 0.1 cm top and bottom), and errors in identifying and digitizing the particles (estimated as plus or minus 0.1 cm). 2.2. Numerical Model Configuration [21] The S-Coordinate Rutgers University Model Version 3.1 [SCRUM, Song and Haidvogel, 1994] developed at Rutgers University was used for the numerical modeling part of this study. SCRUM is a hydrostatic primitive equation ocean circulation model with a free surface using a vertical s (terrain-following) coordinate well suited for use in simulations with variable bathymetry. The model also allows for a curvilinear horizontal coordinate system which enabled us to represent a cylindrical model geometry. [22] The combination of steep topography and strong stratification can cause errors in numerical simulations [Haney, 1991]. The errors, usually in the form of mean along isobath currents, are due to truncation errors in the calculation of the horizontal pressure gradient term in terrain-following coordinates. One current method to assess errors is to initialize the model with horizontal isopycnals and no forcing and evaluate the size of the observed error flows [e.g., Petruncio, 1996]. The resulting flows are due to pressure gradient errors and the boundary condition on density at the bottom which requires isopycnals to have no normal gradient at the bottom. The net effect of pressure gradient errors on strongly stratified flows over steep topography have been successfully assessed using laboratory experiments [Kliem and Pietrzak, 1999]. Steep topography and stratification can lead to other errors in s coordinate numerical models (as we show here). Such errors cannot be identified without forcing and have to date received little attention. [23] The large slopes of the tank canyon were too steep for the model to accurately compute horizontal pressure gradients. Therefore, for the stratification case only, the length scale was increased to ten times that of the tank. To ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL 3-5 temperature was used as an inverted proxy for density. For the stratified simulation, the initial vertical density profile was the same as for the tank (Figure 3) and the starting density was horizontally uniform. [26] At the start of the simulation the flow was zero and the free surface was level. The water was forced to move relative to the canyon by a body force over the entire water column. The forcing was zero at the start, increased to a maximum value over six seconds, held constant until 273 seconds and then set to zero. The forcing was in the angular direction only and increased linearly with radial distance in order to simulate the situation in the tank where the water is initially in rigid body rotation before the tank is accelerated with respect to the rotating water. Linear (Ekman) bottom stress was applied as a body force over the bottom layer (drag coefficient of 7.2 Â 10À5 m sÀ1). 3. Homogeneous Comparison Figure 4. The numerical grid for a cross-section through the canyon. The vertical scale is in cm. The horizontal scale is in grid points with Áx % 6.2 cm for this cross-section. The deepest cell at grid point 37 is used for the error calculation in section 5. use the same bathymetry, stratification and flow velocity as the tank, the Coriolis parameter was reduced (from 0.52 sÀ1 to 0.052 sÀ1) and the vertical viscosity was reduced (from the molecular value of 10À6 m2 sÀ1 to 10À7 m2 sÀ1) by a factor of 10 in order to maintain the same Rossby, Burger and Ekman numbers. The duration of the forcing was increased by a factor of 10 so that the model, like the tank, was accelerated over approximately one pendulum day. [24] The model domain included 242 points in the azimuthal direction, 81 points in the radial, and 16 vertical levels. The grid spacing in the vicinity of the canyon was about 6.2 cm (Figure 4). The domain did not cover the entire 2p radians in the azimuthal direction, but rather 0.6 Â 2p (in order to save computation time) with a periodic continuation in that direction. The duration of the simulations was much less than the time required for disturbances to wrap around the domain. In the radial direction, the model domain was based on a tank with an inner wall 1 m from the center, an outer wall 5 m from the center (10 Â the actual tank dimensions) and a no-slip condition for momentum at both lateral walls. The gridded model bathymetry was derived from the digitized tank depth contours using an objective analysis scheme. [25] The Laplacian horizontal mixing was along geopotential surfaces for both momentum and density with a coefficient of 10À6 m2 sÀ1. Vertical mixing of momentum and density used constant coefficients of 10À7 m2 sÀ1 and 1.5 Â 10À10 m2 sÀ1 (one tenth of the value for molecular diffusivity of salt in water [Weast, 1979]), respectively. A simplified linear equation of state dr = ÀdT means that [27] To compare the numerical and laboratory results, tracer particles were introduced into the numerical simulation at the time (0.063 rotation periods after acceleration was complete) and place tracer particles were observed to start in the laboratory experiment. Similarity and differences between the flows were determined by comparing the tracks of these particles (Figure 5). A value of 1 cm was chosen to delineate ‘‘good’’ and ‘‘poor’’ agreement. This value was determined using the estimate of the error in the laboratory of 0.2 cm sÀ1 and the maximum track durations of 5 s. For most tracks the final position of the particle according to the numerical model is within 1 cm of the final position according to the laboratory model (Figure 5a). The tracks show good agreement upstream of the canyon and fair agreement both across the mouth of the canyon and near the head of the canyon. [28] All tracks which finish more than 1 cm apart are downstream of the canyon inshore of the shelf break. In these cases the numerical tracks show stronger upwelling (stronger flow toward the coast) than the laboratory tracks (Figure 5b). The discrepancy suggests that greater cyclonic vorticity is generated within the canyon in the numerical model than in the laboratory model. All numerical tracks are bent strongly inshore over the canyon by this vorticity. [29] One explanation for the smaller vorticity generation in the laboratory comes from careful inspection of the vertical particle movement. The laboratory particles are color coded based on their depth. Upstream of the canyon, green particles become red particles as the flow drops down into the canyon (Figure 5b). On the downstream side, red particles become green particles as the flow is strongly upwelled. These numerous changes in color illustrate the strong vertical velocities occurring in the homogeneous case. We suspect that because the laboratory model is strongly vertically exaggerated nonhydrostatic flows occur in the laboratory for homogeneous flows. [30] To examine more carefully whether nonhydrostatic effects are indeed occurring in the laboratory model a visualization of the vertical shear downstream of the canyon was performed. In a homogeneous flow the pressure gradient is uniform with depth (barotropic) and vertical shear can only be generated by differential body forcing (not present) or in boundary layers. However, besides a clear and 3-6 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 5. A plan view of the particle tracks in the laboratory experiment (unfilled red and green symbols) and the numerical simulation (filled blue symbols) for the homogeneous case. The depth interval of the laboratory tracers was determined by the color-sheet illuminating them. The laboratory particles here are color-coded by depth, green: 1.2 – 2.2 cm depth, red: 2.2– 3.2 cm depth. The bathymetric contours are shown in blue with depth in cm. The axes are in cm in the tank frame. The last position is shown as a larger box. (a) Tracks that finish within 1 cm of each other. (b) Tracks that finish more than 1 cm apart. Note that these tracks are all downstream of the canyon. The weaker upwelling observed in the laboratory case is due to non-hydrostatic effects (see text). expected bottom Ekman layer, a middepth jet was observed with a strong up-canyon component below a surface layer flowing more slowly downstream (not shown). Thus middepth vertical shear was strong in the laboratory consistent with the presence of nonhydrostatic effects. [31] The numerical model is hydrostatic and thus constrained to have no vertical shear above the bottom Ekman layer. Fluid columns passing over the canyon are stretched to the full depth of the canyon, producing the observed strong cyclonic vorticity. In the laboratory model, columns passing over the canyon are only more moderately stretched. 4. Stratified Comparison [32] To compare the numerical and laboratory results, tracer particles were placed in the numerical simulation at the time and position tracer particles were observed to start in the laboratory experiment. In the stratified case, vertical shear is significant (Figure 6). In the laboratory, particle depths were determined, by color, only to within 1 cm depth. This resolution proved inadequate because of the strength of the vertical shear. Therefore in the numerical simulation, particles were seeded throughout the 1 cm depth interval and the track that finished closest to the laboratory track was chosen for comparison. [33] As discussed in the overview (section 1.1), flow in the canyon can be divided into three layers (Figure 1). Nearsurface flow is not affected or only weakly affected by the canyon. Flow deep within the canyon turns cyclonically. Active upwelling onto the shelf occurs between these two layers. The upwelling layer includes water from both the shelf and slope which turns cyclonically in the vicinity of the canyon, flows up the canyon and exits the canyon near the head. The discussion of the stratified results will be in reference to these layers. [34] The deep flow (below 2.15 cm) shows the expected trapped cyclonic eddy within the canyon (Figure 7d). The size of the eddy appears to be different between the numerical and laboratory realizations: the eddy is wider and longer in the lab, and shorter and narrower in the numerical simulations. Also in the region of the mouth of the canyon, the laboratory flow appears to slow more than the numerical flow. [35] Agreement between the numerical simulation and the laboratory experiment is generally very good in the actively upwelling layer (1.50– 2.05 cm depth, Figures 7b and 7c). The flow shows upwelling flow over the downstream rim and near the head of the canyon. Two systematic differences are apparent. Upstream of the canyon over the slope the laboratory tracks show weak downwelling flow (e.g., Figure 7b, track starting at x = À1, y = 2). This flow is not observed in the numerical simulation where the tracks follow the curvature of the tank. The difference in the tracks is about 0.5 cm and occurs over about 2 s giving a velocity error of 0.2 cm sÀ1. Downstream of the canyon the opposite situation occurs with weak downwelling flow in the numerical simulation but not in the laboratory flow (e.g., Figure 7b, track starting at x = À4, y = 2). The tracks of the numerical floats are consistent with the numerical velocity so these differences are not due to the Lagrangian tracking scheme. Flow above the canyon (between 1.25 cm to 1.35 cm depth, Figure 7a) shows the same systematic difference upstream of the canyon as seen in the active layer. ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL 3-7 Figure 6. A cross-section of the cross-canyon horizontal velocity at about mid-canyon for the stratified case. Positive velocities are to the right. View is toward the head of the canyon. Note the large vertical shear over the canyon. The vertical scale is 5 cm from top to bottom. Units are cm sÀ1. [36] To quantitatively measure the fit between the numerical results and the laboratory results, a normalized error was calculated for each particle at each point: Error ¼ j Positionðnumerical Þ À PositionðlaboratoryÞj LabFluctuationVelocity Ã time ð1Þ where the position is the (x, y) position in lab coordinates of the point, the laboratory fluctuation velocity is the velocity error between runs estimated above at 0.2 cm sÀ1 and time is the time since the particle started (Table 2). [37] The deep point near the mouth of the canyon that is within the deep eddy in the lab and not in the numerical simulation has an error at the end of 6.05. This point is very close to a flow separation point, dividing flow traveling into the canyon (as the laboratory tracer does) and flow traveling along the slope (as the numerical tracer does). Removing this point, the remaining points have an average endpoint error of 0.7 and an average error of 0.8. If the error estimate for the laboratory measurement was accurate and the model is reasonably accurate, a value near 1 would be expected. 5. Discussion [38] The utility of a numerical laboratory model comparison lies in resolving differences between them. In the following, we demonstrate the limitation of the numerical model simulating stratified shear flows over steep topography. The largest errors here are not due to the pressure terms as investigated by previous authors [e.g., Haney, 1991; Kliem and Pietrzak, 1999] because for our grid geometry these errors are small (see below) but are due to errors in vertical advection of the density field. [39] Although the average error in the stratified experiments is less than the estimated laboratory error, the results show systematic differences. These differences are most obvious near the upstream rim of the canyon where the laboratory flow is downwelling and the numerical flow is not. Downwelling flow in this region is observed in the field [Hickey, 1997; Allen et al., 2001]. A detailed search for the origin of the differences was made. Running the model with no forcing showed that pressure gradient errors are 50 times too small to explain the error. The timescale of the simulation is so short that errors due to diffusion are also too small to explain the observed differences. The details of the tracer advection scheme (fourth-order centered or upstreambiased third-order (Gamma scheme)) make little difference. Below we present an estimate of the nonhydrostatic effects in the laboratory and an estimate of the error in the vertical advection of density in the numerical simulation. Nonhydrostatic effects are shown to be too small and of the wrong sign to be the origin of the error. 5.1. Estimate of Nonhydrostatic Effects [40] The homogeneous case presented in section 3 illustrates that nonhydrostatic effects can sometimes explain the difference between laboratory and numerical models. In the stratified case however the strong stratification should be enough to prevent nonhydrostatic flows (conditions calculated by Boyer et al. [2000] are met). Here we will carefully evaluate the possible strength of the nonhydrostatic terms in the momentum equation for the stratified case. 3-8 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 7. Plan views of the particle tracks in the laboratory experiment (unfilled symbols) and the numerical simulation (filled symbols) for the stratified case. The bathymetric contours are shown with depth in cm. The axes are in cm in the tank frame. The ends of the tracks are shown as larger boxes. Note the generally good agreement. However, there are systematic errors upstream and downstream of the mouth of the canyon. In particular in the active layer and near surface layer at x = 1 cm, y = 1 – 2 cm the laboratory tracks show downwelling flow compared to the numerical tracks. (a) Tracks between 1.25 and 1.35 cm depth. (b) Tracks between 1.50 cm and 1.60 cm depth. (c) Tracks between 1.85 cm and 2.15 cm depth. (d) Tracks between 2.15 ( just above rim depth) and 3.15 cm depth. [41] The vertical momentum equation is @w @w @w @w ¼ ro þu þv þw @t @x @y @z 2 @p @ w @2w @2w À þ þ À r0 g þ m @z @x2 @y2 @z2 ð2Þ where ro is a constant reference density, (u, v, w) is the velocity, p is the pressure, r0 is the fluctuation density, g is the gravitational acceleration and m is the viscosity. [42] Scale analysis shows that the third largest term (after the pressure and gravitational terms) is the horizontal advection of the vertical momentum. After taking the horizontal derivative of equation (2) with respect to x and retaining only the three largest terms, equation (2) becomes À1 @ 2 p g @r0 @2w þu 2 ¼ ro @z@x ro @x @x ð3Þ ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Table 2. Initial Positions, Final Depth Position (Numerical) and Errors Between the Numerical and Laboratory Simulationsa Point Number X Start Y Start Z Start Z Finish Aver Error End Error 1 13 8 11 27 12 23 18 25 30 6 7 26 17 2 28 33 21 22 19 24 29 4 20 15 3 14 5 31 9 10 16 À3.8 4.6 À0.7 0.8 À0.7 2.9 4.4 1.8 0.0 1.6 4.5 À1.3 1.6 1.1 0.6 0.5 3.5 3.6 4.5 2.1 5.0 1.3 2.2 3.5 6.1 1.2 5.2 2.3 2.7 3.6 À3.9 À4.2 5.3 0.2 2.5 0.6 À4.4 0.4 À0.1 0.1 À2.0 À5.0 5.5 2.2 À1.1 À0.2 5.2 À3.6 0.8 À0.9 À0.2 0.4 À1.0 À5.1 4.2 À1.3 1.0 6.3 1.4 5.4 À6.2 3.9 1.1 À0.7 3.2 3.2 3.2 3.2 3.2 3.0 3.0 2.5 2.2 2.2 2.2 2.2 2.2 2.0 2.0 2.0 1.9 1.9 1.9 1.6 1.6 1.5 1.5 1.5 1.5 1.4 1.3 1.2 1.2 1.2 1.2 1.2 3.2 3.2 3.1 3.1 3.0 3.0 3.0 2.4 2.2 2.2 2.1 2.1 2.1 2.1 2.0 2.0 1.9 1.8 1.8 1.6 1.5 1.6 1.5 1.5 1.4 1.4 1.2 1.3 1.3 1.2 1.2 1.2 0.99 1.26 2.63 0.82 0.79 1.24 0.32 0.78 5.13 0.42 0.38 0.91 0.66 1.08 0.48 0.71 0.17 0.29 0.53 0.39 0.28 0.69 0.98 0.30 0.63 0.26 1.36 0.50 3.09 0.73 0.69 0.13 0.81 1.24 2.14 1.04 0.28 1.03 0.47 0.55 6.05 0.72 0.19 1.05 0.60 0.41 0.49 0.51 0.14 0.13 0.86 0.28 0.06 0.08 0.92 0.12 0.24 0.32 1.13 0.50 2.77 0.65 0.45 0.14 a Data are sorted in-depth order. 3-9 where the first term on the right is due to the hydrostatic pressure and the second term is the nonhydrostatic effect. [43] The largest terms in the cross-shore momentum equation are the pressure gradient and Coriolis terms (that is, the flow is approximately geostrophic). Substituting equation (3) into the z derivative of the cross-shore geostrophic balance gives @v Àg @r0 u @ 2 w À ¼ @z f ro @x f @x2 ð4Þ So the nonhydrostatic effect on the velocity shear is À u @2w f @x2 ð5Þ [44] To estimate the magnitude of the nonhydrostatic effects in the laboratory where w was not measured we need a magnitude for the vertical velocity. The numerical vertical velocity (Figure 8) is in error (see below). A local minimum in vertical velocity is expected over the upstream rim of the canyon but the maximum should not occur until the downstream rim of the canyon based on observed isopycnals over Astoria Canyon [Hickey, 1997]. Estimating the vertical velocity as the radial velocity (u = 1 cm sÀ1) multiplied by the slope of the canyon edge in the tank gives an upper estimate for the magnitude of the vertical velocity (w % À1 cm sÀ1). Using the width of the canyon as an appropriate scale for x implies equation (5) is of order À0.04 sÀ1. Because v is near zero near the surface, the nonhydrostatic effect gives positive (upwelling) flows, the Figure 8. A cross-section of numerical model true vertical velocity at about mid-canyon for the stratified case. Positive velocities are upward. The vertical scale is 5 cm from top to bottom. Units are cm sÀ1 in the numerical model. Note the small scale of the variations compared to the width of the canyon. 3 - 10 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 9. A cross-section of the density at about mid-canyon. The vertical scale is 5 cm from top to bottom. Units are st. 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 upwards to the right in a region observed to have strongly descending isopycnals. opposite of the observed effect. Assuming 1 cm is an appropriate depth for the strong downward flow gives an up-canyon, nonhydrostatic flow of 0.04 cm sÀ1. This error is one order of magnitude too small (the observed difference is 0.2 cm sÀ1). Since nonhydrostatic effects are the only reason the laboratory results could be ‘‘wrong’’, we now turn to the numerical model to find the origin of the labnumerical difference. 5.2. Vertical Velocity and Vertical Advection [45] The numerical vertical velocity (Figure 8) shows unexpected spatial small-scale features. In particular the upward vertical velocity, just inside the rim of the canyon, is counter to our physical understanding based on observations of isopycnals over Astoria Canyon [Hickey, 1997]. The downward velocities at the downstream rim are also questionable. We expect similar error near both rims but here will concentrate on the upstream rim. [46] Just below the upstream rim the numerical horizontal cross-canyon velocity is about 0.2 cm sÀ1 (Figure 6) the slope is about 1 in 20, the up-canyon velocity is about 0.05 cm sÀ1 and the slope is about 1 in 60. Thus from the bottom boundary condition (no flow through the canyon wall) we would expect a negative vertical velocity of about 0.01 cm sÀ1, not the positive velocities observed. [47] The vertical velocity is not a dynamic variable within the model; it is calculated as a diagnostic variable. Temporal fluctuations can thus be large but here are less than 1%. The vertically important dynamic within the model is the vertical advection of density (Figure 9). The calculated isopycnals reflect the calculated vertical velocity; they tilt gently down before the canyon and then abruptly upward over the canyon. The total deflection of the isopycnals in the numerical model results is very small (about one tenth the magnitude) compared to observations [Hickey, 1997]. If we assume the isopycnals are tilted in the laboratory as they are in the field we can estimate the generated down-canyon flow using the thermal wind equation (N2Di2/fDx) where the isopycnal drop Di observed in the field is about 25 m which would correspond to about 0.2 cm in the lab, and Dx is a lengthscale for the drop taken as one quarter of the canyon width. The resulting flow is 0.2 cm sÀ1 in agreement with the observed flow error. [48] The calculation of the advection of temperature in the numerical model uses the Shchepetkin and McWilliams [1998] Gamma advection scheme to reduce dispersive error. This scheme uses a large number of terms. The truncation error in the vertical velocity can be calculated for this scheme but the algebra is dense. To more easily illustrate the error, we will use the diagnostic calculation for the vertical velocity instead. The derivation for the Gamma advection scheme is in the Appendix. In the following, (1) the link between density advection and vertical velocity will be shown mathematically. Then (2) the truncation error in the vertical velocity will be calculated, followed by (3) the results (only) of the Gamma advection scheme. [49] The calculations will use the numerical values for a 10 m diameter tank. However, unlike the nonhydrostatic calculation above, the calculation of the advection error can simply be scaled to give values for a field simulation. Factors (based on Astoria Canyon) are 23 Â 103, 6.8 Â 103, 17, 5.0, for horizontal distance, vertical distance, 3 - 11 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL horizontal velocity and vertical velocity, respectively. Percentage errors are the same, regardless of scale. 5.2.1. Step One: Density Advection and Vertical Velocity [50] Consider a case in which the density (temperature here) only varies vertically. Then the time variation of temperature T is given by [52] Consider a single grid cell (Figure 10) near the upstream rim of the canyon (Figure 4). The velocity at the lower depth is taken as zero, and the velocity at the upper depth is U where @T @T ¼ Àw : @t @z ". The diagram shows that @z/@x = Àtan a [53] The vertical sigma velocity is w = Hz where = Ds/Dt and Hz = dz/ds. The value of w at the top of the cell is found from ð6Þ The advection of temperature is calculated within the model as @T @T @T @T ¼ Àu Àv À @t @x @y @s ð10Þ r Á ðuHz Þds ð11Þ ð7Þ ð8Þ where variation in y has been neglected and temperature has been assumed to only vary vertically. The problem with equation (8) is that the two terms on the right-hand side almost cancel, so that the calculation of @T/@t is based on the small difference between two large terms. Combining terms we can relate the density advection in the model to the vertical velocity in the following way: @T @z @z @T ¼À u þ @t @x @s @z @u @z Áx @z @x Z w¼À where s is the terrain-following vertical variable and = Ds/Dt. Equation (7) can be rewritten as @T @z @T @z @T ¼ Àu À @t @x @z @s @z U ¼À which, using the expansion formula [Hedstro¨m, 1997], and assuming v = 0 gives wjs ¼ wjðsÀÁsÞ À ðudzÞjðxþÁx=2Þ ÀðudzÞjðxÀÁx=2Þ Áx ¼U ÁZ ð12Þ Áx [54] The vertical velocity is calculated at the centre point as w¼ wjðsþÁs=2Þ þwjðsÀÁs=2Þ 2 þ À Á À @zÁ @z u @x þ u @x ð xþÁx=2Þ ð xÀÁx=2Þ 2 U ÁZ U @z ¼ þ 2Áx 2 @xð xÀÁx=2Þ ð13Þ ð9Þ where the term in brackets is the model vertical velocity as shown in equation (6). The numerical model does not use dT/dz (only @T/@x and @T/@s) and so equation (9) is not computed. However, this equation illustrates how the model calculated vertical velocity is related to the actual advection of a vertical temperature gradient. The vertical velocity is found from the small difference between the two large terms u@z/@x and À@z/@s. The actual calculation of the advection of temperature in the model is complicated by the chosen discretization. However, the error identified here is due to the estimation of the vertical velocity. The small differences due to the discretization will be discussed in step 3. 5.2.2. Step Two: Error in the Vertical Velocity [51] The second step is to estimate the truncation error in the model vertical velocity given by equation (9). The flow over the canyon is about 0.5 cm sÀ1 with a vertical shear of about 1.6 sÀ1. Only the shear contributes to the error in the vertical velocity, so we will set the background flow to zero. This choice reduces the magnitude of the terms in equation (9) so the fact the problem is the small difference between two large terms will be partially obscured. For reference, with the observed background flow u@z/@x % À0.02 cm sÀ1, from the model w = @z/@s % 0.01 cm sÀ1 and so the vertical velocity is about À0.01 cm sÀ1. We will assume that the flow at the bottom is parallel to the topography and that the flow in the up-canyon ( y) direction is zero. Variations in y will be neglected. These two terms are of opposite sign near the upstream canyon rim. The first, ÁZ/Áx is positive whereas the second, @z/@x is negative. Substituting for U gives the expression for the numerical w wnum ¼ 1 @u " ðÁZ À Áx tan a‘ Þ tan a 2 @z ð14Þ where a‘ = 1/2 (a1 + a2). [55] Now w at the bottom (here taken on the straight line between the two u grid points) is Àu tan ab where ab = 1/2 (a2 + a4) and as there is no horizontal change in u, w is constant with depth. So wreal ¼ 1 @u ðÁZ 0 À Áx tan ab Þ tan ab 2 @z ð15Þ The correct velocity is also found from the difference of two terms. There is truncation error in both numerical terms which becomes even more significant because the vertical velocity is the difference between them. [56] The four angles for the grid cell at the upstream corner are: a1 = 0.0094, a2 = 0.012, a3 = 0.086, a4 = 0.105. The size of this grid cell is Á x = 6.2 cm by ÁZ = 0.31 cm. The vertical shear was estimated from the numerical results (Figure 6) as @u/@z = 1.6 sÀ1 giving wnum = 0.0103 cm sÀ1. The length ÁZ0 = ÁZ + Áx/2 [tan(a2) + tan(a4) À tan(a1) À tan(a3)] is 0.373 cm and wreal = 0.00047 cm sÀ1. This difference results in an upward velocity error of 0.010 cm sÀ1 3 - 12 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure 10. Sketch of a single grid cell with a sharp change in orientation. This grid cell sits at the upstream rim of the canyon (Figure 4). Four angles (a1 – a4) and the vertical (ÁZ) and horizontal (Áx) grid sizes characterize its geometry. Horizontal flow is sheared with U at the higher left-hand grid-edge and 0 at the lower right-hand grid-edge. The vertical velocity is w and the cross-grid ‘‘vertical’’ velocity is w. in the same direction as the observed error. The error in the vertical velocity is of the same size as the actual velocity. This positive vertical velocity error at the upstream canyon edge gives a positive density anomaly and therefore positive (upwelling) flow error upstream of the canyon. [57] The grid cell at the upstream corner and the corresponding grid cell at the downstream corner are expected to give the largest errors in the domain because they have strong topography curvature (angles on the right are different from the angles on the left) and are in regions of strong vertical shear. [58] Equations (14) and (15) can be combined and assuming Áx and ÁZ are sufficiently small, the error can be shown to be of second order in these quantities. However, it is important to note that for the case considered here Áx and ÁZ are too large for the error to be dominated by the leading terms. 5.2.3. Step Three: Error in Advection [59] We have shown that (1) the vertical velocity and the vertical advection of temperature are linked and (2) that there are large truncation errors in the calculated vertical velocity calculation. However, the vertical velocity calculated above is not directly used by the model in the advection of temperature. One can, however, similarly calculate the equivalent vertical velocity used in the Gamma scheme for advection of a vertical gradient in temperature (see Appendix A). This gives À Á @T @T À ¼ 0:007 cm sÀ1 @t @z ð16Þ ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL 3 - 13 Figure 11. Grids and vertical velocities near the upstream rim for two grids with smaller grid spacing in the horizontal. (a) Grid with no smoothing. (b) Vertical velocity for grid with no smoothing. (c) Grid with smoothing of bottom topography. (d) Vertical velocity for grid with smoothing. Vertical velocity units are cm sÀ1 (positive upward) in the numerical model. Note the smaller vertical velocities in the smoother case. The error is smaller (by 30%) but still very significant at about 70% of the expected real vertical velocity (including the background flow). [60] To illustrate the effects of the truncation error on vertical velocity the model was run with two finer grids, one which does not reduce the size of the changes in bottom slope (Figure 11a) and one which does (Figure 11c). The first finer grid has a horizontal grid size one half that of the original grid (Figure 4) but the same bottom topography. The smaller horizontal grid size leads to smaller scales in the vertical velocity (Figure 11b) and no decrease in the magnitude of the positive velocities over the upstream rim. Close examination of equations (14) and (15) reveals that the difference is primarily dependent on the change of bottom angle. The first fine grid still contains a cell with the same angle change and thus the same large error in vertical velocity. [61] The second finer grid slightly smoothes the bottom topography (Figure 11c). The maximum change in angle is reduced by 25%. The vertical velocity (Figure 11d) shows a similar reduction in error. However, the vertical velocity error is still very large. The density field and particle trajectories are not significantly altered. For the current configuration, reducing the vertical grid size would 3 - 14 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Figure A1. Expanded view of grid in the vicinity of the upstream rim of the canyon. The center of the cell used for error estimation is marked by T. The geometry of the grid is needed farther upwards and upstream for the vertical advection calculation. Values of the angles and grid spacing are given. Note that the vertical scale is exaggerated. actually increase the error. Although the error is related to grid size and the numerical scheme does converge as the grid size approaches zero, reducing the error by reducing the grid size is onerous because both the vertical and horizontal grids must be reduced. Other possibilities for reducing the error are smoothing the topography or using an advection scheme that accounts for strong vertical shear. Calculation of the vertical advection errors for the smooth broad canyons used by Perenne et al. [2001] show that their truncation errors are small, consistent with the good match they found between laboratory and numerical results. However, canyons in the ocean have steep, sharp topography. 6. Conclusions [62] Upwelling over submarine canyons has been investigated using both a numerical and a physical model. Although qualitative results are similar, the results show systematic differences between the models. In the homogeneous case, the differences can be accounted for by nonhydrostatic effects in the laboratory model. [63] In the stratified case, nonhydrostatic effects are too small and have the wrong sign to account for the observed differences. However, the numerical truncation error in the calculation of the vertical advection was shown to be large enough and to have the correct sign to account for the differences. The presence of this error requires (1) strong vertical gradients in density, (2) vertical shear in the horizontal velocity, (3) topography with strong vertical curvature, and (4) the use of terrain-following coordinates. As the first three conditions are typically found in coastal environments, a solution to this error will be important for simulations that use realistic topography. In the model configuration used here, the error in the vertical advection is calculated to be 70%. As the canyon shape used is similar to Astoria Canyon the error for field-scale numerical simulations is the same, 70% in vertical advection. A similar error has been found by the authors in realistic simulations of Astoria Canyon. Because the error is dependent on both the vertical and horizontal grid size, reduction in grid size is an expensive approach to the problem. Higher order estimates of the vertical velocity that account for the vertical shear are likely required. ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Appendix A: Advection of a Vertical Gradient in Temperature where a0 and ÁZ1 are defined on Figure A1. Thus [64] In section 5.2, an error in the vertical velocity was calculated. Here we will calculate the dynamically important error, the error in the vertical density advection. Consider the same grid cell. To calculate vertical advection, we must consider flow further upstream and further up. Consider the extended mesh shown in Figure A1. [65] The change in temperature at the center point (marked T) should be À @T @T ¼ Àwreal @t @z @T 1 @u @T " ¼ tan a @t 4 @z @z Áx 2ÁZ ðtan a2 þ 2 tan a1 þ tan a0 Þ Á ðÁZ þ ÁZ1 À 2 ÁZ þ ÁZ 0 This expression is similar to wnumdT/dz with wnum given by equation (14). There are three differences. The expression for the vertical velocity wnum uses a2 rather than the weighted average of a0, a1 and a2 and uses ÁZ rather than the average of ÁZ and ÁZ1. The last factor does not appear in the calculation of wnum. The first change is a detriment, the second and third are improvements. [68] In particular, for the coarse grid: assuming that the temperature only varies vertically and wreal is given by equation (15). The spatial discretization for the calculation of the density advection is: @ 1 Áx Hz ð x þ Áx; sÞ þ Hz ð x; sÞ À ðHz T Þ ¼ u xþ ;s Á @t Áx 2 2 T ð x þ Áx; sÞ þ T ð x; sÞ 1 Áx Á À u xÀ ;s 2 Áx 2 Hz ð x; sÞ þ Hz ð x À Áx; sÞ T ð x; sÞ þ T ð x À Áx; sÞ Á Á 2 2 1 Ás Ás Á Hz x; s þ þ x; s þ Ás 2 2 T ð x; s þ ÁsÞ þ T ð x; sÞ 1 Ás À x; s À 2 Ás 2 Ás T ð x; sÞ þ T ð x; s À ÁsÞ Á Hz x; s À ðA1Þ 2 2 neglecting the additional ‘‘Gamma’’ term in the Shchepetkin-McWilliams advection scheme. It will be discussed below. To calculate the value for the center of the grid cell we will ignore the uniform horizontal velocity and include only the vertical shear. Thus without loss in generality we assume u(x + Áx/2, s) = 0 and u(x À Áx/2, s) = U. As w(x, s À Ás/2) is on the topography it is zero and as before Hz = w. We note that 0.5 (Hz(x À Áx, s) + Hz(x, s)) = ÁZ/Ás, Hz(x, s) = 0.5 (Z + Z0)/Ás, and that equation (12) gives w(x, s + Ás/2) = UÁZ/Áx. [66] Then equation (A1) becomes ÁZ þ ÁZ 0 @T U ÁZ À ¼À ðT ð x; sÞþT ð x À Áx; sÞÞ 2Ás @t 2Áx Ás U ÁZ ðT ð x; s þ ÁsÞþT ð x; sÞÞ þ 2Ás Áx ðA2Þ [67] Using the geometry shown in Figure A1, we can expand @u U¼ tan a " Áx @z 1 dT Áx ÁZþÁZ1 þ xðtan a2 À tan a0 Þ T ð x; s þ ÁsÞ ¼ T ð x; sÞ þ 2 dz 2 dT Áx T ð x À Áx; sÞ ¼ T ð x; sÞ þ ðtan a1 þ tan a2 Þ dz 2 3 - 15 À Á @T @T À ¼ 0:009cm sÀ1 ; @t @z which is an improvement of 10% reducing the error to 90% of the expected total vertical velocity including that background horizontal velocity as well as the shear. [69] The Shchepetkin-McWilliams advection scheme has an additional term (Àd1) in the compution of a tracer flux through a grid cell surface where À = À1/4 and d1 = T(x) À 2T(x À Áx) + T(x À 2Áx) for U positive. This term only affects the horizontal advection. Also, for the vertical temperature advection a four-grid point calculation is used. Using these expressions instead of the simple second order center difference scheme, equation (A2) for our chosen box becomes ÁZ þ ÁZ 0 @T U 9 1 À ¼ T ð x; s þ ÁsÞ À T ð x; s þ 2ÁsÞ 2ÁZÁs @t Áx 16 16 ! 3 1 1 À T ð x À Áx; sÞ þ T ð x; sÞþ T ð x À 2Áx; sÞÞ 4 8 8 where we note that the code uses the no-heat flux bottom boundary condition. To evaluate this expression it is necessary to consider the structure of the grid further upwards and farther to the right. The geometric parameters are given on Figure A1 and the result is À Á @T @T À ¼ 0:007cm sÀ1 @t @z [70] This scheme is an improvement over the direct calculation of w and reduces the error to about 70% of the estimated vertical velocity. Thus the vertical advection of the velocity field (which should be strongly downward in this region) is reduced to 70%. [71] Acknowledgments. This work was supported by a National Science and Engineering Research Council of Canada Strategic Grant to S. Allen and by National Science Foundation Grants OCE 96-18239 to J. Klinck and OCE 96-18186 to B. Hickey. Kurt Nielson was instrumental in the design and building of the laboratory apparatus. The authors thank an anonymous reviewer for his/her detailed and helpful comments on a previous version of the manuscript. References Allen, S. E., Topographically generated, subinertial flows within a finite length canyon, J. Phys. Oceanogr., 26, 1608 – 1632, 1996. 3 - 16 ALLEN ET AL.: COMPARISON OF A NUMERICAL AND A PHYSICAL MODEL Allen, S. E., M. S. Dinniman, J. M. Klinck, and B. M. Hickey, Dynamics of advection-driven upwelling over a submarine canyon, Eos Trans. AGU, 80, 221, 1999. Allen, S. E., C. Vindeirinho, R. E. Thomson, M. G. G. Foreman, and D. L. Mackas, Physical and Biological processes over a submarine canyon during an upwelling event, Can. J. Fish. Aquat. Sci., 58, 671 – 684, 2001. Batchelor, G. K., An Introduction to Fluid Dynamics, 615 pp., Cambridge Univ. Press, New York, 1967. Boyer, D. L., X. Zhang, and N. Pe´renne, Laboratory observations of rotating, stratified flow in the vicinity of a submarine canyon, Dyn. Atmos. Oceans, 31, 47 – 72, 2000. Butman, B., Long term current measurements in Lydonia and Oceanographer Canyons, Eos Trans. AGU, 64, 1050 – 1051, 1983. Gent, P. R., J. Willibrand, T. J. McDougal, and J. C. McWilliams, Parameterizing eddy-induced tracer transports in ocean circulation models, J. Phys. Oceanogr., 25, 463 – 474, 1995. Haidvogel, D. B., and A. Beckmann, Numerical models of the coastal ocean, in The Sea, edited by K. H. Brink, and A. R. Robinson, vol. 10, chap. 17, pp. 457 – 482, John Wiley, New York, 1998. Haney, R. L., On the pressure gradient force over steep topography in sigma coordinate ocean models, J. Phys. Oceanogr., 21, 610 – 619, 1991. Hedstro¨m, K. S., User’s manual for an S-coordinate primitive equation ocean circulation model (SCRUM) version 3.0, Contrib. 97-10, Inst. of Mar. and Coastal Sci., Rutgers Univ., New Brunswick, N. J., 1997. Hickey, B. M., The response of a steep-sided narrow canyon to strong wind forcing, J. Phys. Oceanogr., 27, 697 – 726, 1997. Kinsella, E. D., A. E. Hay, and W. W. Denner, Wind and topographic effects on the Labrador current at Carson canyon, J. Geophys. Res., 92, 10,853 – 10,869, 1987. Kliem, N., and J. D. Pietrzak, On the pressure gradient error in sigma co-ordinate ocean models: Comparison with a laboratory experiment, J. Geophys. Res., 104, 29,781 – 29,799, 1999. Klinck, J. M., Circulation near submarine canyons: A modeling study, J. Geophys. Res., 101, 1211 – 1223, 1996. Osler, G., Density gradients, Sci. Am., 213, 70 – 76, 1965. Perenne, N., D. B. Haidvogel, and D. L. Boyer, Laboratory-numerical model comparisons of flow over a coastal canyon, J. Atmos. Oceanic Technol., 18, 235 – 255, 2001. Petruncio, E. T., Observations and modeling of the internal tide in a submarine canyon, Ph.D. thesis, 176 pp., Nav. Postgrad. School, Monterey, Calif., 1996. Shchepetkin, A. F., and J. C. McWilliams, Quasi-Monotone advection schemes based on explicit locally adaptive dissipation, Mon. Weather Rev., 126, 1541 – 1580, 1998. She, J., and J. M. Klinck, Flow near submarine canyons driven by constant winds, J. Geophys. Res., 105, 28,671 – 28,694, 2000. Song, Y. H., and D. B. Haidvogel, A semi-implicit ocean circulation model using a generalized topography-following coordinate system, J. Comput. Phys., 115, 228 – 244, 1994. Weast, R. C., editor, Handbook of Physics and Chemistry, Chem. Rubber. Publ. Co., 59th edition, 1979. ÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀ S. E. Allen, D. D. Gorby, and A. J. Hewett, Oceanography, Department of Earth and Ocean Sciences, University of British Columbia, 6270 University Boulevard, Vancouver, B. C., Canada, V6T 1Z4. (allen@ocgy.ubc.ca) M. S. Dinniman and J. M. Klinck, Crittenton Hall, Center for Coastal Physical Oceanography, Old Dominion University, Norfolk, VA 23529, USA. B. M. Hickey, School of Oceanography, University of Washington, Box 357940, Seattle, WA 98195, USA.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- On vertical advection truncation errors in terrain-following...
Open Collections
UBC Faculty Research and Publications
On vertical advection truncation errors in terrain-following numerical models: Comparison to a laboratory… Allen, Susan E.; Dinniman, M. S.; Klinck, J. M.; Gorby, D. D.; Hewett, A. J.; Hickey, B. M. Jan 31, 2003
pdf
Page Metadata
Item Metadata
Title | On vertical advection truncation errors in terrain-following numerical models: Comparison to a laboratory model for upwelling over submarine canyons |
Creator |
Allen, Susan E. Dinniman, M. S. Klinck, J. M. Gorby, D. D. Hewett, A. J. Hickey, B. M. |
Publisher | American Geophysical Union |
Date Issued | 2003-01 |
Description | Submarine canyons which indent the continental shelf are frequently regions of steep (up to 45°), three-dimensional topography. Recent observations have delineated the flow over several submarine canyons during 2–4 day long upwelling episodes. Thus upwelling episodes over submarine canyons provide an excellent flow regime for evaluating numerical and physical models. Here we compare a physical and numerical model simulation of an upwelling event over a simplified submarine canyon. The numerical model being evaluated is a version of the S-Coordinate Rutgers University Model (SCRUM). Careful matching between the models is necessary for a stringent comparison. Results show a poor comparison for the homogeneous case due to nonhydrostatic effects in the laboratory model. Results for the stratified case are better but show a systematic difference between the numerical results and laboratory results. This difference is shown not to be due to nonhydrostatic effects. Rather, the difference is due to truncation errors in the calculation of the vertical advection of density in the numerical model. The calculation is inaccurate due to the terrain-following coordinates combined with a strong vertical gradient in density, vertical shear in the horizontal velocity and topography with strong curvature. An edited version of this paper was published by AGU. Copyright 2003 American Geophysical Union. |
Genre |
Article |
Type |
Text |
Language | eng |
Date Available | 2011-05-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0041900 |
URI | http://hdl.handle.net/2429/34544 |
Affiliation |
Science, Faculty of Earth and Ocean Sciences, Department of |
Citation | Allen, Susan E., Dinniman, M.S., Klinck, J.M., Gorby, D.D., Hewett, A.J., Hickey, B.M. 2003. On vertical advection truncation errors in terrain-following numerical models: Comparison to a laboratory model for upwelling over submarine canyons Journal of Geophysical Research Oceans 108(C1) 3003 |
Publisher DOI | 10.1029/2001JC000978 |
Peer Review Status | Reviewed |
Scholarly Level | Faculty |
Copyright Holder | Allen, Susan E. |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 52383-Allen_AGU_2001JC000978.pdf [ 1.13MB ]
- Metadata
- JSON: 52383-1.0041900.json
- JSON-LD: 52383-1.0041900-ld.json
- RDF/XML (Pretty): 52383-1.0041900-rdf.xml
- RDF/JSON: 52383-1.0041900-rdf.json
- Turtle: 52383-1.0041900-turtle.txt
- N-Triples: 52383-1.0041900-rdf-ntriples.txt
- Original Record: 52383-1.0041900-source.json
- Full Text
- 52383-1.0041900-fulltext.txt
- Citation
- 52383-1.0041900.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.52383.1-0041900/manifest