Flow dynamics around downwelling submarine canyons

Flow dynamics around a downwelling submarine canyon were analysed with the Massachusetts Institute of Technology general circulation model. Blanes Canyon (northwestern Mediterranean) was used for topographic and initial forcing conditions. Fourteen scenarios were modelled with varying forcing conditions. Rossby and Burger numbers were used to determine the significance of Coriolis acceleration and stratification (respectively) and their impacts on flow dynamics. A new non-dimensional parameter ( χ ) was introduced to determine the significance of vertical variations in stratification. Some simulations do see brief periods of upwards displacement of water during the 10-day model period; however, the presence of the submarine canyon is found to enhance downwards advection of density in all model scenarios. High Burger numbers lead to negative vorticity and a trapped anticyclonic eddy within the canyon, as well as an increased density anomaly. Low Burger numbers lead to positive vorticity, cyclonic circulation, and weaker density anomalies. Vertical variations in stratification affect zonal jet placement. Under the same forcing conditions, the zonal jet is pushed offshore in more uniformly stratified domains. The offshore jet location generates upwards density advection away from the canyon, while onshore jets generate downwards density advection everywhere within the model domain. Increasing Rossby values across the canyon axis, as well as decreasing Burger values, increase negative vertical flux at shelf break depth (150 m). Increasing Rossby numbers lead to stronger downwards advection of a passive tracer (nitrate), as well as stronger vorticity within the canyon. Results from previous studies are explained within this new dynamic framework.


Introduction
Submarine canyons are features typical of continental slopes and deeply incise the continental shelf.On a regional scale, physical processes (such as upwelling and cross-shelf exchange) can be modified/enhanced due to the presence of a submarine canyon (Hickey, 1995).In particular, numerical models have shown that submarine canyons can enhance upwelling/downwelling in coastal regions (Klinck, 1996).
High biological productivity is typically associated with upwelling canyons (Bosley et al., 2004;Allen et al., 2001); however, downwelling canyons can also be very productive (Mann, 2002;Skliris and Djenidi, 2006;Flexas et al., 2008).Some of these canyons -for example, the Gully of Nova Scotia, Canada -are probably productive due to strong mixing (Le Souëf and Allen, 2014).However, tides are small in the Mediterranean Sea and yet submarine canyons along the Catalan continental margin (northwestern Mediterranean Sea) support important commercial fisheries (Company et al., 2012).This may be associated with upwelling around these predominately downwelling canyons (Flexas et al., 2008).
The direction of alongshore flow is critical to the circulation over a canyon (Klinck, 1996).For right-bounded coastal flows, the geostrophic pressure gradient is offshore.Away from the canyon, this pressure gradient force is balanced by the Coriolis force due to the alongshore flow.However, within the canyon, the alongshore flow is reduced due to canyon walls.This creates an unbalanced offshore pressure gradient, which tends to drive flow offshore, leading to downwelling (Freeland and Denman, 1982).For leftbounded flows, the unbalanced pressure gradient force is onshore, leading to upwelling (ibid).
The right-bounded current, approaching the canyon along the shelf, slows and descends as it crosses the canyon (Klinck, 1996).Once it crosses the axis of the canyon, it accelerates and begins to rise (ibid).With weak dissipation, water returns almost to its original depth and continues alongshore (Klinck, 1996).The downwelling response of right-bounded flows is generally smaller than the upwelling response of left-bounded flows (ibid), so that oscillatory flow usually leads to weak upwelling from slightly below shelf break depth (Boyer et al., 2004).
Stratification controls the magnitude of a forcing response and limits the influence of a canyon on overlying flow, independent of the direction of alongshore flow (Klinck, 1996).Increased stratification reduces vertical and crossshore transport, as well as depth range over which fluid parcels move in a circuit around a canyon (Klinck, 1996;Skliris et al., 2001Skliris et al., , 2002)).Variations in flow strength impact horizontal and vertical flux exchange.Increasing the Rossby number due to a wind event (occurring in the same direction as alongshore flow) drives stronger net cross-shore and net vertical transports (Skliris et al., 2001(Skliris et al., , 2002)).
Observational cruises near Palamòs Canyon (northeastern edge of Spain) reveal that small-scale variability in the onshore/offshore location of an incoming zonal jet has important impacts on flow dynamics (Alvarez et al., 1996).Transient factors (such as river runoff and the climatology of the area) induce a series of modifications in the permanent frontcurrent of the region, affecting both its vertical extension and offshore location (ibid).When incoming zonal jets are near the head of the canyon, flow is narrower and faster, and vertical velocities are greater (relative to an incoming zonal jet placed further offshore) (Jordi et al., 2005).In addition, in areas where canyon width is narrow and depth variations are strong (i.e.canyon head), vorticity adjustments and associated vertical velocities are induced as part of ageostrophic adjustment, and the core of the front-current is displaced offshore downstream of the canyon (Alvarez et al., 1996).However, in areas where depth changes are not as strong and the canyon is wide (i.e.canyon mouth), flow adjustment is almost geostrophic and vertical velocity allows flow to maintain a thermal wind balance (Alvarez et al., 1996).
Although upwelling canyons have been studied more thoroughly than downwelling canyons, a number of studies have been done on downwelling canyons.Careful evaluation of the results of these studies shows at least two patterns of flow over canyons.Some studies show flow nearly following the isobaths around the canyon, leading to positive vorticity at shelf break depth (Klinck, 1996;Skliris et al., 2001Skliris et al., , 2002) ) (Table 1), whereas other studies show a trapped anticyclone or negative vorticity in the canyon (She and Klinck, 2000;Flexas et al., 2008) (Table 1).
In simulations with flow nearly following isobaths, antisymmetrical upwards (downwards) vertical velocity is seen in the downstream (upstream) region of a canyon, with upwards velocity being less intense than downwards velocity (Klinck, 1996;Skliris et al., 2001Skliris et al., , 2002)).In these cases, negative density anomalies occur everywhere over the canyon (with positive anomalies both upstream and downstream of the canyon in the slope region).
In simulations with a trapped anticyclone or negative vorticity, differences occur between studies.Under constant downwelling winds, a strong anticyclone within a canyon (in the upper 200 m) and small net cross-shore exchanges are driven by vortex compression or frictional coupling to alongshore flow (She and Klinck, 2000).Vertical flux is downwards everywhere over the canyon at shelf break depth (ibid).Observations over Blanes Canyon (BC) reveal that flow near the shelf break follows isobaths along canyon walls, with weak circulation in the canyon head (Flexas et al., 2008).In the upper 100 m, circulation is cyclonic along the canyon mouth, but weakly anticyclonic within the canyon.Vertical velocity estimates reveal net downwards transport in the upper 100 m, but net upwards transport between 100 and 200 m (Flexas et al., 2008).Density sections suggest local downwelling/upwelling occurs along the upstream/downstream canyon walls between 100 and 200 m depth (ibid).
In the studies of downwelling submarine canyons, there does not appear to be a clear agreement on flow dynamics (Table 1).This study attempts to better understand these and other differences between previous studies, as well as resolve the parameters that drive general flow dynamics in downwelling submarine canyons.The specific objectives of this study are (1) to determine whether upwelling occurs in or around downwelling canyons and (2) to determine what parameters affect flow dynamics -in particular, which parameters impact horizontal circulation, vertical transport, density advection, and passive tracer advection.
For the purpose of this study, canyon topography is based on the bathymetry of BC (Fig. 1).BC lies along the Catalan coast, in the northwestern Mediterranean Sea and is one of the few submarine canyons for which there have been multiple observational and numerical studies.The core simulation is a model of BC which replicates observations of Flexas et al. (2008).However, simulations of Klinck (1996), She and Klinck (2000), and Skliris et al. (2001Skliris et al. ( , 2002) ) are also replicated, and nine other simulations with variations in flow, stratification, topography, and boundary conditions were modelled.Section 2 describes the numerical model, domain bathymetry, parameter choices, modelled cases, and analysis calculations.Section 3 details results for all model simulations.Section 4 discusses the impacts of various parameters to modelled dynamics.Section 5 examines the significance of the results.Canyon mouth is the open region along the shelf break, canyon head is the shallowest onshore canyon region, mid-canyon is the region between the canyon head and mouth, and lower canyon is offshore of the canyon mouth.Canyon wall refers to canyon topography between shelf break and bottom depth.Zonal flow is in the alongshore direction, and meridional flow is in the cross-shore direction.Contour intervals are 100 m.

Model description
Simulations were run with the Massachusetts Institute of Technology general circulation model (MITgcm) (Adcroft et al., 2004).The model is rooted in incompressible Navier-Stokes equations; non-hydrostatic terms were used for all simulations.

Domain and canyon bathymetry
The model domain is 120 km in the alongshore direction (x direction), 90 km in the cross-shore direction (y direction), and 1200 m in the vertical direction (z direction) (Fig. 1).Positive x points upstream (eastward), positive y points onshore (northward), and positive z points upwards.Minimum ocean depth is 20 m and stretches for ∼ 20 km in the cross-shore direction (hereafter referred to as the inner shelf).Between the inner shelf and shelf break lies the outer shelf; in this region, depth drops to 150 m over 20 km.The slope extends from the shelf break (150 m) to an abyssal depth of 1200 m, and extends over 25 km in the cross-shore direction.The canyon topography was based on BC bathymetry shown by Flexas et al. (2008).Geometric parameters were kept nearly constant in all model simulations (Table 2).

Parameter specifications
Temperature, salinity, and nitrate stratification in the model were based on data from the National Virtual Ocean Data System (NVODS, http://ferret.pmel.noaa.gov/NVODS/UI.vm) (Fig. 2 temperature and salinity were collected from the World Ocean Atlas 2005 1 • × 1 • monthly means at approximately 40.5 • N, 2.5 • E. Values for nitrate were collected from annual means of the same data set at the same position.For runs with high vertical resolution, values between data points were linearly interpolated.A linear equation of state was applied, with a thermal expansion coefficient of 2.0×10 −4 ( • C −1 ) and a haline contraction coefficient of 7.4 × 10 −4 .
Horizontal resolution varies in alongshore and cross-shore directions; grid spacing is ∼ 1 km along each boundary and decreases linearly to 200 m over the canyon.Overall, there is 200 m of horizontal spacing between 33 and 87 km in the alongshore direction, and 20 km to 80 km in the cross-shore direction.Ninety vertical layers are concentrated around the top of the domain, and vertical spacing ranges from 5 m (in the upper 200 m) to 20 m (everywhere below 200 m).
The Coriolis parameter was assumed constant (f = 1.0 × 10 −4 s −1 ).Bottom friction was parameterised with a quadratic drag coefficient of 2.0 × 10 −3 .A vertical eddy viscosity of 1.0 × 10 −2 m 2 s −1 was applied.The model used non-hydrostatic equation sets, with a time step of 40 s for all runs.Viscous (i.e.no-slip) conditions were applied at the sides and bottom of the domain, and an implicit free surface was used.Heat and salt were laterally and vertically diffused with a Laplacian diffusivity of 1 × 10 −7 m 2 s −1 .A Smagorinsky harmonic viscosity factor (Smagorinsky, 1963) of 2.2 was applied (as recommended in Griffies and Hallberg, 2000).All tracers (i.e.temperature, salinity, and nitrate) were advected in time using a third-order direct spacetime scheme with flux limiting.
All model scenarios had a closed (no-slip) boundary along the coastal boundary.The offshore boundary was open with an Orlanski (1976) radiation condition applied.All but two simulations used periodic alongshore conditions; these two simulations will be explained further in the next section.

Model simulations
All modelled scenarios were forced by applying a wind stress and/or body force over the domain.A body force was applied as an additional forcing to the momentum equations (Dawe and Allen, 2010).Fourteen scenarios were modelled based on minor changes in either domain stratification or forcing (Tables 3 and 4).Two non-dimensional parameters were calculated to highlight incoming velocity (Rossby number, R o ), and stratification (Burger number, B u ), Dynamic parameters are incoming velocity, U ; the Coriolis parameter, f ; and stratification characterised by the buoyancy frequency at shelf break depth (150 m), N sb .Geometric parameters are length of the canyon, L, and depth at the shelf break, H sb .
To better understand the impact changes in stratification have on flow dynamics, a third parameter -a nondimensional measure of vertical stratification, χ -was introduced.This new parameter measured uniformity of stratification and was calculated as the change in buoyancy frequency (N) divided by the average buoyancy frequency near shelf break depth: where N (z) is measured over a length scale of ±75 m from the shelf break and N (z) is the average stratification over the length scale.Negative χ values indicate stronger stratification in the shallower layers.
In addition, placement of incoming coastal jets varied, both vertically and horizontally, in previous studies (Table 4).These were recreated to ensure that the dynamics of the original studies were reproduced.
Core model simulations were based on the Flexas et al. (2008) observations.The first scenario (uniform wind, UW) consisted of a uniform wind stress (τ = −0.0626N m −2 ) to drive a current along the surface and a body forcing (applied near shelf break depth) to drive a current along the shelf break (similar to the Northern Current seen in the Mediterranean Sea).The current was accelerated over the first two model days, and then held at a "steady" state for the remainder of the simulation ("steady" state indicates maximum flow velocity never varied more than 20 %).
The second scenario (opposing wind, OW) consisted of two opposing wind stresses to drive surface flow and a slightly stronger body forcing to drive the shelf break current.This setup was used to reproduce eastward flow seen over the continental shelf (Flexas et al., 2008).To match the offshore distance of the eastward flow, wind stresses were applied such that the offshore two-thirds of the domain had a wind stress of τ = −0.0626N m −2 and the nearshore onethird of the domain had a wind stress of τ = +0.0376N m −2 .Again, the current was increased during the first two model days.
Additional scenarios were modelled to either recreate flow dynamics seen in previous numerical studies (three scenarios) or investigate other impacts to flow dynamics (nine scenarios).
A Klinck-like (KL) scenario was modelled using a uniform stratification (N = 0.0016 s −1 ) and a flat shelf at 150 m (topography everywhere else in the domain remained the same) (Klinck, 1996).A mostly uniform flow was reproduced by removing all wind stress and y dependence in the body forcing.However, flow over the flat shelf was weaker relative to flow along the continental shelf and over the open ocean.Speed of the body forcing was reduced to create a zonal velocity of about 10 cm s −1 .
To simulate the She and Klinck (2000) study, a constant weak body forcing was applied over the upper 40 m, generating a maximum zonal speed of ∼ 13 cm s −1 (She).Stratification was varied over the entire depth, based on the equation for density provided in the original study.Initial fields are temperature (T ), and salinity (S), where z is depth below 0.
To better understand the impact of open vs. periodic boundary conditions, two scenarios with open alongshore boundaries were modelled.In these simulations, Orlanski radiation conditions were applied across both alongshore boundaries and the offshore boundary.The first scenario (open boundary conditions, OBC) has the same geometry as the UW case, but both wind stress and body forcing were increased to recreate a similar zonal flow field as seen in the UW simulation.The second simulation (slanted topography, ST) used geometry that was similar to real-world BC bathymetry (i.e. a slanted coastline and curvature within the canyon).Forcing in this scenario was the same as the OBC case.
Two scenarios of constant surface-forced (SF) flow were modelled, one forced by a wind stress and one forced with a surface body force applied to the upper 30 m. Results from these simulations were very similar, and will therefore be discussed as one example.
A uniform stratification (US) scenario was modelled with the same geometry and forcing as the UW scenario, but stratification from the KL case was used (N = 0.0016 s −1 everywhere).Similarly, a high Burger number (HB) scenario was modelled.This case is almost exactly the same as the uniform stratification scenario, but with a uniform N value of 0.005 s −1 .
Four final scenarios were modelled using various parameter specifications from previous simulations.To generate a simulation with low Rossby and χ values (LRC), SK forcing was applied, but with the same stratification as the HB case.The SK scenario was modelled again, but with a stronger forcing to generate a similar simulation but with a high Rossby number (SHR).Core case forcing and stratification was reduced to generate a simulation with low Rossby and Burger values (BLRB).The barotropic forcing (used in KL) was increased and run with core case stratification to produce high Rossby and Burger values (KHRB).
For all simulations with a wind forcing, the wind stress was linearly ramped over the initial model day, then held constant for model days 1-10 (wind magnitude; Table 5).All but one simulation with a body forcing (She) was linearly increased during the first model day.For these scenarios, a constant force was applied over model days 1-2 (body force magnitude, Table 5), followed by a linear decrease at the same rate as the increase, down to a constant value which was maintained to the end of the simulation (Table 5).For the She simulation, a body force was linearly ramped over the initial model day, then a constant body force was applied for model days 1-10.The two simulations with constant forcing (She and SF) are not steady in time and experience a large time dependence in their flow dynamics.However, none of the conclusions/trends discussed in this study are dependent on the results from these two simulations.

Result calculations
Transport calculations were used to estimate the volume of water exchanged vertically and horizontally in the domain.An initial plane along the canyon axis divides the canyon into an upstream and downstream half (Fig. 3).Zonal flux was calculated across this plane from surface to shelf break depth, and from the canyon mouth to coastal boundary (U 3).Meridional flux was calculated across two regions that lie along the canyon mouth: one in the upstream region (V 2) and one in the downstream (V 3).Again, flux was calculated from surface to shelf break depth.Finally, two planes were used to calculate vertical flux at shelf break depth.These planes extend from the canyon head to canyon mouth and split across the canyon axis (upstream: W 1; downstream: W 2). Net vertical flux was calculated by summing flux across these two planes.Flux across all planes was found by multiplying velocity of each grid cell by area of each grid cell and summing over the entire plane.
Relative vorticity in the basin can be expressed as where ζ is the vertical component of vorticity, V is the meridional velocity, and U is zonal velocity.Absolute vorticity was measured as the relative vorticity divided by the Coriolis parameter, f .Average zonal velocity across the canyon axis at shelf break depth (150 m) was used to calculate a second Rossby number (R U can ).This velocity is calculated as where U is taken as the zonal velocity in each meridional grid point that lies along the canyon axis, and y is the horizontal distance the zonal velocity is applied.A canyon Rossby number was calculated as Density was calculated for all model simulations as average density in the canyon across the shelf break plane (W 1 and W 2).This value was averaged during the approximate advection-dominated phase (averaged from model days 4 to 10).Average density was subtracted from initial density at shelf break depth to give an average density anomaly in the canyon.
An average density anomaly was also calculated in nine regions across the canyon domain.These values were averaged over the approximate advection-dominated phase (averaged over model days 4-10).The nine regions are as follows: -US_can, DS_can: width from canyon axis to upstream rim or downstream rim, respectively; length from the shelf break to the coast; and depth from shelf break depth to bottom depth.
-US_shallow, DS_shallow: same as above, but depth from surface to shelf break depth.
-US_shelf, DS_shelf: width from upstream rim or downstream rim, respectively, to 10 km away from canyon; length from shelf break to coast; and depth from surface to shelf depth.
-Lower_can: width from upstream rim to downstream rim, length from shelf break to 15 km offshore, and depth from surface to shelf break depth.
-US_off, DS_off: width from upstream rim or downstream rim, respectively, to 10 km away from canyon; length from shelf break to 15 km offshore; depth from surface to shelf break depth.
Changes in density difference within the canyon relative to away from the canyon were determined by calculating a density difference anomaly.This anomaly was found by subtracting a background density difference (calculated as a five-grid-point average along the downstream boundary) from the difference at grid points of similar isobaths: where x i and y i are alongshore and cross-shore points (respectively) that are ±5 m of the isobath used to calculate the background density difference (ρ boundary (y i , z)).
Nitrate concentration was used as a passive tracer in the model.An average nitrate concentration was also calculated as the average nitrate value in the canyon across the shelf break plane (W 1 and W 2) during the advection dominant phase (model days 4-10).The average nitrate concentration is subtracted from the initial nitrate concentration at shelf break depth to give an average nitrate anomaly in the canyon.

Flow evolution
All model scenarios show an initial time-dependent response to model forcing, similar to that described by Allen and Durrieu de Madron ( 2009), which lasts approximately 2 to 3 days (Fig. 4).During this phase, zonal and vertical flux exhibit a negative ramping everywhere in the domain.In all but two simulations, vertical flux across the downstream plane (W 2; Fig. 3) reverses at approximately day 1, and continues towards a maximum positive value by day 2-2.5.In these simulations, the magnitude of meridional flux over the canyon gently increases across both planes, being positive (onshore) in the upstream region (V 2; Fig. 3) and negative (offshore) in the downstream (V 3; Fig. 3) until reaching a maximum near the end of the time-dependent phase.In the simulations with a coastal jet (She and SF) across both the upstream and downstream planes until day 1 (Fig. 4b).After this, negative flux across the downstream plane weakens in time, while negative flux across the upstream plane continues to strengthen.For these scenarios, upstream onshore flux and downstream offshore flux strengthen during the model simulation.
The time-dependent phase is followed by an advectiondominated phase.During this phase, zonal flux stays within approximately 2-17 % of the maximum value reached during the time-dependent phase for most simulations.Meridional flux also gently fluctuates, being always positive in the upstream region and negative in the downstream for all model scenarios.The two surface-forced simulations (She and SF) show a zonal flux that continuously increases during the model simulation, with a final zonal flux that is approximately 50 % stronger than flux at the end of the time-dependent phase (Fig. 4b).This indicates that neither scenario may be reaching an advection-dominated phase.
Vertical flux time dependence varies between model simulations, with three primary patterns emerging.Firstly, vertical flux varies between positive and negative transport over both the upstream and downstream plane of the canyon, with flux values being roughly the same in the upstream/downstream region (Fig. 4a).This pattern is seen in the UW and OW simulations.Secondly, flux across the upstream plane is mirrored across the downstream plane, i.e. as the magnitude across one plane increases, the magnitude across the other plane increases well (Fig. 4c).In six simulations (ST, KL, SK, US, HB, SHR) this pattern occurs with upstream transport always positive and downstream transport always negative.In four runs (OBC, LRC, BLRB, KHRB) the above pattern occurs, but flux across the upstream/downstream planes does cross between positive and negative values during model days 4-6.Thirdly, simulations with a coastal jet (She and SF) exhibit strengthening negative flow across the upstream plane and weakening negative flow across the downstream plane (Fig. 4b).In the She case, flux across the downstream plane becomes positive around model day 4.In the SF case, flux over the upstream plane begins to weaken around model day 7.
To ensure that aliasing is not occurring with 12 h model output, another 10-day UW simulation was run with model output written every 3 h (Spurgin, 2014).Small differences in flux estimates are seen during the time-dependent phase.Differences during the advection-dominant phase are less than 10 %.

Circulation in the canyon
Model simulations exhibit three types of horizontal circulation: (1) formation of an anticyclonic eddy within the canyon, (2) cyclonic circulation everywhere within the canyon, and (3) weak circulation everywhere within the canyon.The evolution of the first two horizontal flow patterns is discussed below.Firstly, the anticyclonic circulation is detailed, followed by a description of the cyclonic circulation.
For high Burger number simulations (particularly UW, OW, OBC, ST, HB, KHRB) horizontal flow during the timedependent phase is cyclonic over the canyon (Fig. 5a, left).Toward the end of this phase, flow along the downstream rim becomes stronger relative to the upstream rim (Fig. 5b, left).After one more model day (by day 3), flow in the canyon head becomes anticyclonic, and this pattern persists for the remainder of the model simulation (Fig. 5c, left).Vertical velocity during the first day of simulation is negative everywhere in the canyon, and strongest in the upstream region.As the flow evolves, a region of positive vertical velocity appears in the downstream half of the canyon and moves toward the downstream corner of the canyon mouth.During the time-dependent phase, flow becomes faster over the canyon axis and impinges on the downstream wall.Strong downwelling occurs at and above shelf break depth (Fig. 5, left) but decreases with depth.This leads to compression of the isopycnals as they cross the canyon, and an anticyclonic eddy forms in the canyon, which persists during the advection-dominated phase (Fig. 6a).
For low Burger number simulations (particularly KL, SK, US, SHR), the cyclonic circulation that forms during the time-dependent phase strengthens as zonal flow accelerates and remains cyclonic (Fig. 5, right and 6b) as downwelling in these cases is quite uniform with depth.Similar to the cases with anticyclonic circulation, vertical flow is negative everywhere within the canyon and strongest in the upstream region.However, as the flow evolves, positive vertical velocity begins to appear in the downstream half of the canyon where it remains (it does not get pushed offshore) (Fig. 5, left).For these cases, horizontal flow is strongest along the canyon walls and weaker across the canyon axis (Fig. 6b).
For the purposes of this study, results focus on flow dynamics during the advection dominant phase.Results in the following sections are time-averaged model output.Based on oscillations in the vertical flux time series (Fig. 4), results are averaged from model days 5 to 8.

Comparison to previous studies
Current simulations reproduced canyon circulation features seen in previous studies (Table 1).Circulation features observed in the Klinck (1996) study were reproduced in the KL simulation.In this low Rossby number, low Burger number, and low |χ | simulation, the cyclonic circulation, antisymmetrical vertical velocity, and density change pattern (positive density anomaly outside the canyon mouth) seen in the original study was reproduced in the current model.Similarly, circulation features observed in Skliris et al. (2001Skliris et al. ( , 2002) ) were Circulation features observed in She and Klinck (2000) were mostly reproduced in the She simulation.In this low Rossby number, intermediate Burger number, and intermediate |χ| simulation, net downwards vertical velocity was reproduced.However, a weak cyclonic circulation was seen at all depths in the canyon.This is different from the original study, which saw anticyclonic circulation over the canyon and cyclonic circulation at 300 m and below.Multiple scenarios with varying Rossby, Burger, and |χ | values (UW, OW, OBC, ST, LRC, SHR) exhibit anticyclonic vorticity at shelf break depth and cyclonic vorticity deeper in the canyon, similar to that seen in She and Klinck (2000).Observations from Flexas et al. (2008) were reproduced in the UW and OW simulations.In these high Rossby number, high Burger number, and high |χ| simulations, the anticyclonic circulation and periods of net positive vertical velocity similar to that seen in the observations were replicated.

Vorticity in the canyon
As previously discussed, all model scenarios either form an anticyclonic eddy in the canyon after the time-dependent phase (Fig. 6a) or have cyclonic circulation in the canyon throughout model simulation (Fig. 6b).Looking at shelf break depth circulation, six simulations exhibit an anticyclonic eddy, four simulations show cyclonic circulation, and four other simulations show weak circulation at this depth.
Simulations with anticyclonic circulation display negative vorticity within the canyon, but opposing positive vorticity along the canyon walls (Fig. 7a).Simulations with cyclonic circulation have the opposite feature, i.e. positive vorticity within the canyon, but negative vorticity along canyon walls (Fig. 7b).This reversal of vorticity along bottom topography is due to friction between water parcels and canyon walls.The simulation in which the anticyclonic eddy appears only at a depth below the shelf break (HB) has negative vorticity at shelf break depth (Table 6).

Vertical velocity
Two main flow patterns are seen in plane views of vertical velocity.In the first pattern, enhanced downwards (upwards) velocity is confined to the upstream (downstream) corner of the canyon mouth at shelf break depth (Fig. 5c, left).This pattern occurs in simulations with weak or negative vorticity at shelf break depth (Table 6).One exception is the HB scenario, which exhibits varying positive and negative vertical velocity patterns everywhere within the canyon (Spurgin, 2014).
In the second vertical flow pattern, vertical velocity presents a more antisymmetric pattern, similar to that seen in previous studies (Klinck, 1996;Skliris et al., 2001Skliris et al., , 2002)).In these simulations, regions of positive and negative vertical velocity are split along the canyon axis from canyon head to canyon mouth, with negative velocity in the upstream region and positive velocity in the downstream (Fig. 5c, right).This pattern occurs in simulations with strong, cyclonic vorticity (Table 6).
Net vertical transport values averaged over the advectiondominated phase (model days 4-10) reveal two patterns across the upstream and downstream planes (Table 7).In six model simulations, transport is downwards (negative) across both shelf break planes.In three of these cases (OW, SF, LRC) upstream transport is larger than downstream  Units are 10 3 m 3 s −1 .U , V , W represent zonal, meridional and vertical transport values, respectively (Fig. 3).Negative values indicate downstream, offshore, and downwards transport, respectively.
transport.While the three other simulations (UW, BLRB, KHRB) show the opposite pattern, downstream transport is larger than upstream transport.Eight model simulations exhibit downwards transport across the upstream plane but upwards transport across the downstream plane.In all of these cases, downwards transport is stronger than upwards transport.Meridional (cross-shore) transport shows flow to be onshore in the upstream plane and offshore in the downstream (Table 7).Using the 12-hourly flux time series, periods of positive vertical flux across three planes (100, 150, and 600 m) in the canyon (i.e.everywhere between canyon walls from the canyon head to mouth) were calculated over the 10-day model period (Fig. 8).Net upward flux does occur in various downwelling canyon simulations.Note that only four scenarios do not see net upwelling at any time (She, SF, SHR, BLRB).Net upwards flux most commonly occurs across the 600 m plane, and least often across the 100 m plane.Periods of net upwards flux mostly occur during the advectiondominated phase.The longest occurrence of net upwards flux appears in the OBC simulation: this period lasts 2.5 model days.Overall, the OBC scenario exhibits the most occurrences of net upwards flux across the 150 and 600 m planes.

Density anomaly
Density anomalies were calculated as density variations (averaged over days 5-8) relative to initial density profiles.Similar to vertical velocity, two distinct anomaly patterns appear.In the first pattern, the density anomaly is negative everywhere in the canyon domain (Fig. 9a).At all depths, anomalies are strongest along bottom topography and weaken towards the offshore.This pattern is seen in simulations with a coastal or outer-shelf jet (UW, OW, OBC, ST, She, SF, BLRB; Table 8).
In the second pattern, there are strong negative anomalies along bottom topography, but there are also positive anomalies away from the canyon (Fig. 9b).This pattern is seen in simulations with a shelf break or offshore jet (KL, SK, US, HB, LRC, SHR, KHRB); the depth range of positive density anomalies varies (Table 8).For the majority of these simulations, positive anomalies do not extend down to shelf break depth (SK, US, HB, SHR, KHRB), but, in two simulations, positive anomalies do extend down to shelf break depth and below (KL, LRC).
Average density anomalies (averaged during the advection-dominated phase, model days 4-10) calculated in and around the submarine canyon are negative everywhere within (US_can and DS_can) and over the submarine canyon (US_shallow and DS_shallow), as well as along the shelf in the upstream (US_shelf) and downstream (DS_shelf) regions for all model simulations (Table 9).Four simulations (UW, OW, She, SF) exhibit negative density advection everywhere near the canyon (all nine regions).Eight simulations (ST, KL, SK, US, HB, LRC, SHR, KHRB) exhibit positive anomalies everywhere offshore of the canyon (Lower_can, US_off, and DS_off), while two simulations (OBC, BLRB) display a positive density anomaly only in the upstream and downstream offshore regions.

Density difference anomaly
As shown in vertical velocity and density difference, regions with downwelling canyons exhibit a background downwelling flow, with both negative vertical velocity and negative density changes away from the canyon region.To determine what extra effect a canyon has in a downwelling region, a density difference anomaly was calculated.Downwards density advection is enhanced in all canyon scenarios (Fig. 9c and d).Two patterns are seen; however, these patterns do not line up exactly with density anomaly patterns.Pattern 1 exhibits downwards density advection which is strongest in the canyon head and along the canyon axis (Fig. 9c); this pattern is seen in half of the simulations (UW, OW, OBC, ST, HB, SF, and LRC).Pattern 2 exhibits downwards density advection which is strongest along bottom topography (KL, SK, US, She, SHR, BLRB, and KHRB; Fig. 9d).Independent of density advection pattern, eight simulations exhibit positive density difference anomalies away from the canyon.In these cases, weaker downwelling (positive density difference anomalies) occurs along either the upstream (OW) or downstream (SF, LRC, and BLRB) corner of the canyon mouth, or along both corners of the canyon mouth (KL, US, She, and SHR).

Nitrate anomaly
All model simulations included a passive tracer (nitrate concentration) that was initialised with the same vertical variation for all runs.This provided one parameter that was the same in all simulations, allowing for better comparisons between model runs.Again, anomalies were calculated as changes between the initial nitrate profile and the model output nitrate profile averaged for days 5-8.
Similar to vertical velocity and density anomalies, two patterns appear.With the exception of two cases (OBC and ST), simulations that exhibit the pattern 1 density anomaly exhibit the same pattern for the nitrate anomaly, i.e. negative everywhere within the canyon and strongest along bottom topography (Table 8; Fig. 10a).Similarly, simulations that show pattern 2 density anomalies also exhibit negative nitrate anomalies along bottom topography and positive anomalies a few kilometres offshore (Table 8; Fig. 10b).The two simulations with open boundary conditions (OBC and ST) exhibit a pattern similar to pattern 2; however positive nitrate anomalies occur 5-20 km offshore.

Upwelling in downwelling canyons
There are various ways in which upwelling can be defined.Firstly, upwelling can be characterised as the net upwards movement of water in a region.Secondly, upwelling can be described as the net onshore movement of dense, cold (usually nutrient-rich) deep ocean water to the shallower coastal ocean.Coastal upwelling typically involves both processes working together, i.e. as surface waters are pushed offshore, deep ocean water is brought up from the depth to replenish surface waters along the coast.However, upwelling along the coast does not always occur following this same process.
During two observational cruises around BC, Flexas et al. (2008) used velocity and hydrographic samples to calculate vertical flux in the canyon.The authors estimated that, at approximately 100 m depth and shallower, vertical velocities were negative; below the thermocline (∼ 100-200 m depth), vertical velocities were positive.The authors concluded that upwelling did occur in BC, with a maximum near the shelf break depth and extending between 100 and 200 m.Ardhuin et al. (1999) modelled an upwelling cell beneath a trapped anticyclonic eddy.In their study, offshore deep waters from 300 to 500 m depth were lifted at the canyon wall and pulled out to the open ocean in the 200-300 m layer.   .Pattern 1, negative nitrate anomalies at all depths, corresponds to pattern 1, negative density anomalies at all depth (Fig. 9a); similarly, pattern 2, positive nitrate anomalies away from the canyon, corresponds to pattern 2, positive density anomalies away from the canyon (Fig. 9b).
Anomalies are averaged over a period of 3 model days during the advection phase (model days 5-8).
Plan-view images of time-averaged vertical velocity in the current simulations do not directly reveal net upward movement of water in any canyon scenarios.All of the simulations show regions of both upwards and downwards vertical velocity; however, downwards motion always appears to be dominant.Results of snapshots of net vertical velocity across three planes indicate that net upwards displacement of water does occur in the majority of simulations (Fig. 8).In most cases this upwards displacement is brief and commonly occurs at depth.However, the OBC case shows an extended period of net upwards velocity across the 150 m plane from model days 4 to 6, the beginning of the advection-dominated phase.This period of net upwards velocity at shelf break depth is comparable to observations in Flexas et al. (2008).In the current study, a period of net upwelling may be occurring, but the time-mean flow of the advection phase indicates an overall net downwelling.
The prevalence of upwards displacement across the 600 m plane indicates a possible upwelling cell similar to Ardhuin et al. (1999), in which deep water is upwelled along canyon walls but returns to the offshore before crossing shelf break depth.Increased stratification has been found to reduce vertical transport (Klinck, 1996).For all current simulations, stratification is weaker with depth, which is likely the reason why vertical exchange shows more variation at depth.The irregular occurrences of net upwards displacement across vertical planes, even under semi-steady circulation, indicate that observational studies may not be detecting the time-mean flow dynamics occurring in submarine canyons.
Though vertical velocity reveals that upward displacement of water does occur in downwelling canyons, density and nitrate anomalies indicate there is downwards advection of physical properties.Both density and nitrate exhibit regions of positive advection; however these regions occur away from the canyon and positive advection is weaker than negative advection.
Previous studies have found that submarine canyons in downwelling regions enhance coastal downwelling (Klinck, 1996;She and Klinck, 2000;Skliris et al., 2001Skliris et al., , 2002)).Anomalies of density difference in current simulations show that downwelling is enhanced everywhere within the canyon and over the lower canyon, with downwelling being strongest around the canyon axis.In the upper 100 m, relatively weak downwelling occurs over the mid-canyon.This region of weaker downwelling appears as a relative lifting of isopycnals from the downstream to upstream canyon rim (Fig. 11).Lifting isopycnals is often a characteristic of upwelling occurring in a region.However, these are instantaneous profiles of what is occurring in the canyon.Using profiles of density difference and density difference anomaly, it can be seen that this relative lifting of isopycnals is not actually upwelling but in fact a region of relatively weak downwelling.

Parameter effects
Stratification has been found to have significant impacts on vertical and cross-shore transport (Klinck, 1996).Previous studies have compared forcing responses between weakly and strongly stratified domains.However, these studies use either a uniformly stratified domain (Klinck, 1996) or a domain in which stratification varied in only three regions (Skliris et al., 2001(Skliris et al., , 2002)).Other studies have varied stratification in the canyon domain (Ardhuin et al., 1999;She and Klinck, 2000), but the effects of vertical variation in stratification have never before been studied.Thus, the nondimensional parameter χ was introduced in this study to investigate the impacts vertical changes in stratification have on flow dynamics.Weak χ values indicate stratification is more uniform in the domain.Negative χ values indicate there is stronger stratification variation at shallower depths.In this section, Rossby number, Burger number, χ values, and incoming jet location are used to determine which regional parameters impact various flow dynamics.

Circulation in the canyon
Patterns in horizontal circulation reveal the Burger number to be an important parameter in determining vorticity within the canyon (Fig. 12a).Four simulations with low Burger numbers (SHR, US, SK, KL) exhibit positive vorticity and  cyclonic circulation at shelf break depth.Six simulations with high Burger numbers (UW, OW, OBC, ST, HB, KHRB) exhibit negative vorticity and anticyclonic circulation near shelf break depth.Four simulations with varying Burger numbers (BLRB, She, SF, LRC) show weak vorticity and circulation.The Rossby number also appears to have an impact on vorticity magnitude (Fig. 12b), as for each vorticity sign and jet placement, vorticity magnitude increases with increasing magnitude.The importance of the Burger number and its impact on circulation within the canyon is highlighted in the US and HB cases.These two simulations have the same model setup, with the only difference being their buoyancy frequency (N) value.The US simulation has a lower Burger number and cyclonic shelf break circulation, while the HB simulation has a higher Burger number and anticyclonic shelf break circulation.
Anticyclonic circulation scenarios exhibit strong flow across the canyon axis and weaker flow along canyon walls and rims during the time-dependent phase, creating a negative shear in horizontal flow.In these simulations, flow is negligible along bottom topography (including the outer shelf) and flow turns weakly into the canyon but does not follow canyon isobaths.This causes flow crossing the canyon axis to impinge on the downstream wall, and a small portion moves onshore due to negative vorticity in the canyon.This onshore flow, as well as compressing isopycnals, generates a trapped anticyclonic eddy within the canyon (Fig. 6a), which can persist to depths of 500 m.Cyclonic circulation scenarios have weakest horizontal flow across the canyon axis and strongest flow along canyon walls and rims (Fig. 6b).This creates a positive shear in horizontal flow and enhances the tendency for cyclonic circulation.In these simulations, flow across the canyon axis is relatively weak and water parcels follow canyon isobaths, moving onshore in the upstream region and offshore in the downstream (Fig. 6b).
The location of the incoming zonal jet is affected by stratification variations in the domain (Fig. 13).Several model simulations use the same forcing conditions and different domain stratification, e.g.UW (SK) and US or HB (LRC) have the same forcing conditions with differences in χ values.The simulations with uniformly stratified domains (lower χ; US, HB, LRC) have a zonal jet that is located further offshore relative to their counterparts with high χ values (UW, SK).Simulations forced by wind stress have weak coastal flows, regardless of domain stratification.
As the body force is applied to model simulations, an onshore flow occurs and tilts isopycnals downwards.This leads to surface intensification of the zonal jet.Baroclinicity increases with increasing stratification in the upper water column.However, downwelling tends to reduce stratification over the shelf, making the jet more barotropic.With weak near-surface stratification (low χ ), this leads to an almost barotropic jet which feels bottom friction fairly strongly and is therefore reduced in intensity.In simulations with an offshore jet and low χ values, the zonal velocity is intensified at shelf break depth and negligible along continental slope topography (Table 4; KL, US, HB, LRC, KHRB).Strong nearsurface stratification (high χ) allows for baroclinicity of the jet on the shelf, and thus less friction on it.In simulations with an outer-shelf jet and high χ values, zonal velocity is intensified near the surface and negligible along bottom topography (Table 4; UW, OW, OBC, ST, BLRB).Until now, the Rossby number has been based on incoming flow strength.However, a Rossby number based on flow across the canyon axis (R U can ) may be more appropriate for looking at flow dynamics within the canyon.These values are compared in order to determine correlations between incoming and canyon axis flow (Fig. 14).Due to the more complex topography in the slanted canyon simulation (ST), it is suspected that R U can is overestimated for this scenario and is thus marked as a possible outlier in subsequent plots.
For simulations with cyclonic circulation at shelf break depth, there is a relatively strong coupling between incoming and canyon axis flow speed.However, for anticyclonic circulation the coupling is weaker.The same incoming Rossby number results in a weaker canyon Rossby number.For the simulations with weak circulation, both Rossby numbers are small and not strongly correlated with each other.
Flow patterns describe the stronger (weaker) coupling between incoming and canyon Rossby numbers in the cyclonic (anticyclonic) simulations.In the anticyclonic cases, the eddy is focused over the canyon axis and flow in the canyon head is towards the upstream, while flow along the mid-canyon is towards the downstream canyon wall (flow between mid-canyon and canyon mouth is downstream and slightly stronger along the canyon axis) (Fig. 6a).This causes net zonal flow between the canyon head and mid-canyon to be almost negligible and thus weaken overall flow strength across the canyon axis.For cyclonic simulations, zonal flow is weaker across the canyon axis but everywhere towards the downstream canyon wall (Fig. 6b).Therefore, strong incoming flow increases zonal flow everywhere within the canyon.
For all simulations, the Rossby number across the canyon is approximately one-third (or more) lower than the Rossby number based on incoming flow.The incoming Rossby number is based on maximum zonal flow at shelf break depth upstream of the canyon, whereas the Rossby number across the canyon is integrated from the canyon head to canyon mouth.Thus, the incoming Rossby number is slightly overestimated relative to the canyon Rossby number.Downwelling submarine canyons have been observed to modify incoming coastal jets by deflecting the current along canyon walls, with major modifications observed at shelf break depth (Flexas et al., 2008).Current simulations indicate two types of flow deflection occurring, depending on Burger number.With a low Burger number, flow follows canyon isobaths with strongest flow along bottom topography.With a high Burger number, the flow is more baroclinic and a cut-off anticyclonic eddy forms at shelf break depth in the canyon.A zero Burger number has been shown to generate almost symmetric cyclonic flow in previous simulations (Kämpf, 2006).

Vertical flux
Net vertical flux was calculated for all model simulations as net flux in the canyon across the shelf break plane (150 m) averaged during the approximate advection-dominated phase (averaged from model days 4 to 10).Initial errors in net vertical flux were calculated as the difference in flux values during two averaging periods: days 4-10 and days 3-9.However, two other sources of error were taken into consideration: (1) error due to variations in zonal flux (which varied by 2-17 % for most cases and 50 % for the She and SF simulations), and (2) 12 h model output provided an approximate 10 % aliasing error.Therefore, error in all model simulations was taken as the maximum error in (1) the sum of the minimum 10 % aliasing error plus error due to zonal flux variations or (2) errors in averaging period.
Net vertical flux is directly proportional to Rossby number of flow across the canyon axis and inversely proportional to Burger number (Fig. 15).For example, US and SHR have the highest ratios of canyon Rossby number to Burger number, and exhibit the greatest downward flux.Just as the strength of net vertical flux is affected by canyon Rossby number and Burger number, so too is the transport pattern (Table 7).Simulations with a ratio of canyon Rossby number to Burger number of approximately 0.1 or lower exhibit negative transport across both vertical planes.However, simulations with a ratio of canyon Rossby number to Burger number of approximately 0.1 or higher exhibit downwards flow in the upstream plane and upwards flow in the downstream.
Although previous studies have not specifically looked at changes in zonal flow strength, current simulations show that scenarios with stronger flow across the canyon axis lead to stronger downwards flux.This is unsurprising since increasing Rossby number indicates increasing cross-canyon flow, and thus a stronger pressure gradient along the canyon.
Increased stratification has been found to reduce vertical and cross-shore transport, as well as the depth range which fluid parcels move in a circuit around a canyon (Klinck, 1996;Skliris et al., 2001Skliris et al., , 2002)).Current simulations show a similar pattern.For example, US and HB cases have the same forcing and variation in stratification, with the only difference being that HB has an increased buoyancy frequency (N).Net cross-shore transport (not shown) in the US (weaker stratification) case is 2 times greater, and net vertical transport is approximately 6 times greater.
The upwelling flux through submarine canyons can be estimated using a scaling analysis based on the strength of flow, stratification, Coriolis parameter, and topographic shape parameters, including the slope of the continental shelf (Howatt and Allen, 2013).Using this scaling analysis, an approximate upwelling flux was estimated for all canyon scenarios in the present study and was compared to the downwelling flux found from the numerical simulations.The OBC simulation is an outlier because of the very small downwelling flux in this scenario (Fig. 15).Considering the other cases, for low Rossby numbers (She and SF) we find similar upwelling and downwelling fluxes, consistent with a large role of time dependence in these cases.As the Rossby number increases, the downwelling flux increases approximately linearly (Fig. 15).However, the upwelling flux increases much more quickly (Howatt and Allen, 2013) as it is an advectiondominated process.For the high Rossby number cases (UW, OW, ST, US, HB, KHRB), we estimate upwelling of the order of 10× the size of downwelling.As the upwelling and downwelling fluxes are measured differently, we cannot give exact values, but upwelling is clearly stronger than downwelling for Rossby numbers greater than about 0.04.
Circulation and vertical velocity are instantaneous measurements that capture what is occurring during the advection-dominated phase, and both are influenced by Burger and Rossby number.The Burger number drives circulation type and strength of vertical flux: simulations with a high Burger number exhibit a trapped anticyclonic eddy and weak vertical flux, simulations with a low Burger number exhibit cyclonic circulation and strong vertical flux.The Rossby number drives strength of vertical flux: high Rossby numbers generate strong vertical flux.Thus, Burger and Rossby numbers are important parameters during the advection-dominated phase.

Density anomalies
Patterns in density anomalies indicate that the Burger number (and subsequently vorticity) has the largest impact on the magnitude of density anomalies (Fig. 16).All simulations show net downwards advection of density within the canyon.Simulations with lower Burger numbers (cyclonic circulation) have weaker density anomalies at shelf break depth.Simulations with higher Burger numbers (anticyclonic circulation) exhibit stronger downwards density advection.Simulations with weak circulation also appear to be affected by Burger number.
It is unsurprising that simulations with high Burger numbers have stronger density anomalies.These simulations have stronger variations in initial density between vertical layers.Thus, a water parcel that advects the same vertical distance in a simulation with a high Burger number vs. one with a low Burger number will have a stronger density anomaly.The density anomaly was normalised by the Burger number effect, but showed no relationship to Rossby number, χ, or jet location.
Location of incoming zonal jet does affect the occurrence of positive density anomalies in the model domain (Table 8).Simulations with the zonal jet located along the shallow coast or over the outer shelf (directly above the upper canyon) exhibit negative density anomalies (downwards advection) at all depths in the model domain (pattern 1).Simulations which show upwards density advection (positive density anomalies) during the day 5-8 averaged period have a zonal jet located either offshore or along the shelf break.The positive anomalies occur away from the canyon in the upstream and downstream (pattern 2) regions.There is no correlation between pattern 1/2 and the average density anomaly over the whole canyon region.However, simulations with zonal jets that are coastal or on the outer shelf have negative anomalies over the lower canyon, whereas those with shelf break or offshore jets have positive anomalies there (Table 9).
Forces can explain how jet location impacts the occurrence of positive density anomalies.For the zonal jet to turn shoreward along the upstream canyon rim, it needs a centrifugal force.This is provided by a change in the pressure gradient: higher pressure offshore of the jet and lower pressure onshore.This weakens the Coriolis force and pressure gradient force balance and allows flow to turn.In the simulations with an offshore or shelf break jet, this pressure gradient change is provided by upwelling of denser water occurring offshore of the canyon and zonal jet (Fig. 9b).For the simulations with a coastal or outer-shelf jet, this upwelling occurs over the outer shelf, where stronger downwelling is already occurring.This is seen as a reduction of downwelling rather than upwelling (Fig. 9a).
Previous studies have found that submarine canyons, in regions with a right-bounded jet, enhance the downward advection of density properties (Klinck, 1996;She and Klinck, 2000;Skliris et al., 2001Skliris et al., , 2002)).Similar results are seen in plots of density difference anomalies (Fig. 9c and d).Downwards advection of density is enhanced within the canyon for all simulations.

Nitrate anomalies
Comparisons of non-dimensional parameters indicate incoming Rossby number has the greatest impact on vertical advection of nitrate (Fig. 17 10), the nitrate anomaly is weaker and broader in the stronger stratification scenario (UW) relative to the weaker stratification scenario (US).The region of negative anomalies upstream of the canyon is approximately 2.5-3 times larger in the weaker stratification case (UW).The Rossby radius of deformation is 3 times larger in the UW case, so size of the nitrate anomaly is proportional to Rossby radius (and stratification).Therefore, the strength of the nitrate anomaly is inversely proportional to stratification, while the size of the nitrate anomaly is directly proportional to stratification.These have a cancelling effect, and only the Rossby number appears to have an influence on the nitrate anomaly.
With the exception of the two open boundary simulations, density and nitrate exhibit the same anomaly patterns in the same model scenarios (Table 8): negative anomalies everywhere in the canyon domain (pattern 1) and positive anomalies away from the canyon (pattern 2).This indicates that jet location impacts the occurrence of nitrate anomalies following the same reasoning described in the previous section.Therefore, the upwelling occurring away from the canyon includes denser water and higher nitrate concentrations.
Density and nitrate anomalies are integrated measurements and include a strong signal from the time-dependent phase.Anomaly patterns are influenced by jet location, which in turn is affected by vertical variations in stratification.As previously discussed (vertical velocity section, Sect.3.5.1),zonal jet placement effects time dependence of vertical flux.Offshore and shelf break jets generate upward velocity away from the canyon and steadier vertical flux (relative to coastal/outer-shelf jets).Thus, the location of the zonal jet is important during time-dependent phases of flow.
Average diapycnal diffusivity in Ascension Canyon (west coast, North America) has been observed as approximately 3.92×10 −3 m 2 s −1 (Gregg et al., 2011).Diffusion in the current model is small (10 −7 m 2 s −1 ); thus mixing has an insignificant impact on nitrate flux calculations in the model simulations.Comparisons of estimated diffusive flux and model-calculated advective flux provide true comparisons of the separate processes (Spurgin, 2014).In the upper 100 m, diffusion of nitrate is stronger than mean advection of nitrate.At 150 and 600 m, mean downwards advection of nitrate is stronger than upwards diffusion of nitrate.These positive diffusive flux estimates indicate that nitrate anomalies calculated in the previous section (Fig. 17) are likely stronger than what would occur in a real-world downwelling canyon scenario.

Downwelling dynamics and biological productivity
Previous observational studies in the NW Mediterranean Sea have shown that physical transport processes affect planktonic and particle distributions within and around submarine canyons (Alvarez et al., 1996;Granata et al., 1999), respectively.In Palamòs Canyon, an intruding filament of cold, salty oceanic water at 50 m depth correlates to a high density distribution of planktonic larvae at the same depth (Alvarez et al., 1996).In BC, concentrations of total particulate matter are highest in downwelling zones, particularly within anticyclonic cores (Granata et al., 1999).One numerical model coupled a hydrodynamic model and a coastal plankton ecosystem model to further investigate the impacts submarine canyons have on primary production www.ocean-sci.net/10/799/2014/Ocean Sci., 10, 799-819, 2014 (Skliris and Djenidi, 2006).Upwelling of deep water rich in nitrate upstream and downstream of the canyon was found to enhance primary production, while cyclonic circulation in the canyon led to an accumulation of plankton biomass in the canyon (ibid.).All canyon simulations exhibited net downwards nitrate advection, which suggest that, in nitrate-limited regions, steady downwelling canyons lead to a reduction of primary productivity in the region of canyons.Net downwelling occurring in the canyon will focus sinking particulate matter.Reasonably strong vertical flows, particularly on the downstream side of canyons, could lead to a concentration of vertically migrating zooplankton (Ianson et al., 2011).Particle tracks (Spurgin, 2014) showed that simulations with anticyclonic circulation (UW, OW, OBC, ST, HB, KHRB) produced a looping flow for particles, suggesting that, in these simulations, the canyon is generally a retention region.

Summary and conclusions
Our studies have shown that the Burger number (stratification) has the largest impact on flow dynamics in downwelling submarine canyons: (1) cyclonic circulation (positive vorticity) occurs in canyons with low Burger numbers, and (2) anticyclonic circulation (negative vorticity) occurs in canyons with high Burger numbers.Next in importance is the Rossby number; an increasing Rossby number generally increases the magnitude of the vorticity.Weaker stratification (low Burger number) and stronger flow (high Rossby number) lead to greater vertical flux at shelf break depth.Jet placement is of third importance and it is partially determined by the variation in stratification with depth.It is important in density and tracer vertical advection, both of which are mostly determined by the time-dependent phase of the flow.Jet placement impacts the occurrence of positive density/nitrate anomalies away from the canyon with offshore or shelf break jets leading to positive density anomalies away from the canyon.The strength of density anomalies is determined by the Burger number, whereas the strength of the nitrate anomalies is determined by the Rossby number.
Flow dynamics seen in present studies have been found in previous literature (Table 1) using similar forcing conditions.Klinck (1996) and Skliris et al. (2001Skliris et al. ( , 2002) ) had offshore/shelf break zonal jets with small Rossby and Burger numbers.This leads to a cyclonic flow pattern and small patches of weak upwelling away from the canyon.Blanes Canyon has an outer-shelf jet with high Rossby and Burger numbers (Flexas et al., 2008).The high Burger number and jet placement lead to an anticyclonic flow pattern and downwards density advection everywhere.The anticyclonic vorticity leads to a weak coupling between incoming flow strength and flow strength across the canyon axis.Vertical flux is weak and density (nitrate) flux is strong due to the high Burger (Rossby) number.She and Klinck's (2000) jet was near the coast and weakly coupled to the canyon, and thus vorticity, flux, density, and nitrate advection were weak.Our simulations have reproduced the major features of previous studies, including the differences between them.The dynamical explanation gives the reasons for those differences and provides an encompassing explanation of flow dynamics in downwelling canyons.However, this study also illustrates the strong response over canyons to time-varying flow -an area that should see further research in future studies.

FlexasFigure 1 .
Figure1.Canyon bathymetry (grey lines) and reference terminology used in this thesis.Canyon rim indicates the boundary between the continental shelf and the canyon.Shelf break indicates change in slope gradient between the continental shelf and slope.Canyon mouth is the open region along the shelf break, canyon head is the shallowest onshore canyon region, mid-canyon is the region between the canyon head and mouth, and lower canyon is offshore of the canyon mouth.Canyon wall refers to canyon topography between shelf break and bottom depth.Zonal flow is in the alongshore direction, and meridional flow is in the cross-shore direction.Contour intervals are 100 m.

Figure 2 .
Figure 2. Initial (a) density and (b) nitrate profiles for all model simulations.In the density profiles, grey represents SK and SHR; green represents KL and US; light green represents HB, LRC, and KHRB; light blue represents BLRB; brown represents She; and blue represents UW, OW, OBC, ST, and SF.Simulation abbreviations are explained in Sect.2.4.
negligible along bottom topography intensified over mid-outer shelf to shelf break OW outer shelf surface intensified; negligible along bottom topography intensified over mid-outer shelf to shelf break OBC outer shelf surface intensified; negligible along bottom topography; secondary jet at shelf break intensified over shelf break ST outer shelf surface intensified; negligible along bottom topography intensified over mid-outer shelf to shelf break KL offshore uniform offshore uniform offshore of shelf break; weak flow over flat shelf (shelf break to coast) along bottom topography (150-600 m) intensified between shelf break and 10 km offshore BLRB outer shelf surface intensified; negligible along bottom topography intensified over mid-outer shelf to shelf break KHRB offshore uniform offshore; negligible along continental slope uniform 10 km offshore of shelf break; weak flow over outer shelf; negligible over inner shelf

Figure 3 .
Figure 3. Planes used for transport calculations.

Figure 4 .
Figure 4. Time series of horizontal and vertical flux directly over the canyon for (a) core (UW), (b) surface-forced(She), and (c) slope-current (SK) simulations.U 3 indicates zonal flux across the canyon axis (from canyon head to canyon mouth).V 2 and V 3 are meridional flux along the canyon mouth from the upstream rim to canyon axis, and from the canyon axis to downstream rim, respectively.All horizontal fluxes are measured from surface to shelf break depth.W 2 and W 3 are vertical flux across the shelf break depth plane (150 m) everywhere within the canyon, from the upstream rim to canyon axis and from the canyon axis to downstream rim, respectively.Negative U , V , and W values indicate westward, offshore, and downwards fluxes.

Figure 5 .
Figure 5. Horizontal and vertical circulation at shelf break depth during the time-dependent phase on (a) day 1.5, (b) day 2.5, and (c) day 3.5.Pink shading indicates downwards velocity and blue shading indicates upwards velocity.Circulation for simulations with anticyclonic circulation (left) and cyclonic circulation (right) are shown.

Figure 6 .
Figure 6.Horizontal velocity vectors at shelf break depth (150 m) in (a) simulation with anticyclonic circulation (UW) and (b) simulation with cyclonic circulation (SHR).Vectors are averaged over a period of 3 model days during the advection phase (model days 5-8).

Figure 7 .
Figure 7. Cross section of vorticity at mid-canyon in (a) a simulation exhibiting negative vorticity (UW) and (b) a simulation exhibiting positive vorticity (SK).Values are averaged over a period of 3 model days during the advection phase (model days 5-8).

Figure 8 .
Figure 8. Times of net upwards flux across three vertical planes for all model simulations.Net upwards flux is plotted if larger than a minimum value of 1000 m 3 s −1 .

Figure 9 .
Figure 9. Density anomaly at shelf break depth (150 m) in (a) simulation with pattern 1, i.e. negative density anomalies at all depths (UW), and (b) simulation with pattern 2, i.e positive density anomalies away from the canyon (US).Density difference anomaly at shelf break depth (150 m) in (c) simulation with pattern 1, i.e. enhanced density advection everywhere in the model domain (UW), and (d) simulation with pattern 2, i.e. enhanced downwelling within the canyon and weaker downwelling near the canyon mouth (US).Anomalies are averaged over a period of 3 model days during the advection phase (model days 5-8).

Figure 10 .
Figure10.Nitrate anomaly 10 km upstream of canyon axis (left) and along canyon axis (right) in (a) a simulation with pattern 1, i.e. negative nitrate anomalies at all depths (UW), and (b) a simulation with pattern 2, i.e. positive nitrate anomalies away from the canyon (US).Pattern 1, negative nitrate anomalies at all depths, corresponds to pattern 1, negative density anomalies at all depth (Fig.9a); similarly, pattern 2, positive nitrate anomalies away from the canyon, corresponds to pattern 2, positive density anomalies away from the canyon (Fig.9b).Anomalies are averaged over a period of 3 model days during the advection phase (model days 5-8).

Figure 11 .
Figure 11.Isopycnal cross section at mid-canyon in the core simulation (UW).

Figure 12 .
Figure 12.Effect of (a) Burger number and (b) incoming Rossby number on canyon vorticity for all model simulations.

Figure 13 .
Figure 13.Zonal jet location upstream of the canyon.

Figure 14 .
Figure 14.Correlation between Rossby numbers based on incoming zonal flow and zonal flow integrated across the canyon axis.Due to complex canyon topography, canyon Rossby number for ST is likely overestimated and ST is considered an outlier.

Figure 15 .
Figure 15.Burger and canyon Rossby number effect on vertical flux at shelf break depth (150 m).Due to complex canyon topography, canyon Rossby number for ST is likely overestimated and ST is considered an outlier.

Figure 16 .
Figure 16.Burger number effect on average density anomaly in canyon for all simulations (averaged between model days 4 and 10).

Figure 17 .
Figure 17.Incoming Rossby number effect on nitrate anomaly for all simulations (averaged between model days 4 and 10).
). Simulations with higher Rossby numbers have greater changes in nitrate concentration.This indicates that advection of passive tracers strengthens as more flow enters a canyon.Patterns in incoming Rossby number show a stronger correlation than patterns based on canyon Rossby number.Density and nitrate anomalies measure a combination of time-dependent and advection-dominated downwelling.Density and nitrate anomaly time series for the UW case (not shown) indicate averaged anomalies are approximately 85 %

Table 1 .
Features seen in previous studies.

Table 2 .
Constant geometric parameters for model simulations.

Table 3 .
Non-dimensional parameters for all model simulations.

Table 5 .
Temporal variations in forcing for all model simulations.

Table 6 .
Absolute vorticity and canyon circulation for all model simulations.Values are averaged during the advection-dominated phase (model days 4-10).
* Absolute vorticity is taken as maximum vorticity in the canyon head, away from canyon rims, at shelf break depth.

Table 7 .
Average net transport values across shelf break depth for all model simulations.Values are averaged during the advectiondominated phase (model days 4-10).

Table 8 .
Positive density and nitrate anomaly depth range for all model simulations based on advection-dominated phase average (model days 4-10).Pattern 1 anomalies had no positive values.Pattern 2 anomalies had positive values.Positive anomalies are included if their value is greater than 10 % of the negative anomalies in the same horizontal plane.

Table 9 .
Average density anomalies across nine regions for all model simulations.Values are averaged during the advection-dominated phase (model days 4-10).