The purpose of this paper is to understand how long planetary waves evolve when propagating in a subtropical gyre. The steady flow of a wind-driven vertically sheared model subtropical gyre is perturbed by Ekman pumping that is localized within a region of finite lateral extent and oscillates periodically at about the annual frequency after sudden initiation. Both the background flow and the infinitesimal perturbations are solutions of a 2½-layer model. The region of forcing is located in the eastern part of the gyre where the steady flow is confined to the uppermost layer (shadow zone). The lateral scales of the forcing and of the response are supposed to be small enough with respect to the overall gyre scale that the background flow may be idealized as horizontally uniform, yet large enough (greater than the baroclinic Rossby radii) that the long-wave approximation may be made. The latter approximation limits the length of time over which the solutions remain valid. The solutions consist of (i) a forced response oscillating at the forcing frequency in which both stable (real) and zonally growing (complex) meridional wavenumbers are excited plus (ii) a localized transient structure that grows as it propagates away from the region of forcing. Application of the method of stationary phase provides analytical solutions that permit clear separation of the directly forced part of the solution and the transient as well as estimation of the temporal growth rate of the transient, which proves to be convectively unstable. The solutions presented here are relevant to understanding the instability of periodic (including annual period) perturbations of oceanic subtropical gyres on scales larger than the baroclinic Rossby radii of deformation.
a. Background and principal results
Chelton et al. (2004) analyzed four-year averages of 25-km-resolution measurements of near-surface wind speed and direction and showed that they reveal the surprisingly persistent small-scale features in ocean winds. Superimposed on these persistent wind patterns are temporally varying perturbations with magnitudes comparable to the mean winds. The strongest of these perturbations is the annual cycle. Chelton and Schlax (1996) correspondingly documented energetic annual-period Rossby wave sea surface elevation perturbations.
Motivated by these observations, the goal of this paper is to explain in detail how perturbations of a particular background subtropical gyre flow develop in space and in time when forced by winds characterized by scales much smaller than the gyre scale but much larger than the baroclinic Rossby radii of deformation that are started up at some initial instant and are thereafter periodic in time.
The solutions found here thus consist of a transient component and a component harmonic in time at the forcing frequency. Within the confines of the approximations and restrictions noted below, it is found that both the initial transient and the co-oscillating periodic response to spatially localized, suddenly initiated periodic Ekman pumping may be completely described by approximate analytical techniques.
The most important results are as follows. The perturbations excited by a patch of Ekman pumping that is localized within the eastern region of a subtropical gyre and oscillates periodically after being suddenly initiated consist (i) of a directly forced solution that includes a well-defined beam of waves oscillating at the forcing frequency and radiating westward away from the forcing region without spatial growth plus a further harmonically oscillating part growing spatially toward the southwest as well as (ii) of a localized transient structure that grows as it propagates away from the region of forcing. Analytical approximate solutions permit clear separation of the forced and the transient solution as well as estimation of the growth rate of the transient, which proves to be convectively rather than absolutely unstable (Huerre and Monkewitz 1990), meaning that growth occurs only for certain ranges of x/t, y/t as |x|, |y|, |t| → ∞, rather than for all x/t, /y/t.
b. Previous studies
Chelton and Schlax (1996) suggested that purely periodic (annual period) wind-driven long planetary waves can become energized in the western part of ocean basins. Explanations of this in the literature involve either (i) modification by midocean ridge topography of long waves radiated from the eastern boundary (Barnier 1988; Tailleux and McWilliams 2002) or else (ii) baroclinic instability (Liu 1999a; Dewar and Huang 2001; Kubokawa and Nagakura 2002). This paper further explores (ii) the instability. Previous studies in the literature that analyzed stability of long baroclinic waves in the subtropical gyre are as follows. Liu (1999a) found a linearized planetary wave instability in the southwestern part of a subtropical gyre that was associated with coupling, by the background subtropical gyre flow, of the fast and slow modes of a 2½ layer QG model. Liu (1999a) perturbed the steady flow of an eddy-resolving 3-layer QG ocean model with spatially uniform 3-yr period Ekman pumping and correspondingly found the 3-yr period ocean response to be “almost unchanged in longitude until it enters the southwestern part of the subtropical gyre, where linear eigenvalue calculations show the strongest planetary wave instability.” He suggested that this response was related to the possibility that “the temporal instability of an eigenmode (complex frequency with real wavenumber) may correspond to a spatially unstable radiating mode that is forced by a given frequency (real frequency with a complex wavenumber).” Liu (1999b) as well as Dewar and Huang (2001) carried out similar model perturbation experiments but with a more complex and partially ventilated subtropical gyre flow. Dewar (1998) analyzed the effect of the mean flow on the phase speed of long baroclinic waves in a subtropical gyre. Kubokawa and Nagakura (2002) studied the evolution, in both ventilated and unventilated regions, of disturbances generated by more localized initial disturbances and/or regions of Ekman pumping, both by the Fourier method and numerically. Even in the present idealized case, with no horizontal current shear, Walker and Pedlosky (2002) showed that, with nonzero purely meridional background flow confined to the upper layer of a 2-layer QG model, there was no minimum vertical shear required for instability so that instability need not be confined to any particular region of the subtropical gyre.
c. Restrictions and approximations
Ocean circulation varies over a range of time scales from days to many years. Motivated by a desire to understand the effects of periodic forcing on the underlying and much more slowly varying circulation, the numerical model studies quoted above have typically been carried out by generating a numerical solution retaining all terms including nonlinear ones, first with steady forcing and then additionally imposing the periodic forcing (e.g., Liu 1999b; Dewar and Huang 2001). In contrast, the accompanying supporting and interpretive analytical studies linearize the periodic flow about a background flow that is assumed steady [Liu 1999a, b; Dewar and Huang 2001: the single nonlinear analytical exception known to us is Dewar (1989)]. The solutions of this paper are likewise linearized, about the particularly simple background flow described below.
The overriding motivation for linearization is, of course, the resulting great simplification of analysis. Linearization generally fails when the amplitude of the periodic forcing/flow becomes large, yet often provides insight into the structure of more complete solutions. Thus in the model studies quoted above, the authors show that linearized solutions are very useful in understanding the spatial propagation and evolution of the—in principle nonlinear—perturbations of model solutions that are excited by spatially localized periodic forcing. Note however that it is difficult to assess from these studies the limits of validity of the linearized solutions for the ocean, partly because these model studies do not attempt to closely match the strength of the periodic forcing to that acting on the real ocean.
Investigation of the linearized solutions obtained below was initially motivated by the observation of persistent patterns in wind forcing on scales smaller then the gyre scales (Chelton et al. 2004), and furthermore we focus on the annual period signal. However, it is important to note that the annual period solutions presented here are readily scaled to different periods, as explained in section 3, and the effects on the sympathetic part of the solution of changing the forcing period are discussed in section 4.
Because the annual frequency is one of the most strongly energized frequencies in the ocean, linearization may be less appropriate than at other frequencies. In defense of linearization at the annual period, it is appropriate to cite the qualitative correspondence between simple linear Rossby wave models driven by annual winds and the sea level signal documented by Chelton and Schlax (1996) and by Wang et al. (2001) in the Indian Ocean. In the final analysis however, particularly for the annual period in the present baroclinically unstable case, linearization must be viewed as a preparation for more complete analysis rather than as definitive.
In this paper, attention is restricted to unventilated background flows. Kubokawa and Nagakura (2002) compare free and forced solutions in the ventilated and unventilated regions of a 2½ layer model subtropical gyre and find that the stability properties of the linearized long-wave solutions are very different between the two regions. This makes it difficult to simply anticipate the results of extension of the present analysis to ventilated regions; such an analysis must be left for a subsequent study.
The most important remaining approximations are (i) the assumption that the dominant scales of the solutions are much smaller than the lateral scales of the subtropical gyre so that the gyre flow may be idealized as spatially constant and (ii) the long-wave approximation. It is important to note that i precludes exploration of the effects associated with lateral variation in the vertical structure of the modes. The long-wave approximation ii restricts the analysis to scales substantially larger than the largest baroclinic Rossby radius of deformation. To simultaneously satisfy both i and ii, the initial scales are chosen to be a good deal smaller (hundreds of kilometers) than the scales characterizing the gyre (thousands of kilometers) yet a good deal larger than the baroclinic Rossby radii (tens of kilometers) and solutions are presented only for times sufficiently small that scales on the order of Rossby radii are not yet strongly energized. While it is difficult to rigorously satisfy these two requirements in a realistic subtropical gyre, it is meaningful to investigate the solution dynamics in detail since, on one hand, the long-wave model is widely used to discuss basic ocean dynamics and, on the other hand, there exist persistent wind features on scales of several hundred kilometers (Chelton et al. 2004).
It has been shown that the fastest-growing waves in the subtropical gyre flow are at the scales of the Rossby radius of deformation (e.g., Spall 1994). This study provides detailed analysis of baroclinic instability acting on larger scales at which the planetary geostrophic approximation is valid. A subsequent and more complete analysis would include both the gyre scale and the Rossby radius of deformation scales simultaneously. Because of the complexity of this problem, it would have to be addressed numerically. The present detailed analytical study should thus be viewed as necessary preparation for a more complete study that also includes smaller scales.
Theories of large-scale ocean flow commonly make use of the planetary geostrophic (PG) approximation in which the open ocean flow is geostrophic except for a surface Ekman layer. A 2½-layer PG model (appendix A) is here supposed to govern the wind-driven steady circulation of a subtropical gyre in a rectangular midlatitude basin. In this paper we consider small amplitude time-dependent perturbations of such a steady flow that are driven by specified time-dependent surface Ekman pumping w′E(x, y, t). An important feature of the steady solutions considered here is that, since the deeper layers do not surface to be exposed to the wind, the steady flow is confined to the upper layer. The 2½-layer equations for time-dependent perturbations η1 and η2 of the upper and lower interfaces between the three layers of the background flow are (appendix A)
In these equations and subsequently, subscripts x, y, and t denote differentiation. For notational simplicity we have introduced the abbreviations
in which H0j is the unperturbed depth of the base of layer j (indices j = 1, 2, 3 correspond to the upper, middle, and lower layer), the notation H20 ≡ H02 emphasizes that the depth of the bottom of the unperturbed stagnant middle layer is constant, f (y) is the Coriolis parameter, the γj are reduced-gravity parameters defined in appendix A, and
where ug and υg are the mean upper layer geostrophic velocity components −γ1H01y/f and γ1H01x/f in the zonal x and meridional y directions and cR = C − (γ1/γ2)c.
3. The initially forced problem
Suppose that both interfacial displacement perturbations η1 and η2 are initially zero. Spatially localized Ekman pumping w′E is turned on at the initial time t = 0 and thereafter oscillates at the fixed frequency σ0. We call this the initially forced problem. The resulting solution consists of a transient part that satisfies a homogeneous initial value problem and a sympathetic part that oscillates at the frequency of the forcing.
We assume that the gyre scale over which the steady background flow varies horizontally is substantially greater than the horizontal scale of the perturbations. Under this assumption c, C, UR, and VR are taken to be constant. For validity of the underlying PG approximation, the horizontal scale of the perturbations must remain larger than the baroclinic Rossby radii characterizing the stratification. This requirement ultimately limits the length of time over which our solutions remain valid approximations.
The governing equation for the time-dependent interface displacement of the middle layer η2 is then, from (1),
and η1 is given in terms of η2 by either of (1). We ultimately consider a Gaussian patch of Ekman forcing of width L
The initial conditions are
Note that solutions η2(x, y, t) of (4), (5), and (6) are also solutions η2(sx, sy, st), where s is any real scale factor, provided that σ0, W0, and L are scaled by σ0/s, W0/s, and sL. The solutions presented subsequently are thus readily scaled up or down in space and time without change of form (providing the scales remain larger than the Rossby radii).
We Fourier transform (4) in space using the conventions
where Ê(k, l) is the Fourier transform of the Gaussian disturbance (5)
The solutions for η̂1 and η̂2 with Fourier transformed initial conditions (6) consist of a transient (T) part and a sympathetic (S) part that oscillates at the forcing frequency σ0:
where σ+(k, l) and σ− (k, l) are the two solutions of the dispersion relation for plane wave solutions of the homogeneous version of (9),
and, from (1), Σ0 is
The properties of free solutions with dispersion relation (13) have already been discussed by Liu (1999a,b) and by Kubokawa and Nagakura (2002). The two roots of (13) correspond to the fast and slow baroclinic mode, or the “Non Doppler Shift” (N) and the “Advective” (A) modes of Liu (1999a), which reduce to the first and the second internal baroclinic mode in the case of no background flow. When the solutions become unstable, their vertical profiles are no longer necessarily orthogonal, but the decomposition into fast and slow modes implicit in the present method of solution remains unique. The wavenumber dependence of the fast root is summarized in Fig. 1. The properties of the solutions necessary for understanding the evolution of the solutions of the initially forced problem will be invoked as needed in the subsequent text.
The transient solution η̂2T satisfies the homogeneous version of (9), so that η̂1T and η̂2T can be written in the form:
in which Σ+ and Σ− are given, from (1), by
The coefficients A and B are obtained from the initial conditions.
The corresponding transient and sympathetic solutions for η2 in physical space are finally given by
Similar expressions for η1T and η1S are readily obtained from these:
4. Solution of the initially forced problem
a. Numerical solution for transient and sympathetic components
This section discusses solutions η1(x, y, t) and η2(x, y, t) of (1) initiated from rest by Ekman pumping of the Gaussian form (5) initiated at time zero and subsequently periodic with annual frequency σ0. These have been evaluated both by summation of Fourier components in a periodic domain whose lateral extent is much greater than that of the forcing and by time stepping the finite difference versions of (1) forward in the same domain. In the Fourier method, the initial conditions and forcing are decomposed into Fourier components ei(kx+ly), each Fourier component is independently advanced in time according to its dispersion relation σ = σ(k, l) by multiplication by e−iσ(k,l)t, and the Fourier components are then recombined as in (17) through (20) to give the solutions at time t.
The background flow is that prevailing in the eastern part of a steady and unventilated subtropical gyre occupying a midlatitude basin of meridional extent 3000 km and zonal extent 10 000 km, driven by steady Ekman pumping of the form (52). The resting depths (before application of the steady Ekman pumping) of the base of two upper layers of the model H10 and H20 are chosen to be 500 and 1000 m with reduced-gravity parameters (appendix A) γ1 = 0.0245 and γ2 = 0.01 m s−2. At the location at which the time-dependent perturbation wind patch (5) is centered (1000 km west of the eastern boundary of the basin and 1000 km north of its southern boundary), the corresponding wave scale speeds C and c [(2)] are, respectively C = 4.60 and c = 3.23 cm s−1. For the particular example shown here, the coefficients UR and VR appearing in (1) have the values UR = −2.3 cm s−1 and VR = −1.9 cm s−1. The zonal component of the geostrophic velocity at this location is ug = −2.8 cm s−1 and the meridional component of the geostrophic velocity υg = −4.7 cm s−1, so the angle between the mean flow and the zonal direction is 239°. Figure 1 shows that the direction of the wave vector k, l associated with the fastest temporal growth Im(σ+(k, l)) for the root σ+(σ−) of (13) is 59° (239°) so that the mean flow direction determines the direction along which the fastest-growing Fourier component wavenumbers lie. The direction of the wavenumber vector of the fastest-growing modes is determined by the competition between maximal energy release for a wave vector oriented parallel to the mean flow and the β effect, which stabilizes the flow for the perturbations across the mean potential vorticity gradient (Pedlosky 1987). The fact that the direction of the wavenumber vector of the fastest-growing mode is almost parallel to the direction of the mean flow suggest that the vertical shear is strong enough to overcome the stabilizing effect of β. The Rossby radii are 41 and 17 km, values comparable to those in the related work of de Szoeke and Chelton [1999, (2.11)].
Figure 2 shows the Fourier summation solution for η1 and η2 at 500, 900, and 1200 days after the initiation of forcing. The solutions after 2000 days are shown in the upper panels of Fig. 4. The origin of coordinates x = 0, y = 0 in the figure is the center of the region of forcing. The background flow and forcing parameters are those of Fig. 1. After 900 days a well defined beam of stable waves is directed slightly toward the north of west away from the region of forcing, but the full Fourier solution contains an additional transient feature (17) and (19) that travels westward and clearly grows substantially from 900 days to 2000 days. Careful examination of Fig. 2 shows that the in-phase nature of the fluctuations of η1 and η2 that prevails in the beam of the sympathetic solution is replaced by an order 90° phase shift in the growing region. The growth is clearly the result of baroclinic instability of the background steady flow, already noted by Liu (1999a), Dewar and Huang (2001), and Kubokawa and Nagakura (2002).
Numerical evaluations (not shown) of the Fourier solution for the beam indicate that, as the frequency of Ekman pumping is increased, the distance between extrema in the beam decreases and the amplitude of the beam also decreases, but the direction of the beam is not visibly frequency dependent. For the present parameters, radiation in the beam is virtually extinguished for forcing frequencies of 0.5 cpy or higher. It is clear from the dispersion relation (13) that increasing the frequency decreases the wavelength of the sympathetic solution, but simple arguments do not predict the direction of the beam nor do they explain the decrease in beam intensity at high frequencies.
The Fourier summation and finite difference methods are very flexible but, on account of necessarily limited spatial resolution and “wrap around” in the periodic spatial domain, they are not well suited to answer a number of questions raised by the solutions of Fig. 2. Those questions are:
(i) Is the beam appearing in Fig. 2 the entire sympathetic solution?
(ii) Why is the beam so strongly directional?
(iii) Why is radiation at the forcing frequency strongly inhibited at higher frequencies?
(iv) What is the long-time structure of the transient solution?
(v) Is the instability of the transient solution convective in the manner defined by Huerre and Monkewitz (1990)?
To answer these questions it is necessary to develop solutions free of the shortcomings of the Fourier method. We therefore next work out analytical approximations to the solutions valid for long times after the initiation of forcing.
b. The sympathetic solution
The foregoing Fourier component summation solution of the initially forced problem consists of a transient, (17) and (19), and a sympathetic, (18) and (20), part. It was seen in Fig. 2 that this transient solution may become unstable. The fact that the plane wave dispersion relation (13) may have complex wavenumbers k and l for real frequencies σ suggests, as noted by Liu (1999a), that the sympathetic solution might grow spatially. Does it do so?
The full sympathetic solution, with allowance for possible spatial growth, may be found by a highly simplified version of the method of Briggs (1964, chapter 2: 8–46). If we rewrite the Fourier transform (10) of the Gaussian disturbance (5) as Ê(k, l) = Ê(k)Ê(l) with Ê(k) = LπW0e−k2L2/4, and so on, then the sympathetic solution (18) may be rewritten as
Here ηx2S(x, l) may be written in terms of a Green's function G2 as
with similar expressions for ηx1S(x, l) and G1. From (1), G1 and G2 obey the coupled equations
Here G2 is continuous across the origin x = 0 but G1 jumps discontinously in the amount −1/C as x increases across the origin: G1 and G2 are of the form Aeik+x + Beik−x away from the origin, where A and B are constants and k+ and k− are the solutions of the quadratic equation in (13) with real l and with frequency σ0 (the superscript signs correspond to the sign choice in the quadratic formula in solving for k). Once A and B have been chosen, then the integrations in (21) and (23) may be done numerically.
The solutions k+ and k− are plotted in Fig. 3 for real l and for the annual frequency σ0. They are complex when l lies between two critical meridional wavenumbers lup and llo (where lup > llo). We thus consider the three separate cases l > lup, l < llo, and lup > l > llo. For l > lup both k+ and k− are real, and may be shown to correspond to free waves with westward group velocity. The radiation condition thus requires that
The constants A and B are fixed by the continuity of G2 and the jump of G1 across the origin noted above. For l < llo both k+ and k− are again real and may be shown to correspond to free waves with eastward group velocity for k+ and westward group velocity for k−. The radiation condition thus requires that
Again the constants A and B are fixed by the continuity of G2 and the jump of G1 across the origin noted above. For the particular parameters of the solutions of this paper, very little energy appears in meridional wavenumbers l < llo on account of the small amplitude of the Gaussian forcing at meridional wavenumbers l < llo in wavenumber space.
For llo < l < lup both k+ and k− are complex, and the radiation condition no longer determines the form of the solution on either side of the origin. Pedlosky (2005, personal communication) has pointed out that, in this case, a heuristic version of the procedure of Briggs (1964) may be used to obtain the form of the solution away from the origin. The procedure is to add artificially large damping (−iσ0 → −iσ0 + ɛ, ɛ large) and to determine the form of the solution (the constants A and B) on either side of the origin by retaining only those terms that correspond to dissipatively dominated decay away from the origin. Once this has been done, let the damping decrease to zero without changing A or B. If a solution that decays away from the origin when dominated by large dissipation grows in the same direction when the dissipation is sufficiently reduced toward zero, then that solution is spatially unstable and must be retained. This is the frictional equivalent of the causality condition enforced by requiring the radiation condition in the case where k+ and k− were real. In the present case, when the substitution, −iσ0 → −iσ0 + ɛ, is made in the dispersion relation (13), then for large ɛ both k+ and k− correspond to westward decay but as ɛ → 0, k+ switches to westward growth. Consequently, for llo < l < lup,
and the constants A and B are fixed by the continuity of G2 and the jump of G1 across the origin noted above. The result is that the sympathetic solution is spatially unstable and grows westward.
Because meridional wavenumbers for which l < llo make effectively no contribution to the present solution, that solution may also be obtained by making a finite difference approximation to the y derivative in (1) and then integrating the resulting equations with harmonic time dependance westward from the eastern boundary. Parameters such as domain size, etc. enter the finite difference code in a simpler manner than the code that evaluates (21)–(27) so that the solutions for η1(x, y) and η2(x, y) shown in the lower panels of Fig. 4 for the parameters and domain of Fig. 2 were actually obtained by the finite difference method.
Figures 2 and 4 make clear that, as the full solution evolves, ever higher amplitudes are found ever farther west and south of the initially forced region. Comparison of the full solution (Fig. 4, upper panels) after 2000 days with the sympathetic solution (Fig. 4, lower panels) shows that, as the most rapidly growing region moves westward, it “leaves behind” the sympathetic solution, first the narrow stable beam (composed of Fourier components with meridional wavenumbers l > lup) and then the westward growing part (composed of Fourier components with meridional wavenumbers llo < l < lup). The latter is just beginning to be clearly visible by day 2000.
The sympathetic solution η2(x, y) given by (21)–(27) may be thought of as a Green's function convolved in x and y with the Gaussian forcing. In the course of solving (21)–(27) when l > lup, that Green's function (not shown) could be constructed by letting the factors Ê(l) and Ê(k) in (21) and (22) be constant; it was a complicated function of position characterized by very different local wavelengths at different locations. When convolved in x and y with the Gaussian forcing, this complicated pattern became the simple beam pattern of Fig. 4 because the very anisotropic dispersion relation (13) forces the longest waves to be found in only one particular direction, that of the beam, and it is only those largest-scale waves that survive the smoothing associated with the convolution. As the frequency is increased, the dispersion relation (13) indicates that the waves become shorter and are consequently more strongly attenuated by the convolution; that is why increasing the frequency ultimately extinguishes the beam.
c. Stationary phase approximation to the transient solution
To be able to answer the questions posed when analyzing the numerical solution at the end of section 1, in particular the question of whether the instability observed in the solution is absolute or convective, we next construct an analytical approximation to the transient solution.
The initial conditions satisfied by the transient solution (17) must cancel the sympathetic solution (18) evaluated at the initial time. But, the sympathetic solution evaluated at the initial time is a complicated and spatially extended function of x and y. To more clearly understand the evolution of the transient solution, we instead consider more compact initial conditions
in which both η10(x, y) and η20(x, y) are taken proportional to the Gaussian e−(x2+y2)/L2 of (5). The solution for η̂2T(k, l, t) is still of the form (15) in which the coefficients A and B are given by
Inverse Fourier transformation of (15) then furnishes η1T and η2T.
We solve the homogeneous initial value problem asymptotically for large time using the method of steepest descents. The classical asymptotic inverse Fourier transformation (8) for the initial value problem as t → ∞ is given by a sum of terms each one of the form
where σ0 = Ω(k0, l0) denotes the dispersion relation (13), ∇kΩ0 = (x/t, y/t) defines the stationary phase points (e.g., Carrier et al. 1966), and G(k, l) gathers together all the remaining k and l dependence of the integrand. It can easily be shown from (13) that
everywhere in wavenumber space and the form (30) fails in general. (This comes about because the long wave assumption has been made.) Nevertheless, there is an orientation of coordinates that makes the matrix
diagonal and one of the eigenvalues is nonzero; we thus may integrate in the direction of the nonzero eigenvalue. We accordingly work in cylindrical polar coordinates; the wavenumber vector components k and l are then described in terms of the magnitude of the wavenumber vector κ and the angle α that the wavenumber vector makes with the zonal direction. We first perform the integration over the angle α by the method of steepest descents. Thereafter, for the present special case of a Gaussian initial disturbance we may perform the integration over the wavenumber magnitude κ again by the method of steepest descents. The radial profile of the cylindrically symmetric initial disturbance thus influences only the result of the last integration.
The middle layer interface displacement η2T (x, y, t) is then a sum of terms of the form
in which σ is either σ+ or σ−. In polar coordinates κ, α this term becomes
where Φ is the phase given by
and we have introduced
with a = 1 + c/C. Here ξ and ζ are scaled horizontal coordinates that move westward with the frontal feature of the disturbance (defined below).
The α integration in (32) may be carried out using the method of steepest descents for large time t. The stationary phase points α0(ξ, ζ) are solutions of
in which ϑ ≡ tan(α). To solve this transcendental equation for α0(ξ, ζ), make use of the dispersion relation (13) written in terms of H,
Solutions of this are the quantity
evaluated at the stationary phase points α = α0. The value of α0 itself then readily follows from (36). The stationary phase values α0 may be complex, corresponding to growth or decay of their contribution to the solution. The term (32) thus becomes
where Φ0 denotes Φ(α0), ɛ(α0) = π/4 − arg(Φ0αα)/2, and the phase of Φ0αα is chosen so that −π/2 < Φ0αα < 3π/2. The κ integration is now readily carried out by stationary phase when t is large to finally obtain
in which κ0 is the stationary phase point
and in which the sum reminds us that there are generally several stationary phase points for each field point ξ, ζ.
Expression (40) shows that the maximum value of the solution, initially the peak of the initial Gaussian, remains very near the curve in the ξ, ζ plane defined by
for times sufficiently short that the solution is not dominated by growth associated with complex stationary phase values α0. We call this curve the front (Fig. 5). To obtain an equation for the front whose coefficients are not functions of α0(ξ, ζ), use (36) to eliminate Hα(α0) from (35) and combine with (42) in the form ξ + ζϑ0 − H0 = 0 to find ζ = υ/(2H0). Now use this and (36) to eliminate reference to ϑ0 and to H0 from (42) to finally obtain the equation for the front as
in which υ is defined in (36). For large |ξ|, the front asymptotes to ξ = ζ b/υ and to ζ = 0. It has focii at ζ = ±−υ2b/4, ξ = 0. Since the coordinates ξ and ζ move westward at speed (c + C)/2, the frontal pattern of Fig. 5 also propagates westward without change at this speed, even though the full solution changes form as it evolves in the moving coordinates.
The solution (40) displays caustics, shown in Fig. 5. They represent a boundary between a region with a complicated wave pattern where groups of waves interfere and a neighboring region where two fewer groups of waves (possibly none) interfere (Lighthill 1978). At a caustic the first and the second derivatives of the phase vanish:
causing the stationary phase solutions to break down at the caustics. In principle this local difficulty may be overcome by using the Airy integral to produce a “healed” version (Lighthill 1978). That has, however, not been done for the solutions discussed below (Fig. 5).
The first of the equations in (44) may be solved to find the stationary phase values α0 of α as in the previous section yielding (35). For each stationary phase value of α0, (44) constitute two equations for ξ and ζ whose solutions define the caustic curves. Insert ξ obtained from the first of the expressions (44) into the second one to obtain ζ as
The caustics, asymptote at large distances ξ, ζ to π/2 and arctan(−b/υ). The caustics meet in cusps at points of the ξ, ζ plane where the first three derivatives of phase with respect to the wavenumber vector direction α are zero:
Between the caustic curves of Fig. 5 the parts of the solution corresponding to the complex stationary phase points grow in time for the present baroclinically unstable case. In the remainder of the domain however, outside of the caustic curves, all the stationary phase points the contribute to the solution are real and do not result in growth. In this manner, the confinement of complex stationary phase points to just a part of the ξ, η plane causes the growing part of the solution to be restricted to just a part of the x, y plane at any time so that the instability is convective in the manner defined by Huerre and Monkewitz (1990) rather than absolute.
3) 3) Discussion of stationary phase solution
The initial conditions are of the form (28) with the particular choice η20(x, y) = e−(x2+y2)/L2 and η10(x, y) = 0. This choice insures that both fast and slow modes are strongly energized in the solution. The choice of stationary phase points is discussed in appendix B. As noted above, the Gaussian scale is chosen to be L = 500 km, whereas the Rossby radii are 41 and 17 km.
The stationary phase solutions (40) for η2 are displayed in Fig. 5 in the coordinates ξ, ζ at two different times. A point fixed in those coordinates moves westward with speed (c + C)/2 and radially outward at speed C in the physical space coordinates x, y. Consequently much (but not all) of the evolution of the solution in physical space coordinates x, y is just westward drift at speed (c + C)/2 and radial expansion at speed C. These stationary phase solutions are in good agreement with their less well resolved Fourier counterparts except in the vicinity of the caustics. We may now answer a number of questions about the transient solution. 1) Why does the dominant part of the transient solution quickly elongate itself north-of-west to south-of-east? 2) Why does the solution evolve from the initial Gaussian to a pattern with two distinct maxima (Fig. 5a) to a pattern with a single central maximum (Fig. 5b)?
In answer to question 1, initially the distribution of energy in wavenumber space is isotropic and so is the solution. As Fourier components whose wavenumbers lie in the unstable wedges of Fig. 1 grow, however, the solution increasingly consists of waves whose crests and troughs are oriented north-of-west to south-of-east at right angles to the axis of the wedge of unstable wavenumbers of Fig. 1 and its counterpart for σ− (not shown). For question 2, as the solution initially evolves (Fig. 5a) the largest amplitude for η2 occurs approximately at that zonal position in the basin where the slow mode contribution to the solution would have been centered in the absence of background flow, but there is a strong secondary maximum farther to the west at which the fast mode contribution to the solution would have been centered in the absence of background flow. The same two extrema appear in η1 (not shown) but with different relative amplitudes. Kubokawa and Nagakura (2002) obtained a similar pattern by summing Fourier modes; whereas they discuss the solution in different parts of the ocean, we focus on a more complete description of the space and time variation of the solution. The maximum amplitude of η2 for the solution of Fig. 5a is about 1.1m whereas the maximum initial amplitude of η2 was 1.0 m; it is thus clear that baroclinic instability is not yet an important factor in the evolution of the solution.
At some time, when solving by summing Fourier modes, the smallest-scale resolved components are energized. At that point they can no longer interfere with yet-smaller-scale components in the same manner as they do in the fully resolved solution so that the shortest resolved components spuriously dominate the numerical solution, giving it a “checkerboard” texture. The stationary phase solution does not have such a short wave resolution limit. Figure 5b thus shows the stationary phase solution at a time appreciably later than that of Fig. 5a. Now the maximum amplitude is significantly greater than that of the initial condition, indicating that baroclinic instability has become a more important factor in the evolution of the solution.
By whatever method they have been obtained, the solutions remain valid providing that the dominant horizontal scale are larger than the baroclinic Rossby radii. The limiting time scale at which this ceases to be true has not yet been reached in the solutions of Fig. 5.
What is the nature of the different wave groups that interfere to produce the full transient solution? The stationary phase angles α0(ξ, ζ) may be deduced from the quantity H defined in (35). Figure 6 shows H obtained by solving the quartic (38) for ξ, ζ lying on two circles (of radius R = 10 and R = 1.5) centered at the origin of coordinates in the ξ, ζ plane as a function of the angle μ = tan−1(ζ/ξ). The four angles associated with the two caustics shown in Fig. 5 are clearly visible. Also shown are large R approximations to the locations of the caustics (appendix C).
At large R (Fig. 6a), there are four real stationary phase points α(ξ, ζ) for the largest values of μ. As μ decreases through successive caustic angles, one pair of these coalesces into a complex conjugate pair at the first (largest) caustic angle, that pair divides into two real values at the second caustic angle, the other real pair coalesces into a complex conjugate pair at the third caustic angle, and that pair finally divides into two real values at the fourth (smallest) caustic angle. In this way, the caustics divide the ζ, ξ plane into two sectors (roughly vertically above and below the dominant region of the transient solution) in which some local wavenumber components grow exponentially, and two sectors (roughly east and west of the dominant part of the transient solution) in which all local wavenumber components propagate without growth.
At small R (Fig. 6b), the range of μ for which all stationary phase points are real is much smaller than for larger R, and correspondingly the region of the ζ, ξ plane in which all local wavenumber components propagate without growth is much smaller for radii near the central region of the solution than for radii further away.
Is the transient solution absolutely or convectively unstable? This question cannot be definitively answered using the method of Fourier summation because both limited resolution and “wrap around” make it difficult to examine the long time behavior of the solution at a fixed point. With the stationary phase solution, however, it is relatively easy to evaluate, for example, η2(x, y, t) either at fixed x, y as t increases or else at fixed x/t, y/t as t increases. Figure 7 shows the result at a number of different times.
The growth following the disturbance is clearly somewhat more rapid than a simple exponential in time. This reflects the successive instantaneous dominance of ever shorter, and hence (Fig. 1) more unstable, Fourier components. As noted above, the present solutions will no longer be valid when Fourier components whose scales approach the baroclinic Rossby radii are excited. That limit has not been reached in any of the solutions presented here but it ultimately limits their validity.
Since the solution at the origin ultimately tends toward zero even though the solution following the disturbance clearly grows without limit, we may say that the PG solution is convectively unstable, but the PG approximation precludes literally taking t → ∞.
5. Discussion and conclusions
The purpose of this paper is to study the nature of long baroclinic wave propagation and instability in an idealized unventilated subtropical gyre flow with vertical shear. The new results are a detailed analytical description of the spatial structure of the long planetary wave disturbance that results when the flow in the eastern part of an unventilated steady subtropical gyre is perturbed by a region of spatially localized Ekman pumping initiated at some initial time and thereafter oscillating at a frequency on the order of a cycle per year. Both the background flow and the perturbations are solutions of a 2½-layer model, the region of forcing is located in the eastern part of the gyre where the steady flow is confined to the uppermost layer.
Besides linearization and restriction to unventilated background flows, the most important approximations are (i) neglect of relative vorticity so that only long (PG) planetary waves are considered and (ii) neglect of the spatial variation of the background flow so that phenomena such as refraction are neglected. The first approximation is not violated provided that the horizontal scale of the forcing is initially chosen to be significantly larger than the baroclinic Rossby radii characterizing the background flow. Chelton et al. (2004) showed that persistent wind features exist on scales of several hundred kilometers. The solutions studied in this paper are all for times shorter than the time at which the horizontal scales of the disturbance have decreased to the Rossby radii characterizing the background flow. The second approximation requires that the horizontal scales of the gyre be substantially greater than the horizontal scale of the forcing. Both approximations are reasonably well satisfied for long waves propagating through a background flow with weak spatial variations over sufficiently short times.
The most important results of the analysis are as follows.
The result of spatially localized surface Ekman pumping initiated at some initial time and thereafter oscillating at a frequency on the order of a cycle per year is (i) a directly forced part of the solution whose stable wavenumber contributions form a well-defined and narrow beam of stable waves emanating from the forcing region (Fig. 4), co-oscillating at the forcing frequency and not growing spatially plus an unstable contribution that grows spatially toward the southwest, and (ii) an accompanying transient (Figs. 2 and 5) that emanates from the forcing region with a westward component of velocity and grows in time as it propagates.
The direction of the co-oscillating beam is fixed (to be slightly north of westward) by the dispersion relation for free waves and by the degree of spatial localization of the forcing region (forcing regions of large lateral extent radiate much longer waves than those of small lateral extent). For the particular parameters of the examples of this paper, a Gaussian forcing region of lateral scale 500 km and frequency 1 cpy radiates strongly but, if the frequency is increased to several cycles per year the wavelengths allowed by the dispersion relation are not efficiently radiated from the forcing, and the radiation is strongly attenuated.
The structure of the transient is somewhat similar to that found by Kubokawa and Nagakura (2002) in the unventilated part of their subtropical gyre, but the method of solution (stationary phase) yields a more complete description of the space and time variation of the solution—for example identifying the frontal feature (43) and the caustics that propagate westward without change of form even though the form of the full solution evolves as it propagates westward—and also (always within the limitations of the PG approximation) allows a clear determination that the transient is convectively—but not globally—unstable.
Spall (1994) made use of QG theory in a case in which the mean flow was only a function of depth to study possibly unstable perturbations of the background state representative of the real ocean, and found the fastest growing perturbations to be at scales not greatly larger than the largest baroclinic Rossby radius O(100 km). The purpose of the present study is to provide detailed analysis of baroclinic instability acting on the larger scales at which the planetary geostrophic approximation is valid. To generalize these results to more realistic ocean subtropical gyre flows, it would be necessary to allow for the presence of smaller scales and horizontal variation of the background flow. Further analysis is thus required to determine how the horizontal spatial variation of the background flow, not allowed for in this paper, affects these solutions. It would be of particular interest to determine whether the baroclinic instability mechanism identified in this paper at larger PG scales would still play a significant role when dynamics important at scales of the order of baroclinic Rossby radius of deformation are included.
Thierry Huck and an anonymous reviewer suggested important changes in the manuscript. Author IC especially thanks Joe Pedlosky for valuable suggestions and several motivating discussions, as well as Ronald Guenther and Roger Samelson for their helpful comments and continuous interest in my work. This work was supported by the National Science Foundation through Grants OCE-9907008 and 0220471.
The Dynamical Model
The steady flow is modeled as in Luyten et al. (1983). Three layers (upper, middle, and lower, labeled 1, 2, and 3) are considered. The depth of the lower layer will ultimately be allowed to become very large. The momentum balance is geostrophic,
where subscripts x and y denote differentiation, uj and υj are horizontal fluid velocities in the zonal (x) and meridional (y) directions within layer j, and the Coriolis parameter f is linear in y, f (y) = f0 + βy. Pj = (pj + ρjgz)/ρ0 is the Montgomery function, whose gradient depends only on the lateral position (x, y) within each layer; pj is the true dynamic pressure in j layer, ρj is the (constant) density in the layer, and ρ0 is a reference density. The vertical coordinate z is positive upward.
If Hj is the depth of the base of layer j, then the hydrostatic assumption causes the Montgomery function in layer j to be related to that in layer j + 1 by
in which γj = g(ρj+1 − ρj)/ρ0 are reduced-gravity parameters. In terms of the layer thickness hj = Hj − Hj−1, the mass conservation equation for layer j is
in which subscript t denotes the time derivative, J denotes the Jacobian J(a, b) = −aybx + axby, and wE is the wind-induced vertical velocity at the base of the surface Ekman layer. The sum of (A3) with constant H3 gives the Sverdrup transport relation
The steady flow occupies a midlatitude basin of zonal extent a and meridional extent Δ with meridional boundaries at y = b − Δ and y = b and is driven by Ekman pumping of amplitude W0 having the form
in which f0 is the Coriolis parameter at the central latitude of the basin. The steady Ekman pumping is chosen to be sufficiently weak that the deeper layers of the model never surface so that the steady flow occupies only the uppermost layer.
We study time-dependent perturbations linearized about the steady flow that are excited by harmonic forcing turned on at some initial time or by initial perturbations. We separate the variables into steady part and time-dependent part as
The perturbation mass conservation equations (A3) for the middle and the lower layer become
where H02 = H20 is the depth of the bottom of the middle layer (constant when the steady flow in that layer is stagnant) and H30 is the flat bottom of the basin. The perturbation Sverdrup equation, obtained by integrating (A4) westward from the eastern boundary x = a and linearizing about the steady flow, is
Stationary Phase Points
To determine the set of stationary phase points in the complex α plane over which to sum in (40), it is necessary to deform the integration path α = (0, 2π) of (32) into an integration path that passes through each stationary phase point along that point's path of steepest descent. This must be done separately for the integrands involving σ+ and those involving σ−.
The only singularities of the integrands of the terms (32) in the complex α plane are at the simple zeros of H2 (36), whose square root appears in the phase (33). These zeros occur in real pairs, the location of the second pair being the location of the first pair plus π radians. In deforming the α contour of integration the two branch cuts joining the members of each pair separately must not be crossed. It is convenient to choose the branch cuts to lie along the real axis between the two members of each individual pair.
Once H has been determined, solution of (36) for the stationary phase values of α yields four stationary phase points whose real parts lie between (0, π) and another four stationary phase points whose real parts in (π, 2π) are the real parts of the first four plus π. For every physical space point (ξ, ζ), the interval (0, π) contains three stationary phase points (either three real points or one real point and one complex conjugate pair) associated with σ+ and one (real) point associated with σ− or vice versa. A similar statement holds for the stationary phase points whose real parts lie in (π, 2π).
When all stationary phase points are real, no deformation of the integration path α = (0, 2π) is needed, and the sum in (40) must include contributions from all stationary phase points. When the integrand of (32) under consideration has one real stationary phase point and a complex conjugate pair whose real parts is in (0, π) with corresponding points in (π, 2π), it is not possible to deform the contour α = (0, 2π) to lie entirely along the steepest descent paths threading one or more of the stationary phase points without crossing the branch cuts joining the simple zeros of H2 (36). If however the deformed contour is made to parallel the branch cuts where it meets them, then along those cuts the imaginary part of the phase in the integrand is positive and large so that as t → ∞ the contributions to the integral from these portions of the deformed path will be exponentially smaller than those from the stationary phase points. Because the complex stationary phase points come in complex conjugate pairs, the sum over all stationary phase points threaded by the deformed contour is the same as the sum over all the real stationary phase points plus one-half times the sum over all the complex stationary phase points. This then is the general rule for each integrand of (32).
Approximate Stationary Phase Solutions Valid at Large Distances away from the Center of the Disturbance
We seek approximate analytic solutions, valid at large distances away from the center of the disturbance, of the quartic equation (38) for H as defined in (34), which yield the four values of the stationary phase angle α0.
Introduce the nondimensional variables as
and transform to cylindrical polar coordinates to write
We seek approximate solutions of this equation for large R.
First seek O(1) solution by neglecting much smaller terms of the order R−1; one iteration on the two roots thus obtained yields
From the dispersion relation (36)
and the two values of α0 are given by
They are two O(1) solutions for which α0 is close to μ, but the correction term in (C5) diverges at μ = π/2 and μ = arctan(−b/υ).
An O(R) solution is obtained from the two highest-order terms in (C3), which are O(R3), and give H4 − 2H3R cosμ = 0, whose solution is H = 2R cosμ. In terms of the angle μ0 defined by
the once iterated O(R) solution becomes
and yields α0:
The correction of this root blows up at μ = π/2.
An O(R−1) solution is approximated by
It becomes, after one iterative correction,
and yields α0
This correction blows up for μ = μ0. The four roots of α0 given by (C6), (C9), and (C12) are the approximations to the angle which makes the phase stationary and over which (40) must be summed. When they are inserted in (40), the result shows that the asymptotic solutions for large R decay transcendentally away from the center of the front. Solutions (C6), (C9), and (C12) are crude approximations, and the correction terms diverge at μ = π/2 and μ = μ0 = arctan(−b/υ). The more accurate expressions for the position of the caustics at large distances from the center of the disturbance, as well as approximate solutions for large distances that are uniformly valid across the caustics, can be obtained by keeping all the terms of the quartic equation in (38) that interact at the caustic considered.
* Current affiliation: Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, Massachusetts
Corresponding author address: Ivana Cerovečki, 54-1721 Dept. of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139-4307. Email: firstname.lastname@example.org