The impact of shelf slope on the linear stability of buoyant coastal currents and on the nonlinear formation of coastal meanders and eddies is investigated. The authors consider a simplified two-layer stratification in cylindrical geometry where a buoyant surface current flows along the coast above a denser water, with a flat bottom or steep shelves. Simulations were performed using the Nucleus for European Modelling of the Ocean (NEMO) ocean global circulation model. The initial state of these simulations was defined according to laboratory experiments performed in the same configuration. Comparisons between laboratory and numerical results highlight the role of momentum diffusion and of the initial perturbations amplitude. The authors’ results confirm that the topographic parameter To (ratio between the shelf slope and the isopycnal slope of the current) is the relevant parameter to quantify the shelf impact on the linear and nonlinear dynamics of the surface current. When the evolution of the buoyant coastal current is controlled by the baroclinic instability, the increase of To yields a selection of smaller unstable wavelengths and a decrease of the unstable growth rates. For finite values of To, a complete stabilization of the surface current can be reached. The typical radius of the first eddies generated by the coastal current is set by the linear stage of the baroclinic instability. However, secondary nonlinear processes may lead to larger or smaller structures. The authors exhibit a new dynamical sequence, leading to the formation of submesoscale cyclonic eddies over a steep shelf by splitting of mesoscale eddies. These cyclonic eddies trap and transport water masses and may play an important role in the cross-shelf exchanges.
Coastal currents are important features of the regional circulation that control the cross-shelf transport. However, coastal current is a generic term that covers a wide variety of dynamical configurations. This study focuses on buoyant coastal currents. Such geostrophic currents are characterized by a light water mass flowing along the coast above a denser water mass and by an outcropping density front located at the offshore edge of the flow, as shown in Fig. 1. Strait connections between distinct ocean subbasins are the main sources of buoyant coastal currents. For instance, in the Mediterranean Sea, the light Atlantic Water entering through the Strait of Gibraltar forms the Algerian Current in the western Mediterranean (Millot 1987; Obaton et al. 2000) and the Lybio–Egyptian Current in the eastern Mediterranean (Alhammoud et al. 2005; Hamad et al. 2005; Millot and Taupier-Letage 2005b). Other buoyant coastal currents in regional seas can be observed: the Norwegian Atlantic Current along the eastern coast of Greenland (Pickart et al. 2005) or the Bransfield Current along the southern coast of the South Shetlands Islands in the Bransfield Strait (Sangrà et al. 2011). Such currents generally flow over the coastal shelf and the bottom bathymetry has a significant impact on the current dynamics.
Coastal current instabilities may form meanders and lead to the formation of eddies. Because of complex interaction processes like air–sea interactions or bathymetric effects, a large range of unstable wavelengths is possible. Therefore, various cyclonic or anticyclonic eddies differing in size and intensity with various thermohaline characteristics can be generated. These coherent structures can either flow along or across the current or be detached from the coast (Millot and Taupier-Letage 2005a; Jouanno et al. 2008; Carton and Chao 1999). Hence, the coastal eddies play a significant role in the local mixing of biogeochemical properties and in the dispersion of pollutants and the redistribution of nutrient-rich coastal waters toward the oligotrophic open sea (Riandey et al. 2005). The numerical simulations of these eddies, in realistic configuration and without assimilation, is a major challenge. Actually, it is very difficult to forecast, in a regional model, the right eddy at the right location according to the various processes involved, particularly the bathymetric effect.
At a first order of approximation, we can simplify the vertical stratification of a buoyant boundary current as a two-layer system including a light water flowing above a dense bottom water. Hence, the stability of buoyant coastal currents with a flat bottom (Fig. 1a) has often been studied using two-layer models. One of the first attempts to describe the baroclinic instability was made by Phillips (1954) using a simplified two-layer quasigeostrophic model. The most unstable wavelength λB of this idealized baroclinic flow (constant velocity in each layer) corresponds roughly to λB ≃ 2πRd with being the baroclinic deformation radius (Pedlosky 1987; Vallis 2006), where h1 and h2 are the upper- and the lower-layer thicknesses, respectively; g* is the reduced gravity; and f is the Coriolis parameter. The baroclinic instability, due to a resonant interaction between two Rossby waves, is controlled for this case by the two-layer aspect ratio δ = h1/(h1 + h2). When the coastal current tends to be surface advected (i.e., δ goes to low values), the values of the growth rates of baroclinic instabilities decrease and the most unstable wavelength increases. The two-layer shallow-water model (Boss et al. 1996; Sakai 1989; Gula and Zeitlin 2010b) was then used to take into account ageostrophic instabilities. Unlike the intermediate models (i.e., quasigeostrophic or frontal–geostrophic models; Swaters 1993; Reszka and Swaters 1999), the shallow-water model takes into account finite Rossby numbers and fast wave motions. Hence, several new branches of instability may appear because of the unstable resonance between a geostrophic Rossby mode and an ageostrophic Kelvin or gravity wave. These ageostrophic instabilities generate unstable perturbations at smaller scales than the standard baroclinic instability does. The unstable wavelengths of a Rossby–Kelvin λRK or a Rossby–Poincaré λRP interaction are close to the deformation radius λRK ~ λRP ≃ Rd (Sakai 1989; Gervasio 1997; Gula and Zeitlin 2010a,b), which is 5–6 times smaller than the standard baroclinic wavelength selection λB ≃ 2πRd. However, these ageostrophic instabilities have large growth rates only for finite Rossby (or Froude) numbers (Sakai 1989; Gula and Zeitlin 2010b) and they are generally neglected for small Rossby number flows. The validity of quasigeostrophic models to describe unstable modes of outcropping fronts having large isopycnal deviation was studied by Boss et al. (1996). Their linear stability analysis shows that the spatial structure of the frontal modes induced by the outcropping front differs from standard Rossby modes but does not change the characteristics of the low wavenumber instability. Both two-layer quasigeostrophic and the shallow-water models predict the same growth rates for the unstable interactions of Rossby–Rossby or Rossby–frontal modes. Hence, the standard baroclinic instability is expected to be the dominant instability of small Rossby number buoyant coastal current.
Attempts to classify the dynamical interaction of a buoyant coastal current with the shelf slope have been previously undertaken. A classification scheme, which segregates surface-advected current over a steep shelf slope (Fig. 1b) from bottom-trapped current over a gentle shelf slope (Fig. 1c), was proposed by Yankovsky and Chapman (1997), Avicola and Huq (2002), and Lentz and Helfrich (2002). When the coastal current is bottom trapped, most of the light water is in contact with the shelf slope and both its width and stability are controlled by the bottom Ekman layer dynamics. According to numerical studies (Chapman and Lentz 1994; Yankovsky and Chapman 1997) or laboratory experiments (Whitehead and Chapman 1986; Lentz and Helfrich 2002; Wolfe and Cenedese 2006), the effects of bottom and lateral frictions and of the bottom Ekman circulation tend to widen and stabilize the bottom-trapped buoyant coastal current in comparison with the surface-advected configuration.
The present paper focuses on the baroclinic instability of surface-advected buoyant coastal current above a steep shelf slope. We emphasize that, for steep shelf configurations, the shelf effect cannot be neglected and has a significant impact on the linear stability of the surface-advected buoyant coastal current and on the nonlinear formation of large meanders and eddies along the coast. Indeed, a bottom slope affects the growth rates and the wavelengths of the most unstable baroclinic modes (Blumsack and Gierasch 1972; Mysak 1977; Mechoso 1980; Gervasio 1997; Lozier and Reed 2005; Isachsen 2011). When the bottom slope is positive (i.e., same sense as the isopycnal tilt), the potential vorticity (PV) gradient may vanish in the bottom layer and suppress the baroclinic instability in agreement with the Charney–Stern theorem (Pedlosky 1987). These theoretical results are supported by observations indicating that meanders do not grow upstream of the Cape Hatteras where the Gulf Stream flows over a steep continental slope. On the other hand, for buoyant coastal currents the bottom slope is negative (i.e., shelf slope and isopycnals tilt in the opposite sense). The impact of such negative slopes on the stability of coastal current is still under discussion. In the framework of quasigeostrophic models, both two-layer model (Mysak 1977) and the continuously stratified Eady model (Blumsack and Gierasch 1972; Mechoso 1980; Isachsen 2011) show that a negative shelf slope reduces the unstable growth of baroclinic modes. These idealized studies demonstrate that the central parameter of the problem is the ratio of the bottom slope over the isopycnal slope. However, these quasigeostrophic models are oversimplified and their predictions may not be valid for steep slope configurations, outcropping front and ageostrophic current. Hence, recent studies generally used the hydrostatic primitive equations to model the unstable dynamics of coastal current over sloping bathymetry. In this context, the linear stability analysis of Lozier and Reed (2005) shows that a negative shelf slope may amplify the unstable growth. Other works, using primitive equations simulations, study the eddy tracer transport across sloping bottom (Spall 2004; Isachsen 2011). According to these fully nonlinear simulations, devoted to thermally forced marginal sea, negative bottom slopes reduce the eddy diffusivity, in agreement with the stabilization predicted by the linear Eady theory (Blumsack and Gierasch 1972). Conversely, in an idealized model of the Nordic Seas, Spall (2010) shows that an increase of the topographic slope may result in an increased eddy flux. Besides, in recent laboratory experiments (Rivas et al. 2005; Wolfe and Cenedese 2006), the steep slope configuration seems to amplify the meander formation compared to the gentle slope configuration. Nevertheless, for such laboratory experiments, the combined effects of the vertical dissipation and the shelf slope can hardly be distinguished. According to these various approaches, several dynamical processes could be affected by the bottom slope variations. Thus, the impact of negative slopes on coastal fronts still leads to contradictory results in the literature.
To better understand how a steep shelf impacts on the stability of a buoyant coastal current we performed several numerical simulations using the Nucleus for European Modelling of the Ocean (NEMO) ocean general circulation model (Madec 2008) in an idealized two-layer configuration. The choice of this model is motivated by future modeling of the regional circulation in the Mediterranean Sea at a submesoscale ( and ) resolution. In section 2, we present the dynamical parameters governing this idealized configuration and the numerical parameterizations. In a first step (section 3), we study the shelf slope effect on a highly dissipative case equivalent to small-scale laboratory experiments. We then focus on low dissipative cases to study the impact of the shelf slope only (section 4). The baroclinic nature of the unstable coastal front is analyzed in section 4a. We then quantify the impact of steep slopes on the unstable surface current (growth rate and wavelength selection) and identify the relevant parameters describing the bathymetric influence in section 4b. The nonlinear evolution of meanders and eddies formed along the shelf and their role on horizontal transport are discussed in section 4c. Conclusions are given in section 5.
2. Idealized configuration
a. Physical parameters for a circular two-layer configuration
As a first approximation of a buoyant coastal current, we used an idealized two-layer configuration in a circular basin (Fig. 2). This dynamical configuration results from standard rotating lock release experiments (Griffiths and Linden 1980, 1981; Bouruet-Aubertot and Linden 2002; Stegner et al. 2004; Rivas et al. 2005) performed in a cylindrical tank. A fixed volume of light water of density ρ1, which flows above a dense water of density ρ2, is confined along the coast between the external cylinder of radius R2 and the density front, for which the surface outcropping is located at the radius R1 (Fig. 2a). Thus, this buoyant water mass is characterized by a width L = R2 − R1 ≈ 5Rd – 6Rd, which is larger than the baroclinic deformation radius Rd. For a surface-advected coastal current, the vertical aspect ratio of the flow δ=h1/(h1 + h2) is small. In the present study, this parameter remains at the fixed value δ = 0.15, and the baroclinic deformation radius is approximately .
The Burger number,
quantifies the ratio between the kinetic energy (KE) and the potential energy for a flow in geostrophic balance. For a small Burger number (i.e., L larger than Rd), the available potential energy (APE) is much larger than the kinetic energy of the surface flow. In the present study, the Burger number is kept constant in all the experiments (Bu = 0.02).
The width L of the light water mass, confined along the coast, may not coincide with the jet width W induced by the isopycnal tilt. Actually, here the jet width is W ≈ 2Rd (full width at half maximum in Fig. 3). The radial velocity profile satisfies the geostrophic balance in both layers and maximum velocity values are localized at the front (Fig. 2b). Both laboratory experiments (Bouruet-Aubertot and Linden 2002; Stegner et al. 2004; Thivolle-Cazat et al. 2005) and in situ measurements (Avicola and Huq 2002; Obaton et al. 2000) show that the jet width W scales with the baroclinic deformation radius Rd. Hence, we define the Rossby number as
where U is the maximum horizontal velocity of the buoyant coastal current. In this case, the Rossby number corresponds to a Froude number , which measures the ratio of the maximum current velocity U to the maximum phase speed of internal gravity waves . Finite Froude numbers lead to hydraulic jump formation and significant wave breaking. In real oceanic configuration, this parameter is generally small and we used in the following Ro = Fd ≃ 0.35.
The isopycnal tilt is quantified here by the maximum isopycnal slope α of the surface front. We use positive values for an anticlockwise slope direction. We then define a topographic parameter To as the ratio of the shelf slope s to the isopycnal slope α,
Positive values of To are obtained when both isopycnal (α) and shelf (s) slopes have the same direction. These cases correspond to typical upwelling events along the coast or to westward boundary currents. On the other hand, buoyant coastal currents correspond to negative values of To when isopycnal and shelf slopes are in opposite directions (Fig. 2a). Unlike previous works (Lozier and Reed 2005; Rivas et al. 2005; Wolfe and Cenedese 2006), where only the shelf slope s was used to quantify the bathymetric effect, we used in this paper the topographic parameter (i.e., a relative slope parameter). This choice is motivated by the Charney–Stern criterion, which indicates that stabilization occurs when the ratio of the isopycnal slope to the bottom slope becomes larger than unity (i.e., To ≥ 1) for a two-layer Phillips model with topography (Mysak 1977). An identical topographic parameter has been introduced by Blumsack and Gierasch (1972), Mechoso (1980), and Isachsen (2011) for Eady-type models with slopping bottom. Unlike the two-layer Phillips model, the bottom slope only enters as a boundary condition in the Eady problem and does not introduce a PV gradient in the interior. This relative bottom slope parameter controls both the growth rate and the unstable wavelength selection of the baroclinic flow. Hence, the topographic parameter seems to play a crucial role on the coastal flow stability and, in what follows, we will mainly vary To while keeping constant the other parameters (δ, Bu, and Ro). To keep constant the vertical aspect ratio δ at the density front location (r = R1) when the shelf slope is changed, the total open water depth H is adjusted. Note that, when the bottom topography is not flat, the total open water depth H is bigger than the total water depth h1 + h2 above the front at r = R1 (Fig. 2a).
Two dimensionless numbers are used to characterize the diffusion of momentum: the Reynolds number Re on the horizontal and the Ekman number Ek on the vertical. They are defined by
with Ah being the horizontal diffusivity coefficient, and
with Aυ being the vertical diffusivity coefficient. In the following, we will discuss different simulations characterized by high (low) Reynolds numbers corresponding to low (high) diffusion. For viscous laboratory experiments, the vertical (Aυ) and horizontal (Ah) diffusivity coefficients are equal to the molecular viscosity Aυ = Ah = ν. However, for high Reynolds number simulations, the vertical and horizontal “turbulent diffusion” (i.e., the Austausch coefficients) can differ.
b. Laboratory experiments
A few laboratory experiments, carried out with a small-scale setup (R1 = 12.5 cm, R2 = 20 cm, and h1 = 1.5 cm), initiate our work and motivate the full numerical study presented in this paper. The experimental setup is similar to standard rotating lock release experiments (Griffiths and Linden 1981; Bouruet-Aubertot and Linden 2002; Rivas et al. 2005). A fixed volume of light water ρ1, initially confined between a bottomless cylinder (r = R1) and the external boundary of a cylindrical tank (r = R2), is quickly released in a denser fluid ρ2. In contrast with other studies (Thomas and Linden 2007; Obaton et al. 2000; Lentz and Helfrich 2002; Helfrich and Mullarney 2005; Wolfe and Cenedese 2006), the coastal current generated here by gravitational collapse has no starting or ending point, because of the azimuthal symmetry, and therefore no boundary effect is present.
A two-layer salt stratification allows us to fix small values for the reduced gravity g* = g(ρ2 − ρ1)/ρ2 ≃ 0.005 g and to control the rotation speed Ω0 = 8 rpm to adjust the baroclinic deformation radius to . In this experiment, the surface layer thickness h1 ≃ 1.5 cm is of the same order of magnitude as the deformation radius Rd, and the isopycnal slopes α are therefore close to unity.
Standard particle image velocimetry (PIV) was used to measure the horizontal velocity field. Small buoyant particles were put in the upper layer and illuminated by a horizontal laser sheet of wavelength 670 nm, located a few millimeters below the upper free surface. The particle motion was recorded by a 768 × 576 pixel charge-coupled device (CCD) camera rotating with the tank. This camera had a resolution of 15 pixels per cm. The particle velocities were then estimated using LaVision PIV software with successive cross-correlation boxes (Teinturier 2010) yielding a final 63 × 64 or 48 × 49 vector field. Hence, the horizontal velocity field was measured in the surface layer with a resolution of about 2 vectors per deformation radius.
Two laboratory experiments were performed, the experiment EXP0 with a flat bottom and the experiment EXP1 with a bathymetry. Their detailed parameters are given in Table 1.
c. Initial state of the numerical experiments
We consider here an initial value problem. The dynamical evolution and the baroclinic stability of the buoyant coastal current strongly depend on the initial state of the experiment. In the laboratory experiments, after about one rotation period T = 2π/Ω0, the density front reaches a mean adjusted state (mean flow averaged over T) in agreement with the geostrophic or gradient wind balance. The transient and three-dimensional (3D) instabilities that occur during the very first stage of the adjustment (typically the first rotation period) are an efficient mechanism of turbulent dissipation at small scales in the frontal region (Stegner et al. 2004; Stegner 2007). It is then difficult to reproduce these small scales and 3D instabilities with a numerical model. Therefore, for the initial state of the numerical simulations, we defined a geostrophically balanced state using the laboratory measurements. To construct the initial three-dimensional velocity and density fields, we used the mean velocity profile estimated from PIV measurements of the laboratory experiment EXP0 (Fig. 3). This mean velocity profile corresponds to the azimuthal velocity spatially averaged over the circular basin and to a temporal averaging over one rotation period to filter out inertial oscillations. We first fit this mean surface velocity profile with the following analytical function:
This analytical function was constructed as a combination of the Rossby adjustment solution for a uniform potential vorticity front (Flierl 1979), with a hyperbolic tangent regularization of the velocity discontinuity with the scale δx. The values of R1 = 11Rd and δx = 1Rd were tuned to maximize the correlation. Because of the small vertical aspect ratio δ = 0.15, we assume that the velocity of the deep part of the bottom layer is negligible and set its value to zero. Here again we use a hyperbolic tangent profile to mimic the continuous velocity shear from the thin upper layer to the deep lower layer at rest,
where Z0 and δz are the mean depth and thickness of this vertical shear layer. The value for the vertical gradient δz = 5 mm was taken from the experimental measurements of vertical density gradient by Stegner et al. (2004) in a similar setup. The three-dimensional mean density field ρ(r, z) is then constructed from the above velocity Vθ(r, z) according to the thermal wind balance,
The static equilibrium (i.e., ∂zρ ≤ 0) of this mean density field was verified everywhere. To trigger out the unstable modes in the simulations, azimuthal periodic density perturbations are added on the basic state along the front, as in Bouruet-Aubertot and Echevin (2002). The modified radius of the density front is then . The wavenumber ki spans the range 2–30 and ε ranges between 0.5% and 20%. More details are given in section 3. The initial velocity field , which is the sum of the mean velocity profile Vθ(r, z) and the azimuthal variations , is then computed from the perturbed density field according to the gradient wind balance relation.
d. Numerical model
The numerical code is the ocean global circulation system NEMO in version 2 (Madec 2008). It solves the rotating hydrostatic primitive equations within the Boussinesq approximation. The equations are discretized on a Cartesian grid. In our idealized configuration, we used z coordinate with vertical depth levels regularly spaced (Δz = h1/15). The vertical grid resolution Δz is small to accurately reproduce the thin vertical density gradient between the two layers. The Cartesian grid is regular and its horizontal resolution is equal to Δx = Δy = 1.25 × 10−3 m, which is on the order of of the deformation radius, allowing the realistic modeling of mesoscale structures. The grid therefore has 323 × 323 × 100 grid points in the flat bottom case. With a bottom shelf topography, the number of vertical levels increases with the total water depth H and may vary from 100 to 160. Because of the Cartesian grid, the lateral boundary of the cylindrical domain is characterized by spatial irregularities at grid size. Nevertheless, these coastal irregularities are located several deformation radii away from the front, and their impacts on the front dynamics are negligible.
We used a rigid lid boundary condition at the free surface to filter the fast barotropic gravity or Kelvin waves. On the lateral boundary and the flat bottom, a no-slip condition is applied. During the run, a convective adjustment scheme is used to keep static stability (Madec et al. 1991). In the vertical, a harmonic Laplacian operator is used for the diffusion of momentum and salinity. In the horizontal, two types of parameterizations are applied for the diffusion. For the low Reynolds simulations in section 3, an explicit Laplacian diffusion operator is used. The diffusivity coefficients (Austausch coefficients) are the same on both vertical and horizontal directions and correspond to the molecular viscosity. For the case of high Reynolds simulations in section 4, a bi-Laplacian operator is used for the horizontal diffusion to filter out the small-scale structures generated by the turbulent cascade or the computational noise.
Two sets of simulations were performed, the low Reynolds simulations RunLR (Table 2) and the high Reynolds simulations RunHR (Table 3). For the low Reynolds simulations, sensitivity experiments differ according to the value of the topographic parameter To or the value of the added perturbations of the initial density front ε. They are presented in section 3, and their detailed parameters are summarized in Table 2. For the high Reynolds simulations, sensitive experiments differ according to the value of the topographic parameter To or the value of the isopycnal slope α). They are presented in section 4, and their detailed parameters are summarized in Table 3.
3. Baroclinic instability at low Reynolds number
The small scales of the experimental setup induce a low horizontal Reynolds number (Re = 50) and a small Ekman number (Ek = 3 × 10−3). In this case, both the dissipation and the amplitude of the initial perturbations play an important role in the evolution of the buoyant coastal current. In a stable configuration, small nonaxisymmetric perturbations decay and the coastal current remains almost circular. However, for low Reynolds numbers, the circular velocity profile also evolves in time because of the viscous dissipation and the kinetic energy decays. For an unstable current, the averaged circular velocity first decays because of the dissipation while small nonaxisymmetric perturbations are slightly growing. During this first stage the circular symmetry is preserved and the maximum surface velocity can decay significantly (Fig. 4). However, after a given time the unstable perturbations reach a finite amplitude and the buoyant coastal current is fully destabilized (Fig. 6), large meanders appear, and coherent eddies are formed. The time needed to reach this full destabilization depends both on the initial amplitude of the perturbations and their unstable growth rates. Hence, to perform relevant comparisons between laboratory experiments and numerical simulations, we need to start the simulations with the same initial conditions (i.e., same mean flow and same amount of initial noise). If the laboratory experiment is always noisy, we need to add initially some radial perturbations in the simulations to trigger the instability.
The matching of the initial conditions and the dissipation rate are thus discussed in section 3a according to the evolution of integral quantities such as the mean azimuthal velocity and the surface kinetic energy (SKE). The stage of the fully nonlinear destabilization is then quantified and analyzed in section 3b for the laboratory experiments and in section 3c for the numerical simulations. The stability analysis is finally performed in section 3d.
a. Viscous dissipation and initial noise amplitude
The evolution of the mean azimuthal velocity Vθ(r, t), spatially averaged along the azimuth in the whole basin, is given in Fig. 4 for a flat bottom configuration, for both the laboratory experiment EXP0 and the numerical simulation RunLR_s0a. The viscous dissipation controls the first stage of evolution (t < 5T–10T) and the surface velocity decays significantly while the width of the mean current increases. In a second stage (10T ≤ t < 20T), a more drastic change occurs because of the unstable growth of finite-amplitude perturbations, leading to an effective diffusion of the averaged azimuthal velocity profile. In EXP0, the center of the mean azimuthal current shifts toward the center of the basin (Fig. 4); whereas it stays around its initial location in RunLR_s0a. The growth of the perturbations occurs earlier in EXP0 than in RunLR_s0a and may lead to these different behaviors. Both the viscous dissipation and the growth of unstable perturbations control the evolution of the mean flow and we can hardly extract the viscous dissipative time scale from the decay of the mean azimuthal velocity.
The viscous decay or baroclinic growth of unstable perturbations strongly impacts the kinetic energy. We use here the surface kinetic energy, which corresponds to the kinetic energy of the horizontal velocities, measured by PIV in the laboratory or computed in the numerical simulations, at the surface level. The evolutions of the SKE are plotted in Fig. 5 for a weak noise amplitude (3% of the SKE) in RunLR_s0b or a moderate noise amplitude (21% of the SKE) in RunLR_s0a. For the numerical simulations, the SKE first decays because of viscous dissipation. In a second stage, the baroclinic instability induces an energy transfer from potential to kinetic energy and the SKE grows until the nonlinear saturation of the instability is reached. We clearly see in Fig. 5 that the time needed to form large-scale meanders or eddies (i.e., nonlinear saturation of the instability) depends on the amplitude of the initial noise. To simulate an early SKE growth as in EXP0, it was necessary to add a significant amount of initial noise (21% of the SKE) in the run RunLR_s0a. For instance, if the initial amplitude of the nonaxisymmetric perturbation is too weak as in RunLR_s0b, the maximum SKE is reached at tSKE ≃ 40T, 15 rotation periods later than in RunLR_s0a (tSKE ≃ 25T) and 20 rotation periods later than in EXP0 (tSKE ≃ 18T). However, a perfect match of the SKE growth between numerical simulations and the laboratory experiment can hardly be achieved. This is probably due to the ageostrophic and three-dimensional nature of the perturbations (Stegner et al. 2004) induced by the geostrophic adjustment of the lock release experiment. Hence, we can hardly quantify the viscous decay in the early stage of the laboratory experiment when the amplitude of initial perturbations is large.
The viscous dissipation can be estimated according to the initial SKE decay only for a low noise level such as in RunLR_s0b. We fit the decrease of the surface kinetic energy of RunLR_s0b with e-folding decay rate, and we obtain a dissipative time scale τυ ≃ 14T while the rotation period is T = 4π/f = 15 s. We can also quantify a viscous dissipation for a later stage, once the nonlinear saturation of the instability is achieved and eddies are fully formed in the basin. In the final stage, the SKE decay of large-scale eddies should be mainly controlled by the viscous dissipation. We use the e-folding decay rate of the surface kinetic energy of RunLR_s0b and EXP0 when t > 30T to estimate a final dissipative time scale. A good agreement is found between the experiment (τυ = 9.6T for EXP0) and the numerical simulation (τυ = 9.4T for RunLR_s0a). Hence, we conclude that the viscous dissipation of this two-layer rotating flow is accurately reproduced by the hydrostatic NEMO model if an explicit Laplacian dissipation operator is used. One can notice that the dissipative time scale τ is not modified when a bottom shelf slope is present in the lower layer. In low Reynolds (Re = 50) laboratory experiments, the molecular viscosity induces a strong dissipation of the kinetic energy on a typical time scale τυ ≃ 9T–14T, which is much smaller than the nonrotating viscous time scale and slightly larger than the Ekman time .
b. Laboratory results
Typical unstable evolutions of buoyant coastal currents, generated by a lock release setup, are shown in Fig. 6. In the flat bottom case, finite-amplitude perturbations destabilize the initial circular current after 10–15 rotation periods (Fig. 6b). The radial perturbation exhibits an azimuthal wavenumber n = 9 corresponding to an unstable wavelength λ = 2πR1/n ≃ 8.7 cm ≃ 7Rd. At the final stage, these unstable meanders lead to the formation of large-scale eddies that invade the whole basin (Fig. 6c). The relative vorticity ξ/f remains moderate (−0.5 < ξ/f < 0.5), even if some amplification occurs during the nonlinear stage of the instability (Figs. 6b,c) because of frontal stretching of the meanders.
When a moderate bottom slope bathymetry is added, while keeping the vertical aspect ratio δ = 0.15 constant, a strong stabilization of the surface currents occurs. The right panels in Fig. 6 show the dynamical evolution of the surface velocity and vorticity fields in EXP1 with a bottom shelf slope s = 0.25 corresponding to a topographic parameter To = −0.68. In this weakly unstable configuration, the buoyant coastal current exhibits some meanders, but coherent vortices are not generated. The dissipation overcomes the unstable growth of radial perturbations and the surface circulation remains almost circular. Both the maximum velocity and vorticity values decay while the current width increases with time (Figs. 6e,f).
c. NEMO model simulations
The evolutions of the horizontal surface velocity and vorticity fields, from simulations, are shown in Fig. 7 for direct comparisons with EXP0 (Fig. 6). In the flat bottom case, finite-amplitude perturbations destabilize the initial circular current after 10 rotation periods (Fig. 7a). The radial perturbation exhibits an azimuthal wavenumber n = 9–10 corresponding to an unstable wavelength λ = 2πR1/n ≃ 7.8 − 8.7 cm ≃ 6 − 7Rd. Then, the unstable meanders lead to the formation of large-scale eddies that invade the whole basin (Fig. 7c). This numerical simulation RunLR_s0a is in correct agreement with the laboratory experiment EXP0. As is shown in Fig. 5, the numerical modeling is very sensitive to the initial noise perturbation. For a smaller noise amplitude, as in RunLR_s0b, the eddy formation (i.e., nonlinear saturation) occurs much later after 20–30 rotation periods. Hence, a fine tuning of the initial perturbation is needed to obtain a qualitative agreement between the simulations and the laboratory experiments. Indeed, if we run a numerical simulation with a bottom slope s = 0.25 identical to the experiment EXP1, we should reduce the initial noise amplitude to 3% to obtain a qualitative agreement with the laboratory experiment. The right panels of Fig. 7 correspond to the simulation RunLR_s25 with the same topographic parameter as EXP1. As in the laboratory experiment, the surface flow is weakly unstable and an azimuthal wavenumber n = 10 is visible on the vorticity field (Fig. 7e). The meander amplitude is reduced and coherent vortices do not emerge. This is probably the signature of a reduced growth rate induced by the shelf slope. Indeed, if the growth rates and/or the initial noise amplitudes are too small, the dissipation overcomes the growth of radial perturbations and prevents the formation of coherent eddies.
d. Stability analysis from the NEMO model
Unlike the laboratory experiment, the linear stage of the simulations can provide a first estimate of the unstable growth rate. The coarse resolution (64 × 64 grid points) and the weak sensitivity of the PIV prevent an accurate measurement of small-amplitude perturbations in the experiments. However, the model resolution (323 × 323) allows a spectral decomposition of the azimuthal modes, and if the initial noise amplitude is weak enough the linear growth of the unstable perturbation can be quantified. Hence, in the following we look at numerical runs with small initial noise (amplitude less than 5%). During the linear stage, the radial velocity is interpolated, each inertial period Tf = T/2, along a circle of radius R1 located at the maximum velocity of the coastal current. A fast Fourier transform is then performed on this signal to calculate the energy of each azimuthal mode. The temporal evolution of the amplitude of each mode is fitted with an exponential law to estimate the corresponding growth rate σ. The dimensionless unstable dispersion relations are plotted in Fig. 8 for three runs, RunLR_s0b, RunLR_s25, and RunHR_s00, differing by their Reynolds number or topographic parameter. The linear growth rates are rescaled by the initial azimuthal velocity U and the deformation radius Rd, whereas the wavenumber k is rescaled by Rd. For the flat bottom case (RunLR_s0b), we find a maximum growth rate of the unstable perturbations for kRd ≃ 0.8 in good agreement with the simplified Phillips model of baroclinic instability, which predicts the maximum growth for kRd = 0.65 when δ = 0.15 (see appendix). The wavelength selection of the unstable mode seems to be correctly approximated by standard baroclinic instability; nevertheless, the growth rates are strongly overestimated. Indeed, the Phillips model predicts a maximum value around σmaxRd/U = 0.15, whereas the perturbations in the simulations exhibit a much smaller growth rate σmaxRd/U ≃ 0.03. When a bottom shelf is added with a slope s = 0.25 (RunLR_s25a), a shift in the wavelength selection occurs and the growth rates are reduced. The maximum growth rate σmaxRd/U ≃ 0.017 corresponds here to an e-folding time τ ≃ 13T, which is very close to the viscous decay e-folding time τυ = 9T–14T estimated above. Hence, the viscous dissipation strongly reduces the growth of unstable perturbation, especially when the bottom shelf slope tends to stabilize the buoyant coastal current.
4. Baroclinic instability at high Reynolds number
In the ocean, unlike laboratory experiments, the molecular viscosity is neglected and the momentum diffusion is mainly controlled by small-scale turbulent advection. To reproduce the high Reynolds dynamics of real coastal flows and allow the comparison with the inviscid quasigeostrophic theory, we performed several numerical runs (Table 3) using a bi-Laplacian operator for horizontal motions and a standard diffusion on vertical motion corresponding to low Ekman numbers (Ek = 3 × 10−6). In the case of low dissipation such as in RunHR_S00a, the unstable growth rate is increased compared to low Reynolds simulations (Fig. 8). In section 4a, we analyze the geostrophic nature of the instability in high Reynolds simulations. Then, in section 4b, we quantify the impact of the bathymetry on the linear stage development. Finally, section 4c describes the nonlinear saturation regime and the eddy formation.
a. The geostrophic nature of the instability
Both geostrophic or ageostrophic instabilities may destabilize a buoyant coastal current. According to the large vorticity values (−0.5 ≤ ξ/f ≤ 0.5) of the initial outcropping current (Figs. 6, 7, 14) one may suspect some unstable coupling between the geostrophic Rossby modes and the ageostrophic Kelvin or Poincaré wave modes (Sakai 1989; Gervasio 1997; Gula and Zeitlin 2010a,b). Hence, in what follows, we try to characterize the dynamical nature of the buoyant coastal current instability, for the flat bottom and for the steep shelf slope configurations.
According to the Charney–Stern criterion, opposite PV gradients in the upper and lower layers are necessary to allow an unstable coupling between Rossby modes (geostrophic baroclinic instability). The PV profiles Qi of the initial axisymmetric current in the top (i = t) and the bottom (i = b) layers are shown in Fig. 9. Even if the simulations are performed with a continuously stratified model, we choose here the shallow-water formulation Qi=(ξi + f)/hi for the PV due to the specific two-layer stratification we used [Eqs. (2) and (3)]. Indeed, the sharp density gradient between the surface and the bottom water induces a virtual interface between the upper and the lower layers and allows a clear definition of their respective thicknesses hi. We first note that the PV is monotonic in each layer, and therefore the flow is expected to be stable in case of barotropic shear perturbations. However, for all the simulated cases, the horizontal PV gradient between the coast and the open sea is positive in the upper layer and negative in the lower layer. Indeed, the surface PV increases from an almost constant value Qt ≃ Q0 = f/h1 near the coast where ht(x) ≃ h1 and diverges close to the outcropping front where ht(x) → 0, whereas in the bottom layer the PV decreases from the coast Qb ≃ f/h2 to the open sea Qb ≃ f/H because of the increase of the bottom layer depth H ≥ h1 + h2. The steep shelf slope may strongly amplify the PV gradient in the bottom layer but does not change its sign, and the surface current remains potentially unstable according to the Charney–Stern criterion.
The baroclinic instability is characterized by the ability of the flow to convert the APE into kinetic energy. Because of the large width L = 5Rd–6Rd of the buoyant water (i.e., small Burger number), the initial coastal flow configuration corresponds to a large amount of APE. According to Fig. 10, the release of the initial APE induces an increase of the total KE and corresponds to the amplification of the unstable perturbations. Indeed, for the flat bottom simulation (RunHR_s00a), the KE increases at coincides with the current meanders and eddy formation (Fig. 14a). Hence, the APE provides the energy for the growth of unstable modes within the coastal current, as is the case for the standard baroclinic instability.
The wavelength selection generally differs between geostrophic and ageostrophic instabilities; therefore, the analysis of the most unstable wavenumber gives information about the instability. According to the dispersion relation in Fig. 8, for both the high and the low Reynolds regimes, the highest growth rate occurs for kRd ≃ 0.8 (i.e., λ = 2πRd/0.8 ≃ 7.8Rd). This value is close to the prediction of the standard Phillips model describing the unstable coupling between two Rossby modes. This unstable wavelength, for the flat bottom case, is much larger than the wavelength predicted by the ageostrophic coupling between a Rossby and a Kelvin wave (λ ≃ Rd), for instance. Hence, the typical sizes of the unstable meanders do not correspond here to the wavelength selection induced by ageostrophic instabilities.
We then perform a careful analysis of the spatial structure of the most unstable mode in both layers. We first decompose each component a(r, θ, z) of the flow (a stands for the velocity υ or the density ρ) into a mean axisymmetric part and an azimuthal perturbation . In the linear stage of the instability, when the nonlinear coupling between modes can be neglected, the perturbations associated with each azimuthal wavenumber correspond to the unstable eigenmodes. From the perturbed velocities , we compute the vorticity of the velocity perturbations . During the linear stage of the instability, the vorticity perturbations (with a standing for ξ in the definition above) is equal at a first order of approximation to the vorticity of the velocity perturbations . The perturbation vorticity fields computed at a very initial stage (t = 10T) is shown in Fig. 11 for a flat bottom and for a steep shelf slope configuration. Figure 11a shows the azimuthal perturbation associated with the eigennumber n = 9 (the most unstable mode) for the flat bottom case. For both the upper and the lower layers, the unstable perturbations are localized in the region of strong PV gradient, which corresponds here to the core of the coastal current. There is no signature of unstable perturbations close to the coast such as Kelvin wave modes. Besides, the perturbed velocity and density fields satisfy the thermal wind balance. Hence, these azimuthal perturbations are geostrophically balanced in both layers. Figure 11b shows the azimuthal perturbation associated with the eigennumber n = 12 for a steep shelf slope configuration s = 50% and To = −1.3. In this case, the unstable perturbations in the bottom layer extent on a wider area along the shelf slope. This spatial structure is similar to a topographic Rossby wave pattern. Here again, the azimuthal perturbations are geostrophically balanced in both layers.
Hence, according to the energy budget, the wavelength selection, and the spatial structure of the unstable modes, we can conclude that this coastal current instability corresponds to a standard baroclinic instability: that is, to the coupling of geostrophic Rossby modes between the surface and the deep lower layers.
b. Topographic impact on the linear unstable growth
To quantify the impact of the bathymetry, we use the topographic parameter To = s/α, where s is the shelf slope and α is the isopycnal slope of the buoyant coastal current. For such a current, shelf and isopycnal slopes are in opposite directions and the parameter To is therefore negative. In the following, we vary To while keeping the other parameters (δ, Bu, Ro, and Ek) constant. The impacts of the relative shelf slope on the most unstable growth rate σm and the corresponding wavenumber km are shown in Figs. 12 and 13. The maximum growth rate decays when To becomes negative, and below a critical value Toc ≃ −2.7 the coastal current is stable. Hence, a clear stabilization of the baroclinic instability occurs when the relative shelf slope increases. This stabilization is due to the strong increase in phase speed of the bottom Rossby mode. Indeed, the shelf slope induces a topographic Rossby mode in the lower layer (Fig. 11b). When the phase speed in the lower layer is too large, the unstable coupling between the upper and the lower Rossby modes cannot occur. A similar stabilization is found with the Phillips model when a linear bottom slope is added (cf. the appendix). However, the critical value needed to stabilize the baroclinic instability is much larger Toc ≃ −23 in this simplified quasigeostrophic model. As far as the wavelength selection is concerned, the increase of the shelf slope shifts the unstable mode to smaller wavelengths (i.e., larger wavenumbers). According to Fig. 13, in a steep shelf configuration (To = −2.5), close to stabilization, the dimensionless wavenumber reaches a value of kRd = 1.3, whereas it was only kRd ≃ 0.8 in the flat bottom configuration. The baroclinic Phillips model (solid line in Fig. 13) exhibits the same trend in the wavelength selection.
To check the relevancy of the topographic parameter To, we perform several runs where the shelf and the isopycnal slopes are changed while To is kept constant. To reduce the slopes, we increased the horizontal scales without changing the vertical ones. The Coriolis parameter f is changed accordingly to keep the Burger (Bu) and the Rossby (Ro) numbers constant. According to Fig. 13, the shelf slope could be varied from s = 0.65 to s = 0.12 without noticeable changes in the unstable growth rate or the most unstable wavelength if To remains unchanged. Hence, the shelf slope does not impact directly on the linear stability of the coastal current. We confirmed here that the relevant parameter that controls the stability and the wavelength selection of a buoyant coastal current over a steep shelf is the topographic parameter To as was suggested by quasigeostrophic theory, in the framework of the two-layer Phillips model (Mysak 1977) or generalized Eady models (Blumsack and Gierasch 1972; Mechoso 1980; Isachsen 2011). We also checked that the nonlinear evolution of the coastal current and the eddy generation are only controlled by the relative shelf slope parameter To = s/α (and not the absolute slope values) for the hydrostatic NEMO runs.
c. Nonlinear saturation and eddy generation
The nonlinear saturation, leading to the formation of meanders and mesoscale or submesoscale eddies, is a key process of the cross-shelf transport. The shelf slope may have a strong impact on the trajectories or the robustness of these eddies (Sutyrin et al. 2003, 2009). The nonlinear evolution of the instability and its impact on the intrusions of dense water mass in the coastal zone are shown in Fig. 14, where both a flat bottom (To = 0) and a steep shelf slope (To = −1.3) configurations are presented. The formation of finite-amplitude meanders is shown in Figs. 14a,d. The typical scale of these nonlinear meanders is controlled by the linear wavelength selection and decreases when the shelf slope gets steeper as shown in Fig. 13. Indeed, 9 meanders are formed in the flat bottom case (To = 0), compared to 12 when To = −1.3. In the deep bottom layer (not shown), dipolar structures are formed just below the surface meanders. The coupling induced by the linear baroclinic instability is still active during the nonlinear stage, and the bottom layer dipoles induce a radial stretching of the surface meanders. Hence, the density front is shifted toward the coast and leads to the formation of mesoscale cyclonic eddies (Figs. 14b,e). Once the eddies are formed and detached from the initial density front, they are submitted to secondary nonlinear processes that affect their size and shape. For the flat bottom case (Fig. 14c), the cyclonic mesoscale eddies merge together, following an inverse cascade, and lead to larger mesoscale vortices. On the other hand, for a steep shelf slope configuration, mesoscale eddies tend to be stretched and split into smaller submesoscale cyclones (Fig. 14f).
To quantify more precisely the impact of the shelf slope on the formation of coastal eddies, we quantify the intensity, the horizontal size, and the thickness of the cyclonic structures that are formed in the surface layer. The evolution of the maximum vorticity is shown in Fig. 15. The initial value corresponds to the intensity of the cyclonic shear at the edge of the axisymmetric buoyant coastal current. This initial value is relatively high ξ/f ≃ 0.5 because of the outcropping configuration. Then, because of the weak dissipation in the high Reynolds simulations, a moderate decay of the front vorticity is induced until the unstable perturbation grows sufficiently to form large meanders or eddies at for the flat bottom case or for a very steep configuration (To = −2.6) close to the stability threshold. At that stage, the baroclinic instability generates finite-amplitude perturbations and the potential to kinetic energy transfer induces a strong amplification of the vorticity. At the final stage, when cyclonic eddies are detached from the coastal front, the relative vorticity at the edge or in the core of the cyclones may reach values up to ξ/f ≃ 1.3, which are much larger than the initial values. According to Fig. 15, this vorticity amplification is not affected by the amplitude of the shelf bathymetry, and the cyclonic eddies reach the same intensity over a flat bottom or a steep shelf slope.
Various dynamical criteria could be used to quantify the location and the size of a coherent vortex (Pasquero et al. 2001; Isern-Fontanet et al. 2004). However, as far as cross-shelf transport is concerned, we define a vortex as a coherent structure able to trap and isolate a water parcel in its core. Hence, to identify the water parcels, we plot the density field at a given depth (z = −h1/2), and, according to the threshold density value (ρc = (ρ1 + ρ2)/2), we can separate the light coastal water (ρ ≤ ρc) from the dense water (ρ ≥ pc) coming from the central basin (i.e., the open sea).
The initial front between dense and light water is clearly visible in Figs. 16a,d, and the coastal cyclonic eddies corresponding to the inflow of dense water parcels along the coast can be detected in Figs. 16b,e. Once they are formed, we can estimate an averaged radius of cyclonic eddies from the area A of dense fluid parcels, or in other words the surface area of the trapped region in the cyclonic core. Figure 17a shows the impact of the shelf bathymetry on the sizes of the cyclones generated by the buoyant coastal current instability. The cyclonic eddy radii rci are estimated just after the nonlinear saturation when the first eddies are detached from the coastal front. The mean cyclonic radius rc is computed by averaging the radius with the number of detected cyclones Nc. The error bars correspond to the standard deviation of the identified eddies, and we stop the computation of mean values when the number of detected eddy is too small (i.e., Nc ≤ 4). Figure 17a shows that the typical radius rc/Rd of the first cyclones generated by the baroclinic instability decays when the topographic parameter reaches negative values. These cyclonic radii follow the evolution of the most unstable wavelength λ according to the linear stability analysis (Fig. 13) and correctly match the relation rc = ¼λ. Hence, the eddy sizes are initially controlled by the linear instability process. However, at later times, secondary nonlinear processes may strongly affect the eddy size. Different backgrounds are used in Fig. 18 to distinguish between three dynamical stages. The white region corresponds to the exponential growth of infinitesimal perturbations. During this linear stage, the front meanders are small and eddies are not formed yet. Then, the nonlinear saturation (light gray) occurs when the meanders reach a finite amplitude and lead to the formation of isolated eddies (Figs. 16a,d). Finally, in the gray regions, the eddies are fully developed and secondary nonlinear processes occur. For the flat bottom case, the secondary merging process is characterized by an increase of the mean radius while the total number of eddies decrease. This inverse cascade leads to mesoscale cyclones with a mean radius rc much larger than ¼λ, the radius predicted by the linear wavelength selection (Fig. 18). However, for steep shelf configurations (e.g., for To = −2.6), once the vortices are fully developed and detached from the density front (gray region) some stretching and splitting induced by the bottom topography lead to much smaller eddy sizes (Figs. 16e,f). The mean cyclonic eddy radius rc decreases toward submesoscale values such as rc/Rd ≃ 0.6 − 0.7. These values are much smaller than the value predicted by the linear stability analysis ¼λ ≃ Rd (Fig. 18).
The depth of cyclonic eddies (i.e., dense water intrusions) can be estimated from the three-dimensional density field. Once a patch of dense fluid (ρ ≥ ρc) with closed contours is detected on the surface density field, we can calculate the maximum depth dc of this three-dimensional lens. The mean relative depth dc/h1 of cyclonic eddies are plotted in Fig. 17b as a function of the topographic parameter To. Unlike, the mean radius (Fig. 17a), the typical thickness of these mesoscale cyclones remains almost constant and seems to be weakly affected by the bottom bathymetry.
Considering the cyclonic intrusions of dense water in coastal areas induced by the baroclinic instability of buoyant coastal currents, the main impact of a steeper shelf bathymetry is to induce smaller eddies. Although the intensity or the vertical extent of these cyclonic lenses is weakly affected by the bottom slope, the mean areas of these dense water parcels are nevertheless strongly controlled by the bathymetry. The typical cyclonic eddy radius first decreases because of the linear instability process, and afterward nonlinear processes amplify this tendency and lead to the formation of smaller submesoscale vortices over the shelf slope.
5. Discussion and conclusions
The stability of buoyant coastal currents above a steep shelf slope was investigated by both laboratory experiments and numerical simulations. Unlike previous papers (Lozier and Reed 2005; Rivas et al. 2005; Wolfe and Cenedese 2006), where the shelf slope s was used to quantify the bathymetric effect, we use here the topographic parameter To = s/α, the ratio of the bottom shelf slope s over the surface isopycnal slope α. We follow here the simplified quasigeostrophic studies (Blumsack and Gierasch 1972; Mysak 1977; Mechoso 1980; Isachsen 2011) that demonstrate that To is a central parameter that controls the impact of the bottom shelf slope on the surface current stability. Moreover, in the framework of two-layer stratification, we separate the influence of the topographic parameter To from the vertical aspect ratio parameter δ (the ratio of the upper-layer thickness over the total water depth below the front), which controls the baroclinic instability over a flat bottom (Pedlosky 1987; Vallis 2006). When both parameters are varied together, the impact of the shelf slope on coastal fronts seems unclear with contradictory results. Hence, to clarify the situation, we mainly vary the topographic parameter To while keeping the vertical aspect ratio δ and other dynamical parameters constant.
The hydrostatic NEMO model was first used with a standard dissipative operator (Laplacian) and a molecular viscosity coefficient to perform quantitative laboratory–numerical comparisons. If the initial noise amplitude of the numerical simulation is accurately fixed, the unstable evolution of the buoyant coastal current, measured by PIV in the rotating tank experiment, is correctly reproduced by the idealized NEMO model. The unstable growth of the surface current meanders and their characteristic wavelengths are significantly affected by the bottom shelf slope. The laboratory and the numerical results both show a strong stabilization of the buoyant coastal current and a smaller wavelength selection when the shelf slope becomes larger than the isopycnal slope of the surface front. We note that the impact of the topographic parameter To on the wavelength selection is correctly predicted by an oversimplified two-layer quasigeostrophic model. Hence, for a geostrophic (small Rossby number Ro ≃ 0.2–0.3) and surface-advected (small δ = 0.15) coastal front, the standard baroclinic instability appears to be the leading instability. Nevertheless, we show that the viscous dissipation of a small laboratory setup strongly reduces the growth of unstable perturbations, especially when the bottom shelf slope tends to stabilize the buoyant coastal current.
To reproduce the high Reynolds dynamics of real coastal flows, we performed several numerical runs using a bi-Laplacian operator for horizontal motions and a standard diffusion for vertical motions corresponding to a very low Ekman number (Ek = 3 × 10−6). The energy budget corresponds to a standard baroclinic instability where the release of the available potential energy induces an increase of the total kinetic energy and drives the amplification of unstable perturbations. The linear stage of the instability was also quantified from the numerical runs. The linear growth of the perturbation spectrum leads to the formation of mesoscale meanders λ ≃ 2πRd in agreement with the wavelength selected by the unstable coupling between two Rossby modes (Phillips 1954). The unstable coupling with ageostrophic Kelvin or gravity modes would have selected a much smaller wavelength here (Sakai 1989; Gervasio 1997; Gula and Zeitlin 2010a,b). Besides, we have shown that the unstable modes in both the upper and the lower layers satisfy the geostrophic balance. Hence, according to the energy budget, the wavelength selection, and the spatial structure of the unstable modes, we can conclude that the coastal front instability studied is driven by the standard baroclinic instability.
An important result of this study is the confirmation, with a fully nonlinear primitive equation model and laboratory experiments, that the topographic parameter To is the relevant parameter to quantify the impact of a steep bottom slope on the stability of buoyant coastal current. We show that a complete stabilization of the coastal front can be reached for finite negative values: for instance, Toc ≃ −2.7 when δ = 0.15, Ro = 0.3, and Bu = 0.02. These results are in agreement with the recent studies of Spall (2004) and Isachsen (2011), who show a strong decrease of the eddy heat flux (i.e., stabilization) of a thermally forced coastal current when the topographic parameter reaches similar values.
Another important result of this study is the evidence that a steep shelf slope affects the nonlinear development of the baroclinic instability of a buoyant coastal current and may lead to the formation of submesoscale eddies. In a flat bottom configuration, the baroclinic instability generally leads to mesoscale eddies (λ ≃ Rd and larger) and only ageostrophic instabilities are known to generate smaller eddies. In this study, we exhibit a new dynamical sequence leading to the formation of submesoscale structures, in the context of weakly unstable geostrophic modes. At the linear stage of the instability, a steeper shelf slope (increase in To) tends to stabilize the coastal buoyant current. For this stabilization process, smaller wavelengths are selected. For intermediate cases, when the coastal front is weakly unstable, the nonlinear saturation leads to smaller eddies. We have seen that, once cyclonic intrusions of cold seawater are formed above the coastal shelf, a secondary process tends to stretch these coherent cyclones and to split them into smaller submesoscale eddies. Then, these small cyclones may reach a characteristic radius rc ≃ 0.6Rd − 0.7Rd over a steep shelf (To ≃ −2.6) while for a flat bottom configuration (To = 0) only large mesoscale eddies can be formed rc ≃ 2Rd − 2.5Rd.
This direct cascade to small-scale structures could be compared to the spectral energy fluxes of a thermally forced coastal current calculated by Isachsen (2011). Following the procedure used by Scott and Wang (2005), the authors found that there is a negative flux (inverse kinetic energy cascade) at scales roughly larger than the deformation radius Rd, whereas a weaker positive flux occurs at scales smaller than Rd. In our case, without any external forcing, a finite amount of kinetic energy is initially released by the baroclinic instability at a given scale. For the flat bottom configuration, this scale is larger than the deformation radius (Fig. 18a) and the nonlinear eddy–eddy interaction leads to an inverse energy cascade in agreement with calculations of Isachsen (2011). However, for the steep shelf slope configuration, the initial release of kinetic energy occurs at a smaller scale (almost Rd). This initial spectral distribution of the kinetic energy may constrain the nonlinear evolution of the flow and emphasis on the direct energy cascade that occurs at submesoscale. However, the exact nature of this direct energy cascade that occurs at submesoscale and the role of steep bottom slope on the splitting process are not explained yet and should be studied in the future.
The two anonymous referees are gratefully acknowledged for their valuable remarks and comments. We thank Samuel Teinturier for helping in the initiation of this study. This work was granted access to the HPC resources of the Institut du Dévelopement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS) under allocation 2009 and 2010 (project 010227) made by Grand Equipement National de Calcul Intensif (GENCI). This study was mainly supported by the TopIECC project funded by the French program LEFE-INSU of the CNRS (http://www.insu.cnrs.fr/) and the ANR-Flowing project funded by the ANR “Programme Blanc” (http://www.agence-nationale-recherche.fr/). R. Pennel acknowledges the financial support of the grant from Direction Générale de l’Armement (DGA).
Phillips Instability Problem over a Sloping Bottom
This simplest baroclinic instability problem considers a two-layer system where the velocities in each layer Ui (i stands for the layer 1 or 2) are independent of y but differ in magnitude. We note hi, the depth of each layer, and α, the mean slope of the interface between the two layers. A sloping bottom of height ηb is added in the lower layer, and we note s = ∂yηb, the mean slope of this bathymetry.
We scale the basic variables according to
and we introduce the nondimensional parameter,
The dimensionless quasigeostrophic potential vorticity qi in each layer can be written as
with ψi the geostrophic streamfunction (∂yψi = −Ui).
The potential vorticity in each layer is advected by the geostrophic velocity [∂tqi + J(ψi, qi) = 0], and therefore the dimensionless two-layer quasigeostrophic model is written as
where J(a, b) = ∂xa∂yb − ∂ya∂xb is the Jacobian operator.
To study the linear stability of that flow, we decompose the streamfunction ψ as follows: ψ = Ψ(y) + φ(x, y, t), where φ is a small perturbation, and we use U1 = Us = U2 − U1 and U2 = 0. The linearized equations are then
The perturbation φ may be decomposed into normal modes,
and we define K2 = k2 + l2
For nontrivial solution, the determinant of coefficients must be zero. This gives a quadratic equation in c,
and solving this we obtain
Finally, we exhibit the growth rate, σ = k × Im(c). We find, for instance, that the maximum growth rate and the most unstable wavelength for the flat bottom case (To = 0) are σ = 0.136 and k = 0.6. This model gives a complete stabilization (σ = 0) of the instability for To = −23.