## Abstract

It is well known that upwelling and downwelling currents are unstable to perturbations. Less is, however, known about the physical mechanism responsible for the observed and modeled instabilities. It is shown that the origin of the long-wave barotropic/baroclinic instability observed on upwelling currents has to be sought among diabatic or thermobaric mechanisms. In particular, the role of mixing associated with Kelvin–Helmholtz instability and of wind forcing is investigated. Low Richardson numbers occur when the pycnocline outcrops at the sea surface. The criterion for instability (Ri ≤ 1/4) can be reached in a narrow region close to the upwelling front, permitting Kelvin–Helmholtz instability and mixing. This can precondition the current for long-wave instability by transforming the current's potential vorticity. A constant wind can likewise modify the potential vorticity. The resulting potential vorticity anomaly is always negative for both upwelling and downwelling currents, and this anomaly interacts with the outcropped front, destabilizing it. Examples are provided via numerical calculations using an idealized front. A wind stress is an effective means of inducing the negative PV necessary for instability; with wind, Kelvin–Helmholz instability, when present, merely modifies the instability characteristics. In addition, upwelling fronts are always less stable than comparable downwelling fronts.

## 1. Introduction

### a. Foreword: Potential vorticity properties

In this paper potential vorticity (PV) is the key quantity, and we analyze its creation and evolution for upwelling currents. It is thus useful to recall the basic properties of PV for geophysical flows:

In an incompressible fluid the PV is conserved on fluid parcels, in the absence of diabatic processes and thermobaric effects (Ertel 1942; Marshall et al. 2001).

A geostrophic velocity field can be inferred from a given PV field (see, e.g., McWilliams and Gent 1980; Hoskins et al. 1985; McIntyre and Norton 1990).

Because of these fundamental principles, many geophysical fluid processes can be interpreted by analyzing the PV evolution (see McWilliams and Gent 1980; Hoskins et al. 1985; McIntyre and Norton 1990; Birkett and Thorpe 1997). The PV distribution is also important for flow stability. Charney and Stern (1962, hereinafter CS) famously showed that a flow can be baroclinically or barotropically unstable only if the PV gradient changes sign somewhere in the current (Ripa 1991; Cushman-Roisin 1994; Morel and McWilliams 2001).

An additional principle, of relevance to upwelling, is Bretherton's notion of a generalized PV. The intersection of isopycnals at the sea surface implies that potential density varies along the surface. Bretherton (1966) showed that “a flow with potential [density] variations over a horizontal and rigid plane boundary [as is the case of the outcropping of upwelling systems at the sea surface] may be considered equivalent to a flow without such variations, but with a concentration of PV very close to the boundary.” The consequences of this were studied by Boss et al. (1996) in the context of an upwelled front. They showed that the region in which the layer vanishes behaves like a region of very high PV, with the intersection point acting like a strong PV gradient; this in turn affects the front's stability.

### b. Subject and plan of this study

The Coastal Ocean Dynamics Experiment (CODE; see Beardsley and Lentz 1987 and references therein) and the Coastal Transition Zone program (CTZ; see Brink and Cowles 1991 and references therein) described the dynamics of upwelling currents along the California shelf. Since then, observations of upwelling systems have been performed in many other regions [see Haynes et al. (1993) and Relvas de Almeida (1999) for coastal upwellings and Buckley et al. (1979), Johannessen et al. (1983), Johannessen et al. (1987), Sandven and Johannessen (1993), Gammelsrod et al. (1975), Clarke (1978), Roed and O'Brien (1983), and Lygre et al. (2002) for upwelling developing along the marginal ice zone]. All of these observations suggest that upwelling currents exhibit long-wave instabilities. As a result, the fronts continually generate vortices and filaments. Observations also often show short-wave instabilities.

The instability of upwelling fronts has been studied both numerically and analytically (Paldor 1983; Barth 1989a, b; McCreary et al. 1991; Barth 1994; Shi and Roed 1999; Roed and Shi 1999; Bouruet-Aubertot and Echevin 2002) and in rotating tank experiments (Narimousa and Maxworthy 1985, 1987; Bouruet-Aubertot and Linden 2002). These studies concur that the long-wave instability is mixed barotropic/baroclinic in character. There is no consensus however on the short-wave instability, which evidently is physically distinct from the long-wave instability. Boss et al. (1996) suggest it is the result of an interaction between a vortical wave and a gravity wave (see also Sakai 1989) while others noted its dependence on vertical shear (Barth 1989a, b) and connection with Kelvin–Helmholtz instability (Paldor and Ghil 1991; Barth 1994).

In this paper we focus on the long-wave instability. As noted, this barotropic/baroclinic instability occurs when the PV gradient changes sign, as was the case in all of the aforementioned studies on frontal instability. Many configurations with two layers indeed assumed no flow in the lower layer, which ensured a PV gradient in this layer, the opposite sign PV gradient being associated with the upper-layer current or front (see, e.g., Boss et al. 1996). However, the dynamical origin of the PV gradients, and more generally the mechanisms establishing the PV structure of the current, remain obscure.

Imagine starting with an ocean at rest (Fig. 1a); the PV field is constant in each layer and there are no unstable PV gradients. A wind blowing parallel to the coast (see Fig. 1b) will generate an upwelling or downwelling current. Without diabatic processes, PV is conserved and no gradient is formed along isopycnal surfaces. However, once the front outcrops at the surface, a high positive PV gradient occurs there, following Bretherton's principle. But we still require a negative PV gradient, to verify the CS criterion and set off barotropic/baroclinic instabilities.

Of course, the initial state could have an existing negative PV gradient.^{1} Barring that, however, one must invoke diabatic processes to initiate instability. Several such processes are

thermobaric effects;

diapycnal mixing associated with Kelvin–Helmholtz instability, internal gravity waves breaking, or convection in mixed layers;

isopycnal mixing and viscous boundary layer effects;

effect of the wind stress; and

mixing associated with atmospheric forcing (including wind stirring).

In this paper, we focus on two: diapycnal mixing associated with Kelvin–Helmholtz (KH) instability and wind stress. In the next section, we define the general framework and equations used in this paper. In section 3, we examine analytically the possibility of Kelvin–Helmholtz instability for upwelling or downwelling currents. The effect of the wind stress is examined in section 4 and numerical results are presented and analyzed in section 5. Our results are summarized and discussed in the final section.

## 2. General equations and model

### a. Equations

The equations that we consider are the primitive equations written in isopycnal coordinates. The momentum and continuity equations along isopycnal surfaces are (Bleck and Boudra 1986; Bleck and Smith 1990; Bleck et al. 1992)

Here **U** = (*U*, *V*) is the horizontal velocity field,^{2 }*f* is the Coriolis frequency (which we assume to be spatially constant in this paper), *h* = −∂* _{ρ}z* is the differential thickness of an isopycnal layer (Cushman-Roisin 1994),

**T**

^{w}= (

*T*,

^{w}_{x}*T*) represents the wind forcing,

^{w}_{y}**F**= (

*F*,

_{x}*F*) and

_{y}*S*are nonconservative terms representing mixing or diffusion (of numerical origin in this paper), and

*M*= (

*P*+

*ρgz*)/

*ρ*is the Montgomery potential (pressure along an isopycnal surface), which can be related to

_{o}*h:*

where ∂_{ρρ} represents the second derivative with respect to *ρ*. When Eqs. (1)–(2) are discretized vertically, *h* can be interpreted as the layer thickness and the Montgomery potential in each layer *k* is a linear combination of all the layer thicknesses [*i* = (1, . . . , *N*), where *N* is the number of layers] and can be written

In the following, Eqs. (1) and (3) will be solved numerically using a version of the Miami Isopycnal Coordinate Ocean Model (MICOM; see Bleck and Boudra 1986; Bleck and Smith 1990; Bleck et al. 1992). The MICOM model has been modified to include a fourth-order scheme for the nonlinear advective terms in the momentum equations and a biharmonic diffusion operator and **F** is then written (see Herbette et al. 2003) as

where *ν* = max(*C*^{1}_{s}Δ*x*^{3}|**U**|, *C*^{2}_{s}Δ*x*^{2}*σ*), Δ*x* represents the grid spacing, and *σ* = [(∂_{x}*U* − ∂_{y}*V*)^{2} + (∂_{x}*V* + ∂_{y}*U*)^{2}]^{1/2} is the deformation tensor. Coefficients *C*^{1}_{s} = 1/32 and *C*^{2}_{s} = 1 were chosen to overcome numerical noise associated with dispersive effects and the development of shocks in the vicinity of the upwelling front.

### b. Potential vorticity

The PV is defined by

where *ζ* = rot(**U**) = ∂_{x}*V* − ∂_{y}*U* is the curl of the velocity field along isopycnal surfaces (calculated with a constant density *ρ*). As in Morel and McWilliams (2001), we define the PV anomaly (PVA) by

where *H* = *h*_{rest} is the stratification of the fluid at rest, and Δ*h* = *h* − *H*. On the *f* plane, Δ*Q* is also a tracer and is conserved for each fluid particle, as PV, but it has the dimensions of vorticity, which makes it easier to analyze.

### c. Richardson number

For Kelvin–Helmholtz instability to develop, the Miles and Howard criterion states that the Richardson number has to be less than a critical number: Ri ≤ 1/4 (see Miles 1961; Howard 1961; and appendix A). It thus seems important to study whether the latter criterion can be reached when upwelling occurs.

The Richardson number is a nondimensional number comparing the local stratification (Brunt–Väisälä frequency) with the vertical gradient of the velocity field. In isopycnal coordinates, it can be written (see appendix A)

## 3. Kelvin–Helmholtz instability

If the primary source of PV anomaly generation is associated with Kelvin–Helmholtz instability, then we can consider that the upwelling development is adiabatic and that PV is constant along isopycnals as long as Ri has not reached the critical value of 1/4. In this case, according to Eq. (6), the vorticity associated with the current (horizontal gradient of the velocity field) is related to *h*:

In addition, upwelling currents are very close to geostrophic equilibrium. The gradient wind balance then relates the vertical derivative of the current to the horizontal gradient of density:

### a. Analytical results for simplified configurations

To evaluate the Richardson number of an upwelling current, let us calculate the Richardson number of a general configuration for a coastal flow in geostrophic equilibrium and with no PV anomaly. We also first hypothesize that there is no outcropping of isopycnals at the sea surface or at the bottom, which yields

where *ρ*_{surf} and *ρ*_{bot} are the surface and bottom densities, respectively (which are constant when there is no outcropping of isopycnals), and *H _{o}* is the ocean depth (which is assumed constant, too).

where *M** = *M* − *M*(*ρ*) is the Montgomery potential anomaly and ∂_{ρρ}*M* = −*g*/*ρ _{o}H*(

*ρ*).

Equations (11) and (12) can be solved using a splitting technique. The general solution can then be written

where Γ* _{n}* is a vertical eigenmode (and −

*γ*

^{2}

_{n}is its corresponding—negative—eigenvalue) of the operator 1/

*H*(

*ρ*)∂

_{ρρ}() and satisfies

and

where *δ _{n,p}* is the Kronecker symbol (=1 if

*n*=

*p*and =0 otherwise).

This yields for the velocity and layer thickness fields

The coefficients *h*^{o}_{n} are determined by the choice of a velocity field at the coast (*y* = 0), subject to the “static stability” constraint *h* = *H* + Δ*h* ≥ 0, in order to get physical solutions.

If we consider a constant background stratification and a perturbation consisting of a single mode, say *n*, then *H*(*ρ*) = const = *H _{o}*/Δ

*ρ*with Δ

_{o}*ρ*=

_{o}*ρ*

_{bot}−

*ρ*

_{surf}, and

In that case, the Richardson number can be written as

with

Note that the latter constraint on *α* ensures static stability.

When the Richardson number is given by Eq. (20) subject to the constraints in Eq. (21), it is easy to see that the minimum value is attained when *X* = 1 (this means that it is attained at *y* = 0 when the velocity is at its maximum) and *α* = ±1 (limit of static stability). Then Ri can be written as

Equations (21) and (22) show that the Richardson number is always larger than 1/2 and above the critical Richardson number of 1/4 for KH instability. As a result, Kelvin–Helmholtz instability cannot develop.

This result can be extended to more complex stratifications and currents. To do so, Eq. (11) can be solved given a (baroclinic) vertical profile for the velocity field at the coast. We found that the minimum Richardson numbers were obtained when the maximum value of the velocity field at the coast reaches the limit of static stability (in agreement with the analytical results *α* = ±1). Figures 2 and 3 show some numerical results with configurations where the vertical structures of the current consist of many different eigenmodes (see Figs. 2b and 3b) and more general background stratifications (see Figs. 2a and 3a). The numerical calculations show that the minimum Richardson number is again larger than 1/2 (see Figs. 2c and 3c).

Our results thus indicate that the Richardson number of a geostrophic flow without PV anomalies does not undergo Kelvin–Helmholtz instability. The reason why Ri cannot reach low values is associated with the fact that there is no PV anomaly. As stated above, the horizontal and vertical gradients of the velocity and density fields are not independent [see constraints in Eqs. (8)–(9)]. In fact, if one decreases (increases), the other must decrease (increase) too, and it is then not possible to reach the critical Richardson number for the current configurations discussed here (no PV anomaly, no outcropping, and constant stratification). This result is, however, only valid for configurations where interior isopycnals do not outcrop.

### b. Upwelling development in a two-layer model

To evaluate the possibility of obtaining Kelvin–Helmholtz instabilities when there exists outcropping of interior isopycnals, we look for the possibility of development of shear instabilities in a two-layer model when an upwelling develops. To do so, we first reevaluate the Howard and Miles criterion for this simplified configuration.

#### 1) Instability criterion

As in Howard (1961) and Miles (1961), we assume *f* = 0, *V* = 0, and all fields to be independent of *y*. In this case Eqs. (1) can be linearized around an initial state for which the fields *U*_{1}, *U*_{2}, *H*_{1}, and *H*_{2} are constant. Assuming the rigid-lid approximation (*h*_{2} + *h*_{1} = const in the continuity equations) and with perturbations of the form *ϕ* = *ϕ _{o}* exp[

*ik*(

*x*−

*Ct*)], we get (see also Ripa 1991)

where *g*′ = *g* Δ*ρ*/*ρ _{o}*, Δ

*ρ*=

*ρ*

_{2}−

*ρ*

_{1}is the density difference between the layers, and Δ

*U*=

*U*

_{2}−

*U*

_{1}. Thus, the system develops unstable waves if and only if

Note that this criterion is a necessary and sufficient condition for instability and is quite different from the general Howard (1961) and Miles (1961) criterion of Ri < 1/4 discussed above. It is, however, associated with a Kelvin–Helmholtz instability in a two-layer configuration. It has been assumed that the current is constant in each layer and the Coriolis effect has been neglected. These are strong assumptions that are not verified for upwelling currents (and the model developed further) so that criterion Eq. (24) is probably modified in the general case. It however remains representative of an instability associated with a vertical shear of the current and is the equivalent of the criterion obtained in the continuous stratification case. Our analysis will thus be based on criterion Eq. (24).

#### 2) Analytical solution for the upwelling

To analyze under which circumstances shear instability can grow during the development of an upwelling current, we now calculate the velocity fields when an upwelling develops and study whether criterion Eq. (24) is met.

For the dynamics of the upwelling development, we can again assume that all fields are independent of *x*. Considering an along-slope wind forcing, Eqs. (1) can be written

for *k* = 1, 2 where *M _{k}* is given by Eq. (3), PV

*= (*

_{k}*ζ*+

_{k}*f*)/

*h*is the PV in layer

_{k}*k*,

*T*

^{w}_{1}=

*T*is assumed constant, and

^{w}*T*

^{w}_{2}= 0.

Equations (25) can be solved provided

the along-slope current is in geostrophic equilibrium [In this case, the second equation in Eq. (25) simply becomes

*fU*= −∂_{k}. This approximation actually filters out the inertia–gravity waves and is well verified.] and_{y}M_{k}the PV is not modified during the upwelling development (PV

=_{k}*f*/*H*), which is true if_{k}*T*is constant.^{w}

In this case, Eq. (25) is a linear system of the unknowns (*U _{k}*,

*h*,

_{k}V_{k}*h*). With the rigid-lid approximation, the solution is (see appendix B)

_{k}where *δ* = *H*_{1}/*H*_{2}, *R _{d}* =

*g′H*

_{1}

*H*

_{2}/(

*H*

_{1}+

*H*

_{2}) is the Rossby radius of deformation, and

where *t _{o}* =

*fR*(l +

_{d}*δ*)/

*T*is the time at which the interface intersects the sea surface (at

^{w}*y*= 0).

#### 3) Kelvin–Helmholtz instability

From Eqs. (26) and (27), it is possible to show that the maximum vertical variation of velocity is reached at the outcropping front when the interface between layers 1 and 2 has reached the sea surface and is given by

When an upwelling develops, the shear instability criterion (24) can be met if

Thus, for an upwelling current in a two-layer model, Kelvin–Helmholtz instability only depends on the stratification at rest: it can develop when the initial upper-layer thickness is larger than the lower layer.

The physical explanation behind this result is that, as long as there is no PV transformation, the PV anomaly is zero (Δ*Q* = 0) and the current vorticity in each layer is then proportional to its thickness variation: *ζ* = *f* Δ*h*/*H* [see Eq. (6)]. In the vicinity of the outcropping front, layer 1 vanishes and its vorticity is then *ζ*_{1} ≃ −*f*, whereas layer 2 increases from *H*_{2} to *H*_{2} + *H*_{1}, gaining vorticity of the opposite sign, *ζ*_{2} ≃ *fδ*, that is all the greater in magnitude as *δ* is large. If the lower layer is initially associated with a small layer thickness, it will be dramatically stretched when the interface reaches the sea surface. This gives rise to very strong velocities in this layer, and strong vertical gradients can thus be achieved. These results are associated with our “PV thinking” approach and our hypothesis of no PV anomaly, which strongly constraints the vorticity field in both layers and the vertical shear. This is not the case in previous studies where it has generally been assumed that the lower layer was at rest (Barth 1989a, b; Paldor and Ghil 1991; Boss et al. 1996). In this case, the vertical shear is no longer driven by the stretching effects discussed here and criterion Eq. (29) does not apply.

This result can be extended to a more realistic model with more layers: small Richardson numbers can be achieved when two adjacent layers are respectively compressed and stretched, leading to the development of opposite-sign vorticities. The positive vorticity, and the vertical shear, is all the more important since the initial thickness of the layer that is stretched is originally small, which corresponds to the existence of a pycnocline. The most favorable moment for Kelvin–Helmholtz instability during the development of an upwelling is thus when the pycnocline reaches the sea surface (maximum compression of the upper layer, and maximum stretching of the layers within the pycnocline). This analysis is confirmed with numerical results, shown in Fig. 4 where the Richardson number along each isopycnal level has been plotted (third column) together with the density profile at rest (first column) and the stratification (layer thickness field, second column) on which contours of the Richardson number have been superimposed (within the range [0, 1] with contour step 0.1). Different configurations have been taken into account: in the first experiment, a constant stratification has been used with five layers outcropping at the sea surface; in all other experiments, the stratification is the same and a pycnocline is represented (in a layered model, this means the layer thickness at rest is small for some layers). In the latter case different stages of the upwelling development have been represented (note that the number of layers outcropping at the sea surface varies). Some of these flows have Richardson numbers lower than 1/4 and are thus likely to undergo KH instability. The region over which Ri is low is, however, pretty narrow (usually located in a narrow band close to the outcropping front, of a few kilometers width according to Fig. 4c).

## 4. Effect of the wind stress on PV

Calculating the PV evolution from Eqs. (1), one arrives at

Neglecting **F** and *S* yields the following equation for the influence of the wind stress on the PVA evolution:

Figure 5 represents a schematic side view of the upwelling system. The lightly shaded zone in the upper part of the ocean represents the region over which the wind stress extends: its depth is *H ^{w}* and

*T*is not zero in this region, but—for constant wind speed—it is constant along the geopotential surface

^{w}*z*= const (see Fig. 5a; also notice that as we have chosen a linear decrease of the wind stress,

*T*is actually constant over the whole shaded area). However, the dynamics of PV anomaly—as a material quantity—has to be followed along isopycnal surfaces (see McWilliams and Gent 1980; Hoskins et al. 1985). As long as isopycnals remain flat, no PV anomaly can be created but, as the upwelling develops, isopycnal surfaces bend upward and isopycnals initially located below

^{w}*z*= −

*H*will eventually enter the area affected by the wind stress (see Fig. 5b). Along such surfaces,

^{w}*T*is not constant and increases toward the coast:

^{w}*T*is negative in the shaded region for upwelling and is zero below. Thus, rot(

^{w}**T**

^{w}) = −∂

_{y}

*T*is negative, which leads to the creation of a negative PV anomaly.

^{w}_{x}For a downwelling system, a similar analysis shows that, again, a negative PV anomaly is created near the coast [the wind stress is of opposite sign, but the isopycnals bend downward, so that we recover the same sign for rot(**T**^{w})]. Also notice that the rate of creation of the PV anomaly is proportional to the inverse of *h*, so we can expect the PV anomaly formation rate to be very large for upwelling and moderate for downwelling.

Last, let us mention that in our simplified two-layer configuration the previous mechanism is qualitatively well represented. In fact, specifying that the wind stress acts over a given depth leads to

for the first layer (as long as its thickness is larger than *H ^{w}*). Then rot(

**T**

^{w}) is indeed negative for both upwelling and downwelling currents, and the rate of PV anomaly creation can also be roughly estimated. Indeed, Eq. (31) yields

assuming that *δh* ≃ *H* (well-developed upwelling or downwelling) over a length scale *δy* ≃ *R _{d}*, which then yields

It is interesting to compare Δ*Q* with the Coriolis parameter *f*; with *τ*^{surf} ≃ 0.3 N m^{−2} and *ρ* ≃ 1000 kg m^{−3}, Eq. (34) becomes

where *δt* is expressed in days. For upwelling currents *h* < *H* and the PV formation rate is then quite large when the upwelling is well developed. For instance, in the region where *h* is only one-third of the layer depth at rest, *H*, the rate at which the PV anomaly is created is at least *O*( *f* ) within one day of wind forcing. For downwellings, the formation rate is about two orders of magnitude weaker because *h* > *H*. As a result, fairly intense negative PV anomalies can be formed for upwelling but stay moderate for downwelling. This negative PV anomaly is then able to interact with the positive PV anomaly associated with the outcropping front when the interface reaches the sea surface (or bottom for downwelling currents), and leads to barotropic/baroclinic instabilities.

## 5. Numerical results

### a. Generalities

Both Kelvin–Helmholtz instabilities and the effect of the wind stress can modify the PV field and provide the negative PV anomaly required to explain the observed instability of upwelling currents. To test these results and evaluate the most efficient mechanism, we use a numerical model with two layers (see Fig. 6). The density difference between the layers is fixed and (*ρ*_{1} − *ρ*_{2})/*ρ*_{1} = 1‰. The layer thicknesses at rest are *H*_{1} and *H*_{2}. Numerical simulations are performed in an east–west periodic channel of 300 × 300 km^{2} (see Fig. 6) with a grid step Δ*x* = 2 km. The Coriolis parameter is *f* = 7 × 10^{−5} s^{−1}. In most of this paper, we will analyze results from one configuration where we have set *H*_{1} = 500 m, *H*_{2} = 400 m, and *H ^{w}* = 10 m, for which

*H*

_{1}>

*H*

_{2}and we can expect the upwelling to develop Kelvin–Helmholtz instabilities. These parameters yield an internal Rossby radius

*R*≃ 21 km. The ocean is initially at rest and a wind stress is applied over a depth

_{d}*H*(see Fig. 6) at

^{w}*t*= 0. The wind stress at the surface is (

*τ*,

_{x}*τ*) = (−0.3, 0.0) N m

_{y}^{−2}(corresponding to a westward wind with speed

*W*≃ 5–10 m s

^{−1}). It decreases linearly over

*H*and is zero below

^{w}*z*= −

*H*. An upwelling current thus develops along the southern coast, and a downwelling current along the northern one.

^{w}Figure 7 represents a horizontal map of the first-layer thickness *h*_{1} at time steps *t* = 46, 47, 50, 52, 54, and 56 days. Note the different stages during the evolution. Initially the upwelling and downwelling currents develop along the southern and northern boundaries, respectively, and both currents keep a two-dimensional structure. At *t* = 46 days, the interface between layers 1 and 2 reaches the surface at the southern coast, and an outcropping front is formed, while at the northern coast the downwelling has not led to outcropping at the bottom at this time. Then a short-wave instability with horizontal scale *L* ≃ 5–10 km and a typical growth rate *τ* ≃ 2–4 *h* (estimated from model output at 2-h intervals that are not shown) appears along the upwelling front. The induced perturbations are clearly visible at *t* = 47 days in Fig. 7 when they have grown to a stage where they are associated with steep waves. Shocks form that are then smoothed by viscosity and diffusion, and seem to disappear. A long-wave instability then appears with a growth rate *τ* ≃ 2 days, close to that found by Boss et al. (1996). The latter gives rise to eddies that are clearly visible and detached from the coast at *t* = 56 days, with a typical diameter of *D* ≃ 50 km. Note also the filament linking the upwelling current and the detached eddies. In this simulation, the filament is rather short, but other configurations can lead to longer filaments whose extension can reach 200 km.

### b. Kelvin–Helmholtz instability

Figure 8 shows the evolution of Δ*U* in the vicinity of the southern (solid line) and northern (dashed–dotted line) coast; the critical shear Δ*U*_{crit} = 3 m s^{−1} is also plotted (dashed line). Note that the instability criterion is met for the upwelling current (solid line) at about *t* = 36 days, whereas it is never met for the downwelling current. It thus seems that the initial short-wave instability that we have observed is of Kelvin–Helmholtz type. It does, however, not closely follow criterion Eq. (24) (notice the instability only starts when the vertical shear strongly exceeds the critical level), which is not surprising considering that the Coriolis term, the layer thickness variations, and horizontal shear of the current were neglected. Analytical calculations of the growth rate of short-wave (and long wave) instabilities have been performed with a frontal current with no PV anomaly in the upper layer and no velocity in the lower layer. The association of the short-wave instability with Kelvin–Helmholtz instability has been identified in previous studies (see, in particular, Paldor and Ghil 1991) and the influence of the vertical shear and total depth is in agreement with criterion Eq. (24): weak vertical shear and large total depth both limit the development of the instability.

Other instability types may also account for the growth of these short-wave instabilities: barotropic/baroclinic instabilities and inertial resonance (see Klein and Tréguier 1993; Tréguier and Klein 1994). We believe that the former is unlikely because this instability usually results in the formation of eddies and generally has smaller growth rates and also because its trace is much less well marked in the PV field (see Fig. 11), which is a feature of inertia–gravity waves. Klein and Tréguier (1993) (see also Tréguier and Klein 1994) have shown that, when there exists a current with horizontal shear, inertial oscillations can extract energy from the wind and are amplified if they propagate against it. The energy source for inertial resonance is the wind, and the mechanism leads to the generation of small-scale instabilities: it is inhibited by diffusive effects and is favored when the mixed layer is very shallow. It could thus also explain the observed results. We have repeated the previous numerical experiment but with the energy source (wind) turned off at *t* = 46 days. The same instability still develops in this experiment, which proves that inertial resonance does not account for the observed instability and that the nature of the short-wave instability observed in our numerical results is of Kelvin–Helmholtz type.

### c. Baroclinic instability

Figure 9 shows the PV anomaly profiles across the channel at different times (*t* = 0, 10, 20, 30, and 40 days) and for the top layer, and Fig. 10 also pertains to times *t* = 45 and 46 days but represents a magnified view of the southern coast. Negative PV anomalies develop at the southern and northern coasts, with the strongest values for the upwelling current, along the southern boundary. As soon as the interface reaches the surface (*t* = 46 days), a strong positive anomaly appears (note that its value has been kept finite by imposing a lower boundary for *h* when calculating Δ*Q*). The PV anomaly in the lower layer is much smaller than the upper-layer one and does not seem to play a significant role in the dynamics of the instability. This is in agreement with the previous analytical results and proves that the wind stress is responsible for the development of the negative PVA. Indeed, it is created before the Kelvin–Helmholtz instability develops, and modifying the viscosity coefficient does not change the results; so isopycnal mixing or viscous boundary layer effects (see Morel and McWilliams 2001) cannot be responsible for it.

Figure 11 represents the evolution of Δ*Q* in the first layer at the same times as in Fig. 7 (*t* = 46, 47, 50, 52, 54, and 56 days). Negative values of PV anomaly are associated with dashed contours and positive ones with plain contours. From *t* = 47 to 50 days Kelvin–Helmholtz instabilities develop, with inertia–gravity waves that become quite steep, forming fronts that are then smoothed by viscous effects. This process alters the PV anomaly field and actually mimics (although very imperfectly) the mixing associated with Kelvin–Helmholtz instability.

At *t* = 50 days, the positive PV anomaly associated with the outcropping front is well developed and larger-scale perturbations are clearly visible (notice that, even though they remain weak at this stage, the long-wave instability can be seen as soon as day *t* = 47). They grow and lead to the formation of a dipole detaching from the coast at *t* = 56 days, which is typical of barotropic/baroclinic instability. Notice that the positive PV anomaly is initially infinite as it is associated with the outcropping front. Again, viscous effects lead to the formation of a cyclonic eddy with finite PV anomaly where the interface between layers 1 and 2 no longer outcrops at the surface.

### d. Sensitivity to Kelvin–Helmholtz instability

The previous experiment confirms the analytical results obtained above but does not permit us to examine the relative importance of both the Kelvin–Helmholtz instability and the wind stress on the transformation of the PV field and subsequent instability. Kelvin–Helmholtz instability indeed leads to mixing, which can also rapidly modify the PV structure of the upwelling current and make it baroclinically unstable (see Morel and McWilliams 2001). On the other hand, we have shown that, in general, the region where Kelvin–Helmholtz instability is possible in upwelling systems is along a narrow band located near the upwelling front, and one may thus wonder whether mixing in such a restricted area can lead to substantial changes in the PV structure and have some impact on the barotropic/baroclinic stability properties of the current.

Choosing a configuration where *H*_{1} < *H*_{2} should lead to the development of a current that is stable with respect to Kelvin–Helmholtz instability. This result was based, however, on the hypothesis of a constant wind stress in the first layer so that PV was conserved, which, as seen above, is not the case in our configuration. Tests with *H*_{1} slightly below *H*_{2} have shown that a Kelvin–Helmholtz instability still develops in these cases so that the PV transformation associated with the wind stress, which we have described, increases the vertical shear and the possibility of Kelvin–Helmholtz instability [Eq. (6) indeed shows that adding a negative PV anomaly in layer 1 leads to stronger negative vorticity when the layer is compressed]. However, when *H*_{2} is much larger than *H*_{1}, the Kelvin–Helmholtz instability is inhibited when upwelling develops.

Figure 12 represents the PV anomaly in the first layer for different times for a configuration identical to the previous one except that *H*_{2} = 1500 m now. In this case the critical velocity shear necessary for Kelvin–Helmholtz instability to develop is Δ*U*_{crit} = 4.5 m s^{−1}, which is never achieved during the simulation. No short-wave instability appears to grow, and we thus believe that any Kelvin–Helmholtz instability was absent from this simulation. A strong negative PV anomaly however develops rapidly in the upwelling current, thanks to the mechanism described above, and when the interface between layers 1 and 2 eventually outcrops, a long-wave barotropic/baroclinic instability develops. This is clear from Fig. 12, which presents many features common to Fig. 11. Kelvin–Helmholtz instability, and the subsequent mixing associated with this process, is thus not necessary to achieve an unstable PV structure. The effect of the wind stress and the formation of an outcropping front is sufficient to create an unstable structure and generate eddies.

Kelvin–Helmholtz instability can, however, modify the PV structure of the flow and change the instability characteristics when it exists. An indication of this is given in Fig. 10, which shows the PV anomaly structure in the first layer just before and just after the development of the Kelvin–Helmholtz instability and formation of the outcropping front. The latter should only be associated with an additional high-PV anomaly near the coast, and the negative part should roughly be conserved (or become even stronger). This is not the case, and the negative PVA has been drastically modified (it becomes weaker) in the vicinity of the upwelling front where the Kelvin–Helmholtz instability has yielded strong perturbations. As a result, the PV gradients are transformed and the subsequent barotropic/baroclinic instability is thus likely to be affected.

## 6. Discussion

Many previous studies on the instability of upwelling currents were based on the hypothesis of no flow in the lower layers. In such cases the Charney–Stern criterion for instability is always met; in a two-layer system for instance, the front is associated with a positive PV gradient and the lower layer with a negative PV gradient. The current thus develops long-wave barotropic/baroclinic instabilities.

We have adopted a different perspective and addressed what establishes the PV structure of the current. We have shown that thermobaric or diabatic processes are required to explain the observed barotropic/baroclinic instabilities, and have examined two particular mechanisms more closely: 1) Kelvin–Helmholtz instability and 2) wind stress. With regard to Kelvin–Helmholtz instability, we calculated the Richardson number of geostrophic currents with Δ*Q* = 0 (no PV anomaly) and no pycnocline. We showed that the Richardson number is always above 1/2, indicating stability. A two-layer formulation suggests that low Richardson numbers occur near upwelling fronts when *H*_{1} > *H*_{2}. In multilayer configurations, this corresponds to the existence and outcropping of the pycnocline. These results are a consequence of assuming a zero PV anomaly, which constrains the vertical shear. Our results with two layers contrast with previous studies in which short-wave instabilities were observed when *H*_{2} was larger than *H*_{1}. In fact, these studies usually assumed a lower layer at rest, and the lower-layer thickness thus has no influence on the vertical shear.

We observed that the Richardson number is subcritical in a narrow region near the front (a few kilometers wide at the most). The Kelvin–Helmholtz instability in this region will manifest itself as short-wave instabilities, leading eventually to diapycnal mixing. This modifies the PV, thus allowing for the possibility of long-wave instability (e.g., Morel and McWilliams 2001).

We then studied the effect of wind stress. Upwelling causes isopycnal surfaces to bend upward closer to the surface, bringing them into the region forced by the wind. This causes wind stress variations along the isopycnal surface, permitting modification of the PV. Note that this occurs even when the wind has no curl. We have shown that a negative PV anomaly is created, for both upwelling and downwelling; the anomaly however is generally much stronger for upwelling. When the front outcrops, the CS criterion is met and the current destabilizes, leading to eddy and filament formation. The wind thus appears to be an efficient means of destabilizing upwelling currents.

This is of course valid only for upwelling or downwelling currents generated by local wind forcings along a coast. In the ocean there exist upwelling systems that are remotely controlled and established through the propagation of coastal Kelvin waves so that the local wind plays little role in their formation. For instance, in equatorial regions, equatorial waves sometimes generate upwelling currents when they reach a coast (see, e.g., Adamec and O'Brien 1978; Clarke 1979; Houghton 1989). In this case, the wind forcing is absent, and the previous conclusions no longer hold. Baroclinic/barotropic instabilities are thus only possible if a PV source can be identified that is different from the wind stress. Kelvin–Helmholtz instability could then be the primary source of instability for this type of upwelling. It would be interesting to compare the dynamics of such upwelling systems with the dynamics of systems controlled by local winds as we could expect very significant differences as far as their (in)stability characteristics are concerned.

Short- and long-wave instabilities are frequently observed with developed upwelling currents (Beardsley and Lentz 1987; Brink and Cowles 1991; Kosro et al. 1991; Strub et al. 1991; Dewey et al. 1991; Ramp et al. 1991; Haynes et al. 1993; Relvas de Almeida 1999). Early studies sometimes associated the small-scale fingers along the upwelling front with the longer filaments developing later. But recent numerical and theoretical studies have shown that these instabilities are different in nature (see McCreary et al. 1991; Barth 1994; Boss et al. 1996; Roed and Shi 1999; Shi and Roed 1999). The small-scale instability has a rapid growth rate and has sometimes been associated with Kelvin–Helmholtz instability (see, in particular, Paldor and Ghil 1991). It causes fingerlike variations along the upwelling front that have weak extensions. The long-wave instability is associated with barotropic/baroclinic instability and generates deformation-scale eddies and filaments. These structures can transport upwelled waters hundreds of kilometers offshore.

Our results also apply to downwelling currents. For the two-layer case, it can be shown that Kelvin–Helmholtz instability develops if *H*_{1}/*H*_{2} < 1. In the open ocean, the bottom layer is generally too deep to intersect the ocean floor. Kelvin–Helmholtz instability for downwelling currents thus seems likely only in shallow areas. In addition, the PV anomalies generated by wind forcing are weaker for downwelling currents than for upwelling currents. As such, downwelling currents are the more stable. However, a continental shelf will act to overcome this difference, enabling the intersection of the isopycnal with a boundary.

Among the other mechanisms that could influence the stability properties of these currents is isopycnal mixing associated with viscous boundary layer effects. In our experiments we minimized this by an appropriate choice of viscosity coefficient and boundary conditions. Viscous parameterizations remain a delicate issue in numerical models and require more study (see Morel and McWilliams 2001).

The evolution of fronts is a delicate numerical issue in isopycnal models, as it deals with infinitesimal layers (Bleck and Boudra 1986; Bleck and Smith 1990; Bleck et al. 1992). We thus ran identical experiments with another isopycnic model [R. Hallberg's isopycnic model (HIM; see Hallberg 1995, 2000)]; this model treats vanishing layer thickness differently than does MICOM. Prior to the development of barotropic/baroclinic instability, the two models yielded the same results (very similar PV evolution in particular). After the onset of instability, the details differed, but the overall evolution was qualitatively the same. So the models produced eddies with similar growth rates, size, and PV structure.

The two-layer configuration is not realistic and increasing realism by taking into account more layers, bottom topography, a better parameterization for the wind stress, or mixing associated with Kelvin–Helmholtz instability, etc. has to be undertaken. In particular, as one of the major source of PV generation and destabilization of upwelling currents is associated with wind stress effects, the present results can—at least quantitatively—be sensitive to the parameterization chosen to take these effects into account. This again calls for more realistic simulations and comparison of different parameterizations of the wind stress effect.

Last, the present results apply in principle to other situations, such as the effect of wind on the PV structure of jets or vortices.

## Acknowledgments

The present study was made possible by the use of the MICOM and HIM (Dr. Robert Hallberg's isopycnic model) codes. The authors are very grateful to Drs. Rainer Bleck, Eric Chassignet, Linda Smith, and Robert Hallberg for providing their codes. The authors are also indebted to Dr. Joseph Lacasce whose suggestions and comments helped to improve the paper. This paper is a contribution to the research program Modélisation Océanique d'un Théâtre d'Opérations Navales conducted by SHOM (the French Navy's oceanographic institute) and funded by the French Ministry of Defense (Grant PEA012401).

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Kelvin–Helmholtz Instability Criterion

In their papers Miles (1961) and Howard (1961) considered fairly general equations with possible nonhydrostatic effects. In this appendix we simply reproduce Howard's calculations in the framework of the primitive equations written in isopycnal coordinates to show that the latter system of equations contains the physics of the early stages of Kelvin–Helmholtz instability development. For this instability, we can neglect the Coriolis effect; if we linearize Eqs. (1)–(2) around a stationary state *U* = *U _{o}*(

*ρ*),

*h*=

*H*(

_{o}*ρ*), we get

If we seek solutions of the form *ϕ* = *ϕ*(*ρ*) exp[*ik*(*x* − *Ct*)], we can derive an equation for the velocity perturbation:

Using *N* ^{2} = *g*/*ρ _{o}H_{o}*, and defining

*W*= (

*U*−

_{o}*C*) and

*u*=

*W*

^{−1/2}Φ, we obtain an equation for Φ:

Multiplying (A3) by Φ*, the complex conjugate of Φ, and integrating over the whole density range ([*ρ*_{surf}, *ρ*_{bottom}]) with appropriate boundary conditions, we get

Defining *C* = *C _{r}* +

*iC*and retaining only the imaginary part of the latter equation, we get

_{i}Thus, if *C _{i}* ≠ 0, necessarily

which yields Miles's and Howard's necessary condition for instability written in isopycnal coordinates:

### APPENDIX B

#### Upwelling Development in a Two-Layer Model

With the hypothesis given in section 3b(2) and the rigid-lid approximation, Eqs. (25) become

This system is linear if we consider the unknowns (*U*_{1}, *h*_{1}*V*_{1}, *h*_{1}, *U*_{2}, *h*_{2}*V*_{2}, and *h*_{2}). Using Eq. (3) and recalling that *T _{w}* = const here, the system can be transformed into

with *R _{d}*

^{2}=

*g*′

*H*

_{1}

*H*

_{2}/

*f*

^{2}(

*H*

_{1}+

*H*

_{2}). Equation (B2) can be solved distinguishing two different time periods:

*t*∈ [0,

*t*] and

_{o}*t*∈ [

*t*, ∞], where

_{o}*t*is the time at which the interface reaches the sea surface (see Fig. B1a). We also define the position where the interface outcrops as

_{o}*y*=

*Y*(

*t*) (see Fig. B1b).

For *t* ∈ [0, *t _{o}*], the general solution for

*U*

_{1},

*U*

_{2},

*h*

_{1}, and

*h*

_{2}for

*t*∈ [0,

*t*] can be written

_{o}where *δ* = *H*_{1}/*H*_{2}. The first and second equations in Eq. (B2) then permit one to find *U _{c}*(

*t*) and

*U*(

_{b}*t*) with the additional constraint that

*V*

_{k}_{=1,2}(

*y*= 0,

*t*) = 0. We thus find

From these results *t _{o}* can be calculated, as the interface reaches the sea surface when Δ

*h*

_{1}(

*y*= 0,

*t*=

*t*

_{0}) = −

*H*

_{1}, which yields

To solve Eq. (B2) for the time period *t* ∈ [*t _{o}*, ∞], we repeat the previous steps, distinguishing between the interval

*y*∈ [−∞,

*Y*(

*t*)] and

*y*∈ [

*Y*(

*t*), 0]. After some manipulations, we get

Here *U _{c}*(

*t*),

*U*(

_{b}*t*), and

*Y*(

*t*) are then calculated noting that

*h*

_{1}

*V*

_{1}[

*y*=

*Y*(

*t*),

*t*] = 0 and

*h*

_{2}

*V*

_{2}(

*y*= 0,

*t*) = 0. We then get

## Footnotes

*Corresponding author address:* Yves Morel, EPSHOM/CMO, BP 30316, Brest Cedex 29603, France.. Email: morel@shom.fr

^{1}

If instead of an ocean at rest, the upwelling develops above an initial current with negative PV gradient; then, when isopycnals intersect the sea surface and form a positive PV gradient, the CS criterion is verified and the current can become unstable. The influence of a preexisting current with negative PV gradients is likely to play a role in some regions (such as the Portuguese coast where upwellings develops above the Mediterranean Sea outflow, which is a current associated with low PV values).

^{2}

Notice that, in isopycnal coordinate models, the velocity field **U** is along geopotential surfaces, while derivatives are taken along isopycnals.