Parametric Subharmonic Instability and the β-effect by Ian Chan B.Sc., University of British Columbia, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Mathematics) The University of British Columbia (Vancouver) August 2010 c© Ian Chan, 2010 Abstract Parametric subharmonic instability (PSI) is a nonlinear interaction between a res- onant triad of waves, in which energy is transferred from low wavenumber, high frequency modes to high wavenumber, low frequency modes. In the ocean, PSI is thought to be one of the mechanisms responsible for transferring energy from M2 internal tides (internal gravity waves with diurnal tidal frequency) to near-inertial waves (internal gravity waves with frequency equal to the local Coriolis frequency) near the latitude of 28.9◦. Due to their small vertical scale, near-inertial waves generate large vertical shear and are much more efficient than M2 internal tides at generating shear instability needed for vertical mixing, which is required to maintain ocean stratification. The earlier estimate of the time-scale for the instability is an order of mag- nitude larger than the time-scale observed in a recent numerical simulation (MacK- innon & Winters, 2005) (MW). An analytical model was developed by Young et al. (2008) (YTB), and their findings agreed with the MW estimation; however as YTB assumed a constant Coriolis force, the model cannot explain the intensificaiton of PSI near 28.9◦ as observed in the model of MW; in addition, the near-ineartial waves can propagate a significant distance away from the latitude of 28.9◦. This thesis extends the YTB model by allowing for a linearly varying Coriolis parameter (β-effect) as well as eddy diffusion. A linear stability analysis shows that the near-inertial wave field is unstable to perturbations. We show that the β-effect results in a shortening in wave length as the near-inertial waves propagate south; horizontal eddy diffusion is therefore enhanced to the south, and limits the meridional extent of PSI. The horizontal diffusion also affects the growth rate of the instability. A surprising result is that as the horizontal diffusion vanishes, the system becomes stable; this can be demonstrated both analytically and numerically. Mathematically, the β-effect renders the spatial differential operator nonnor- mal, which is characterized with the aid of pseudo-spectra. Our results suggest the possibility of large amplitude transient growth in near-inertial waves in regimes that are asymptotically stable to perturbations. ii Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 IGWs and the dynamics of the ocean . . . . . . . . . . . . . . . . . . 3 1.2 Parametric subharmonic instability . . . . . . . . . . . . . . . . . . . 8 1.3 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Mathematic Description of PSI . . . . . . . . . . . . . . . . . . . . . 12 2.1 Derivation of the PSI model . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.1 The pumpwave train . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.2 The simplified amplitude equation . . . . . . . . . . . . . . . 18 2.2 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Dissipative Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 An initial value problem . . . . . . . . . . . . . . . . . . . . . . . . . 23 iii 3.1.1 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Normal modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.1 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.2 Analytical results . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2.3 Energetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 More initial value problems . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.1 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.2 Boundary modes . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4 PSI in the Inviscid Limit . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1 Stability of the inviscid problem . . . . . . . . . . . . . . . . . . . . 60 4.2 Non-normal operators and transient growth . . . . . . . . . . . . . . 62 4.2.1 A simple example . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2.2 Non-normal operators and fluid dynamics . . . . . . . . . . . 66 4.3 Pseudospectra and PSI . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.1 Pseudospectra . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.2 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3.3 Numerical difficulties . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.1 Horizontal dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.3 Outstanding issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 iv List of Figures 3.1 An initial value problem . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Effect of dissipation on the spectrum: l = 0 . . . . . . . . . . . . . . 30 3.3 Effect of dissipation on the spectrum: l = 0.1 . . . . . . . . . . . . . 31 3.4 Effect of computational domain on the eigenvalues . . . . . . . . . . 33 3.5 Effect of the parameters on the growth rate. . . . . . . . . . . . . . . 36 3.6 Energetics of PSI: l = 0 . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.7 Energetics of PSI: l = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . 49 3.8 Energetics of PSI near the stability boundary . . . . . . . . . . . . . 52 3.9 An initial value problem with transient growth . . . . . . . . . . . . 55 3.10 An initial value problem with stable normal modes . . . . . . . . . . 56 3.11 An initial value problem with boundary modes . . . . . . . . . . . . 57 4.1 Transient growth for a simple dynamical system . . . . . . . . . . . . 65 4.2 Pseudospectra for a stable system . . . . . . . . . . . . . . . . . . . . 69 4.3 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4 Initial condition and transient growth . . . . . . . . . . . . . . . . . 72 4.5 Effect of non-normality on calculation of eigenvalues . . . . . . . . . 75 5.1 Pseudospectra for a stable system . . . . . . . . . . . . . . . . . . . . 81 v Acknowledgements First and foremost, I would like to thank my supervisor Neil Balmforth for being an inspiring mentor over the last three years. It has been a pleasure working with him and I have learnt a great deal from him, which will no doubt be invaluable to my career. I would also like to thank Kraig Winters and Bill Young for their comments for this manuscript. Navigating myself through the degree would’ve been much more difficult without the help from the faculty members at UBC. I would like to thank Douw Steyn, Philip Loewen, Eric Cytrynbaum, Mark Maclean, George Bluman, and Richard Anstee for their constant guidance, which helped me grow immensely in my research and teaching. I would also like to thank the staff of the Department of Mathematics and the Institude of Applied Mathematics, especially Lee Yupitan for answering the countless administrative questions I had and Marek Labecki for his help on solving computer related issues. Finally I would like to thank my family for their support during my studies. I would like to thank my wife Giorgia especially for her constant encouragements and love. The research is partially supported by the NSERC PGS-M scholarship. vi Chapter 1 Introduction The ocean is constantly covered by waves, which can range from small ripples gen- erated by a gentle breeze, to giant rogue waves found in storms. These waves exist due to gravity: if the water surface is perturbed upwards, gravity will act on the water to pull the surface downwards to restore the surface back to its equilibrium; however most often it will overshoot, causing a depression to form, and gravity will now act on the water surrounding the depression to fill in the depression. This sequence will occur repeatedly, which results in the familiar waves we observe on the ocean surface. Gravity waves are not confined to surface of the ocean, as they also exist underneath the surface. Due to differences in temperature and salinity, the density in the ocean typically increases with depth, and therefore when a parcel of water is perturbed downwards, the denser ambient fluid provides a positive buoyancy that acts as a restoration force. The resulting waves are called internal gravity waves (IGWs), and these waves are characterized by the up and down motion of isopycnal (constant density) surfaces. These internal gravity waves have several features that differ from their coun- 1 terparts on the surface. Firstly, the sharp density contrast between air and water results in a much higher frequency of oscillation for the surface waves; the typical time scale for surface waves observed near a beach is measured in seconds. In con- trast, the restoration force for IGWs is provided by a more gentle density gradient, and the frequency of vertical oscillations can be characterized by the Brunt-Väisälä frequency N : N = √ − g ρ0 dρ dz , (1.1) where g is the gravitational acceleration, ρ0 is the characteristic density, and z is the vertical coordinate. Below the ‘mixed layer’ on the surface of the ocean, which is uniform in density and extends up to 50m, the buoyancy frequency decreases from 10−2s−1 to about 10−4s−1 in the deep ocean (about 4000m) (Thorpe, 2005); these frequencies correspond to time-scales ranging from minutes to hours. The IGWs also have a much larger spatial scale compared to surface gravity waves. For example the internal waves associated with the diurnal tide can span the entire depth of the ocean (103m), while having horizontal scales measured in tens of kilometres, which is in stark contrast to surface waves that typically have heights of 0.1 − 10m and horizontal wavelengths of perhaps a tens of meters. A consequence of the large horizontal scale of IGWs is that IGWs are susceptible to the Coriolis effect, which is reflected in the dispersion relation: k2 + l2 m2 = ω2 − f2 N2 − ω2 (1.2) which stipulates a relationship between (k, l,m) ≡ ~k (the wavenumbers in the x, y, and z direction respectively1), and the frequency of the wave ω. f denotes the Coriolis parameter, which depends on the earth’s rate of rotation (Ω = 7.27 × 1Throughout this thesis we will adhere to the convention in oceanography by orienting the positive x−axis towards east, positive y−axis towards north, and positive z− upwards 2 10−5s−1) and latitude (φ): f = 2 Ω sinφ. (1.3) Since the wavenumbers are real, and N f in general, an implication of (1.2) is that an IGWs can only exist if N ≥ ω ≥ f . In addition one can use (1.2) to show the direction of the group velocity makes an angle of α with the horizontal plane, where α satisfies: tanα = √ k2 + l2 m2 = √ ω2 − f2 N2 − ω2 . (1.4) Since the energy propagates along the direction of the group velocity, this shows that in general IGWs transports energy both in the horizontal and vertical direction, which is in contrast to the surface waves where energy is propagated along the surface. Two limiting cases can be considered: 1. ω → N : in this case the angle of propagation tends to 90◦, which corresponds to vertical oscillations with frequency N . 2. ω → f (The near-inertial limit): in this limit the angle of propagation tends to 0, and energy is propagated horizontally. In addition, as the RHS of (1.2) vanishes as ω → f , m2 must be large relative to k2+l2 in order for the equation to hold. Since the wave number are inversely proportional to the spatial scale, this implies that the vertical scale m−1 must become small compared to the horizontal scale (k2 + l2)−1/2. 1.1 IGWs and the dynamics of the ocean As warm water has a lower density than cold water, the constant heating at the surface will in principle result in a warm layer near the surface, and the rest of the ocean will consist of cold dense water underneath. In reality the density profile is 3 much smoother due to mixing, and one can find circulations that span the entire depth of the ocean (the thermo-haline circulation for example ). Molecular diffusion itself cannot account for the mixing, as the diffusivity of heat and salt are in the order of 10−6 − 10−7m2s−1, while the diffusivity needed to maintain the stratification is estimated to be 10−4m2s−1 (Munk & Wunsch, 1998). It is believed that the breaking of internal waves is responsible for the turbulence, and the diapycnal (across density surfaces) mixing that follows (Laurent & Garrett, 2002). Munk & Wunsch (1998) estimated that 2.1TW of energy is needed to main- tain the stratification in the abyssal ocean via mixing; as an additional 2.6TW of energy is dissipated through turbulence mixing in the shallow oceans, a total energy source of 4.7TW needed. Wind induced mixing, either directly near the surface, or indirectly via the excitation of geostrophic motions and internal gravity waves, pro- vides 1.2TW of energy. The remaining 3.5TW is supplied by the various modes of tides due to the sun and the moon; amongst the modes, the semi-diurnal component of the lunar tide (denoted by M2) is the dominate source, accounting for 2.5TW. The effect of the moon’s gravity field on the earth is not uniform everywhere on the surface, and the minute difference is sufficient to cause a deformation in the ocean: the result is that the ocean bulges out along the earth-moon axis, both towards and away from the moon. As the earth rotates around its axis, a point on the surface over the ocean will experience a rise and a subsequent fall in the sea-level as it passes through a bulge. As most locations will pass through both bulges, the result is the semi-diurnal tide that occurs approximately twice daily, or every 12 hours and 25 minutes. Since the earth’s rotation rate (7.292×10−5s−1) is much higher than the rate at which the moon orbits around the earth (2.662×10−6s−1), we can view the earth 4 as a sphere moving through a viscous film of fluid held in place by the gravity of the moon. The friction between the earth and the ocean causes the earth’s rotational rate to slow down, as well as a misalignment of the tidal bulge along the moon- earth axis. The misalignment in turn imparts a torque on the moon, causing it to gain rotational momentum and increase its distance from the earth. The rotational energy is partially dissipated in the ocean, and partially transferred into the orbital energy of the moon, causing it to move further away from the earth. The modern value of lunar recession, obtained via laser ranging, is determined to be 3.82 cm/year (Dickey et al., 1994), which corresponds to a tidal dissipation rate at 2.5 TW (Munk & Wunsch, 1998). Tidal dissipation in the ocean The subject of tidal dissipation can be dated back to Sir G.I. Taylor (1920) and Sir Harold Jeffreys (1921). Taylor estimated the the frictional dissipation over the Irish sea alone to be 0.051 TW, while Jeffreys estimated the global dissipation over shallow seas to be 1.1 TW; later Miller (1966) pegged the value at 1.6TW, while recent values inferred by satellite data is for the M2 mode is 1.8TW (Egbert, 1997; Munk & Wunsch, 1998). The discrepancy between the estimated dissipation over shallow seas and the total dissipation calculated from astronomical data can be attributed to the excitation of internal gravity waves due to bottom topography in the deep ocean. While the conversion over continental shelves is long known, the conversion over deep ocean has been thought to be insignificant until recently (Garrett, 2003). Part of the reason for the revision is due to the availability of satellite altimetry data, which suggests that dissipation over the deep ocean is responsible for converting 1TW of energy (Egbert & Ray, 2000, 2001). 5 The mechanism behind the tidal conversion is simple: as the tide moves around the earth, the water in the ocean is moved back and forth; in the presence of submarine ridges, the density surfaces are perturbed and internal gravity waves are launched in a fashion similar to lee waves generated by mountains. These internal gravity waves (called internal tides) then propagate from the bottom up into the open ocean, carrying energy upwards that ultimately leads to overturning and mix- ing (Garrett, 2003). These internal tides can be observed with the aid of satellite altimetry: for example, Ray & Mitchum (1997) reported internal tides propagating away up to 500km away from the Hawaiian Islands, primary as mode-1 baroclinic tide. The baroclinic tides often have the same frequency as the underlying barotropic tide due to the weak background tidal current (Garrett & Kunze, 2007), while the horizontal structure is around 20-100km (Laurent & Garrett, 2002). The Richard- son number for these internal tides exceeds unity, suggesting that the internal tide cannot lead to overturning directly; it is believed that energy is cascaded down via non-linear wave-wave interaction into smaller scale internal waves, which are more effective at causing over-turnings (Laurent & Garrett, 2002).. Non-linear wave interactions in internal gravity waves are thought to be responsible for shaping and maintaining the energy spectra of internal waves in the open ocean into the so-called Garrett-Munk (GM) spectrum (Müller et al., 1986; Thorpe, 2005). Garrett & Munk (1972, 1975) proposed a model for the observed energy spectra of internal waves that is in remarkable agreement with experimental data from different sites; this suggests that there is a single universal mechanism responsible for cascading energy from large scale motion to small scale mixing. As abyssal mixing is an essential component to the dynamics of large scale ocean circulations that in turns play a pivotal role in our climate, the GM spectrum and 6 ocean mixing have received plenty of attention over the last forty years. Non-linear wave interactions have been explored theoretically since the 1960’s. Phillips (1960, 1961) considered two waves with frequencies and wave-numbers (ω1,~k1) and (ω2,~k2), and considered a third mode whose frequency and wave-number satisfy the resonance condition: ω1 ± ω2 = ω3 , (1.5) ~k1 ± ~k2 = ~k3 . (1.6) In this case, the third wave will be forced resonantly, and energy is transfered from the first two modes into the third. Bretherton (1964) considered the interaction between a triad of discrete modes, where ω1 + ω2 = ω3 and ~k1 + ~k2 = ~k3, and provided a heuristic argument for the evolution equation of the amplitudes ai: iω1 ∂a1 ∂t = D̂a∗2a3 , (1.7) iω2 ∂a2 ∂t = D̂a∗1a3 , (1.8) iω3 ∂a3 ∂t = D̂a1a2 , (1.9) where D̂ is a real coefficient that describes the interaction, and generally depends on the wave-numbers. Hasselmann (1966, 1967) studied weak non-linear wave interac- tion for a continuous field of waves that is assumed to be Gaussian and statistically stationary (‘the random phase hypothesis’). McComas & Bretherton (1977) used a similar method to investigate the resonance interaction in the GM spectrum; they numerically computed the energy transfer in the spectrum, and concluded that the process is dominated by three classes of interactions: induced diffusion, elastic scat- tering, and parametric subharmonic instability. For Induced diffusion, wave action (wave energy divided by frequency) diffuses from a high-frequency, high-wavenumber 7 mode to one that has similar frequency and wavenumber. Elastic scattering occurs between two high frequency modes, with similar horizontal wave numbers but oppo- site vertical wave number; in this case there is rapid energy transfer between the two high frequency modes leaving the third low frequency mode largely unchanged. The last class of interaction is parametric subharmonic instability (PSI), where a large scale wave decays into two small scale waves with half the frequency. An important observation is that among the three pathways PSI is the only mechanism where energy is transferred into waves with a much different wave number. For a thorough discussion, please see McComas & Bretherton (1977) or Müller et al. (1986). 1.2 Parametric subharmonic instability Parametric subharmonic instability allows a wave to spontaneously excite daughter waves with higher wave numbers with half the frequency; it is a special case of (1.5) where ω1 = ω2 = ω3/2. Since the semi-diurnal internal tide has frequency angular frequency M2 = 2pi/12.42hr = 1.405 × 10−4s−1, the daughter waves will have frequency ω = M2/2 = 7.026×10−5s−1. As mentioned earlier, internal gravity waves can exist only if the frequency is greater than the local Coriolis frequency; thus for the daughter waves to exist, the latitude must satisfy the inequalityM2/2 ≥ 2Ω sinφ, or φ ≤ 28.9◦. The implication is that PSI can only occur between the equator up to the critical latitude 28.9◦. At the critical latitude the daughter waves have frequency that is equal to the local Coriolis frequency, and therefore they are near-inertial oscillations (NIO). While PSI have received some attention over the 1970’s, the interest waned as it was believed that the energy transfer via PSI occurs over a time-scale of 100 days for the mode-1 internal tide, which is too slow compared to other processes to 8 remove energy efficiently (Olbers & Pomphrey, 1981); however the interest in PSI have been renewed in recent years due to several new studies. Hibiya et al. (1998, 2002) carried out a set of numerical experiments, where a quasi steady internal wave field was perturbed by injecting energy into the first vertical mode at the M2 freqency. For the experiment at 49◦N, the spectra did not differ significantly over 10 inertial period; on the other hand, the energy was quickly transfered to the high wave-number modes in the f < ω < 2f range within 10 inertial periods at 28◦N. MacKinnon & Winters (2005) (henceforth MW) solved the full 3-D equation of motions numerically with a linearly stratified background; a mode-1 internal tide was excited at the southern boundary, which propagated northward in the presence of background noise. After the initial spin-up, instabilities started to develop near the critical latitude. The instability had much finer vertical structure compared to the mode-1 internal tide, and a plot of temporal spectra plot reveals a high spectral content over the M2/2 frequency. At the same time the energy flux associated with internal tide decreased by 62.5% as it crossed the critical latitude. In addition the time scale of the energy transfer was estimated to be 10 days, which was an order of magnitude smaller than the estimation of Olbers & Pomphrey (1981). MW attributed the disparity to the fact that the estimation of Olbers & Pomphrey was based on the random phase assumption, but from the numerical simulation the interaction were clearly much more coherent. The efficiency of PSI led the authors to describe this as the ‘subtropical catastrophe’. Simmons (2008) ran a global ocean model with realistic geometry, and observed energy transfer into the M2/2 modes near the critical latitude as well. A longitudinal cross section near a internal tide generation site showed that PSI intensified near the surface, possibly due to the non-linear stratification. 9 Recent field measurements also pointed to the presence of PSI activity in the ocean. For example, Rainville & Pinkel (2006) compared the vertical profiles at a site near Hawaii with strong internal tide activity (nearfield) to a site that is further away, and found a strong link between the semi-diurnal energy flux at the nearfield site and the diurnal energy flux at the farfield site; this suggested that as the internal tide radiates away from the generation site it transfers energy down to the diurnal mode via PSI. Carter & Gregg (2006) came to a similar conclusion based on another set of data from the Hawaii area. Through a bicoherence analysis they were able to demonstrate that the strength of the M2/2 diurnal mode was correlated to the strength of the M2 internal tide. Alford et al. (2007) estimated the reduction in energy flux as the internal tide crosses the critical latitude to be 10-20%, which is significantly smaller than the 62% in MW; they suggested that the presence of higher modes de-tuned the resonance and leads to a reduction in efficiency. Young et al. (2008) (henceforth YTB) constructed an analytical model based on a multiple time scale approach, where the time scale over which evolution of the PSI instability takes place was slower than the inertial time scale f−1. With this assumption the Boussinesq equations were simplified into a PDE that describes the evolution of the PSI wave. For simplicity the latitude was restricted to the critical latitude 28.9◦, so that the daughter waves are near-inertial. The stability of a plane wave internal tide in a linearly stratified ocean was considered, and using realistic parameters the growth rate of the unstable near-inertial mode agreed with MW in order of magnitude. YTB also considered the stability problem in a more realistic density stratification; using Gill’s model for the ocean stratification, the most unstable normal modes for the near-inertial field were found to be concentrated just under the mixed layer. 10 1.3 Goals The goal of this thesis is to expand on the model of PSI derived by YTB; in particular we will allow the Coriolis force to vary with latitude, and assess its effect on PSI. The rational behind this modification is that over the time-scale for PSI (10 days), there is significant propagation of the near-inertial waves, and thus the constant Coriolis parameter assumption may not be justified. In addition, we would like to investigate the effect of dissipation to determine if it is crucial to the instability, in particular to the understanding to vertical scale selection that appeared in the numerical simulations. 11 Chapter 2 Mathematic Description of PSI A familiar example of resonance occurs when a simple harmonic oscillator is forced at its natural frequency. Here we will first consider the following equation, which describes an oscillator forced by a weak ( 1) periodic driver with frequency 2: y′′ + ω2y = y∗ e2it. (2.1) The presence of y∗ (the complex conjugate of y) in the forcing term indicates that the forcing is moderated by feedback from the oscillator. Since the natural frequency for the oscillation is ±ω, the homogenous solutions are eiωt and e−iωt. These two solutions will result in forcing terms that are proportional to ei(2−ω)t and ei(2+ω)t respectively on the RHS of (2.1). As long as ω 6= 1, the forcing will be at a frequency that differs from the natural frequency, and hence no resonance will occur. On the other hand, when ω = 1 both the homogenous solution and the corresponding forcing will contain eit, and in this case resonance will occur. One can analyse (2.1) through a multiple time scale approach. As the coupling is weak due to the small parameter , we would expect the amplitude of y to grow over a slow time scale T = t. If we seek an ansatz to (2.1) of the form y(t, T ) = A(T )ei(t−pi/4), 12 we will have: 2iA′(T )ei(t−pi/4) +O(2) = A(T )ei(t+pi/4), (2.2) Equating the O() terms leads to a differential equation for A(T ): dA dT = 1 2 A (2.3) The above equations can be considered as a toy model for PSI, as the natural frequency of the oscillator is half the frequency of the driver, and therefore the resonant condition ω1 +ω2 = ω3 (where ω1 = ω2 = ω/2) is satisfied. In addition the model predicts that the amplitude of y(t) will grow as et, indicating that the driver is transferring energy into the oscillator, and the amplitude grows exponentially over the slow timescale. With this example in mind, we will now derive the model equations for PSI. The PSI model we will be using here is identical to the one derived by Young et al. (2008) except for two modifications. The first is the inclusion of the latitudinal dependence of the Coriolis force via a linearization of the Coriolis parameter f = f0 + βy. In addition, we will include dissipation into our model, which was not considered by Young et al. As these two modifications can be easily included in the rigorous asymptotic derivation of the governing equations outlined in Young et al, we will not repeat the derivation here; instead we will take a more heuristic approach by assuming a priori that the PSI modes will have frequency M2/2, and derive the governing equations. 2.1 Derivation of the PSI model The starting point is the hydrostatic Boussinesq equations linearized around the background flow fields (U, V,W,B, P ), which denote the velocities, buoyancy, and 13 pressure respectively. Since we are interested in the PSI in a mode-1 internal tide with frequency 2f0 = M2, the background fields will be proportional to e−2if0t and have amplitudes depending only on the spatial coordinates (x, y, z), which respec- tively denote the coordinates in the east-west, north-south, and vertical directions. Notice that our model is not a complete triad resonance, since the amplitudes of the background fields are time independent, and do not diminish as energy is being extracted from the internal tide. It can be viewed as a linearized model for the triad resonance, where the amplitudes of the daughter waves are small, and therefore the feedback on the internal tide (also referred to as the ‘pump wave’) can be neglected. With (u, v, w, b, p) denoting the perturbations associated with the NIO, we have: ut + Uux + V uy +Wuz + uUx + vUy + wUz − fv + px + µ(−∇2)nu = 0, (2.4) vt + Uvx + V vy +Wvz + uVx + vVy + wVz + fu+ py + µ(−∇2)nv = 0, (2.5) −b+ pz = 0, (2.6) ux + vy + wz = 0, (2.7) bt + Ubx + V by +Wbz + uBx + vBy + wBz + wN20 = 0, (2.8) where ∇2 = ∂2x + ∂2y + ∂2z . The last terms in (2.4) and (2.5) describe the dissipative process, and the exponent n = 0, 1, 2 corresponds to the linear drag, eddy diffusion, and hyperdissipation model respectively. In the ocean, the presence of stratification implies that the turbulence is not necessarily isotropic, and horizontal mixing occurs much more readily than mixing vertically due to stable vertical stratification. For example, the measured eddy horizontal diffusivity in a dye-release experiment in open ocean ranged from 10−2 to 103m2s−1 depending on the scale (Ledwell et al., 1998), while the vertical eddy diffusivity is of the order of 10−5m2s−1 (Ledwell et al., 1993, 1998); therefore in the eddy diffusion model, the dissipation term is split into 14 its vertical and horizontal components: µ(−∇2)nu = µH(−∇2H)nu+ µV (−∂zz)nu , (2.9) where ∇2H = ∂2x+∂2y , and µH and µV are the horizontal and vertical eddy diffusivity respectively. We will ignore the diffusion of buoyancy b. We can further simplify the equations using the assumptions stated in YTB. Since the background flow has a much larger vertical scale, any z-derivatives of the background variables can be ignored in (2.4) and (2.5). YTB assumed the buoyancy perturbations to be small compared to N2, and therefore the coupling terms (prod- uct of pumpwave and near-inertial variables) in(2.8) are small compared to wN2, and it will be the time derivative bt providing the balance in (2.8). Introducing the complex variables U = u+ iv, ξ = x+ iy, (2.10) the horizontal derivatives can then be written as: ∂ξ = 1 2 (∂x − i∂y), ∂ξ∗ = 12(∂x + i∂y). (2.11) We can then simplify (2.4) to (2.8) to: Ut + UUx + V Uy +WUz + (Ux + iVx)u+ (Uy + iVy)v + ifU + 2pξ∗ + µH(−∇2H)nU + µV (−∂2zz)nU = 0, (2.12) −b+ pz = 0, (2.13) ux + vy + wz = 0, (2.14) bt + wN20 = 0. (2.15) Since PSI primarily transfers energy to waves with half the frequency, we seek os- cillatory solutions with frequency f0; more specifically, we will assume that the 15 solutions are of the form: U = LA(x, y, z, t)e−if0t, (2.16) where L = ∂zf20N −2∂z. The other near-inertial fields can then be obtained from (2.13) to (2.15): w = −f20N−2Aξze−if0t + c.c., (2.17) b = if0Aξze−if0t + c.c., (2.18) p = if0Aξe−if0t + c.c. (2.19) It is worth noting that while obtaining b though (2.15), we have made use of the assumption that the near-inertial period f−10 is much larger than the time scale on which b evolves, and the variation of the amplitude can be ignore. b can then be obtained by integrating wN2 with respect to time. We will also linearize the Coriolis parameter f around the critical latitude: f = f0 + βy. (2.20) Once the pressure p is found, it can be substituted, together with the expres- sion of U into (2.12). The first term in (2.12), Ut, becomes: Ut = −if0LAe−if0t + LAte−if0t. (2.21) The main contributors to the slower tendency LAt arise from the terms that have the same phase (otherwise there are rapid cancellations over the diurnal cycle) i.e. the resonant terms. The others, referred to as non-resonant terms (NRT) can then be discarded. Among the coupling terms, only (Ux + iVx)u + (Uy + iVy)v are resonant: (Ux + iVx)u+ (Uy + iVy)v = 12U(Ux + Vy + iVx − iUy) + 12U∗(Ux − Vy + iVx + iUy) = (U + iV )ξ∗LA∗e−if0t +NRT, (2.22) 16 We will denote (U + iV )ξ∗ by 12Υ. Ignoring the NRT, the amplitude equation (2.12) then becomes: LAt + iβyLA+ 12 if0∇2HA+ µH(−∇2H)nLA+ µV (−∂zz)nLA+ 12ΥLA∗ = 0. (2.23) 2.1.1 The pumpwave train As in Young et al, we will consider a pumpwave that is periodic in the horizontal directions with horizontal wave numbers (k, l), and frequency ω = 2f0 + 2σ, where 2σ 2f0 is defined as the de-tuning frequency. Assuming that the pressure field associated with the pumpwave is given by: P = a cos(kx+ ly − ωt)p1(z) (2.24) where p1(z) is the eigenfunction of the first vertical mode. We can then retrieve the corresponding horizontal velocities through the linear Boussinesq equations. Defin- ing φ = kx+ ly − ωt, we have: U = [ akω ω2 − f20 cosφ− alf0 ω2 − f20 sinφ ] p1(z) , (2.25) V = [ alω ω2 − f20 cosφ+ akf0 ω2 − f20 sinφ ] p1(z), (2.26) Plugging U and V into the definition of Υ, and retaining only the parts that results in resonance: Υ = 2(U + iV )ξ∗ = ia(k + il)2 2(ω − f0) e i(kx+ly−2σt)p1(z). (2.27) We can ignore the argument of (k + il)2 since it just represents a constant phase shift of a. Defining υ = a(k2 + l2) 2(ω − f0) ≈ a(k2 + l2) 2f0 (2.28) 17 as the amplitude of the pumpwave, we have: Υ = iυei(kx+ly−ωt)p1(z) . (2.29) In the numerical model of MacKinnon & Winters (2005), the ocean was assumed to be linearly stratified and have a constant buoyancy frequency: N(z) = N0. In this case the first vertical mode will be: p1(z) = cos ( pi H z ) , (2.30) where H is the depth of the ocean. The dispersion relationship for internal gravity waves then results in a constraint on the horizontal wavenumbers k and l: k2 + l2 = ω2 − f20 N20 − ω2 pi2 H2 ≈ 3pi 2f20 N20H 2 , (2.31) where the buoyancy frequency N0 is assumed to be much larger than f0. 2.1.2 The simplified amplitude equation We now look for solutions of the form A ≈ cos(mz) e−iσt [ eik1xA1(y, t) + eik2xA∗2(y, t) ] , (2.32) where m is the local vertical wave number, and the horizontal wave numbers satisfy the resonant condition: k1 + k2 = k . (2.33) If we assume that the vertical scale of the instabilities, m−1, is small compared to scale over which the buoyancy frequency and the vertical structure of the pump wave vary, N and p1(z) will be approximately constant. Taking N(z) to be some charac- teristic frequency N0 and p = 1, we can introduce the scaled vertical wavenumber m̃ ≡ f0m N0 , (2.34) 18 and the differential operator L becomes: LA = −m̃2A , (2.35) Substituting (2.29) and (2.32) into the amplitude equation (2.23) and separating the two harmonic components, we obtain a pair of coupled partial differential equations for the amplitudes Aj(y, t). With the horizontal length scale (λ) and time scale (τ) defined as λ = ( f0 2m̃2β )1/3 , τ = 1 λβ , (2.36) the non-dimensionalized equations are: iA1t̃ = −A1ỹỹ + ( ỹ − b̃− ã− σ̃ ) A1 + υ̃eil̃ỹA∗2 − iµ̃H [ κ1A1 + (−1)nA1ỹ2n ] , (2.37) −iA2t̃ = −A2ỹỹ + ( ỹ − b̃+ ã− σ̃ ) A2 + υ̃e−il̃ỹA∗1 + iµ̃H [ κ2A2 + (−1)nA2ỹ2n ] . (2.38) In the above equations, the t̃ = λβt and ỹ = y/λ are the non-dimensionalize in- dependent variables, and k̃j = λkj and l̃ = λl are the non-dimensionalized wave numbers. The non-dimensionalized pump wave amplitude υ̃, de-tuning frequency σ̃ and dissipation parameter µ̃ are given by: υ̃ = υ 2βλ , σ̃ = σ βλ , µ̃ = µH λ2n+1β . (2.39) The parameters ã and b̃ are defined as: ã = −λ 2(k21 − k22) 2 , b̃ = −λ 2(k21 + k 2 2) 2 . (2.40) The effective wave number is defined by: κi = k̃2ni ( 1 + µVm 2n µHk2ni ) (2.41) 19 (2.37) and (2.38) can be further simplified by removing some parameters, as their effects on the system can be understood easily. For example b̃+ σ̃ can be removed via a shift in spatial coordinates ỹ − b̃ − σ̃ → ỹ together with a suitable constant phase shift in A1,2, and therefore the effect of these two parameters is limit to a shift of the origin. If we assume that k1 = k2, the dissipative term −iµ̃Uκi can also be removed by a change of variable: Ai → e−µ̃HκitAi ; this suggests that the effect of dissipation in the zonal and vertical direction is to suppress the amplitude of the PSI waves without altering the meridional structure, and hence we can ignore it. Finally we can also remove the parameter ã via (A1, A2) → ( eiãtA1, e−iãtA2 ) , and therefore the effect of ã is limited to a frequency shift. With the above simplifications, as well as dropping the tildes, we have: iA1t = −A1yy + yA1 + υeilyA∗2 − iµH(−1)n∂2ny A1, (2.42) −iA2t = −A2yy + yA2 + υe−ilyA1 + iµH(−1)n∂2ny A2. (2.43) (2.42) and (2.43) will serve as the basis for the rest of this study. 2.2 Parameter estimation We will close this chapter by providing a summary of the oceanographically relevant values for the physical parameters used in our scaling. In table 2.1, the pump wavenumber l and amplitude υ are similar to the ones used in YTB. The buoyancy frequency of 0.0087s−1 is a value that is typically found in regions just underneath the mixed layer. These values, together with a vertical wave number of 100m, give a time scale of 30.5 days and a length scale of 19km; this length scale corresponds to a change of 0.171◦ in latitude. 20 M2 tidal frequency 1.405× 10−4s−1 Coriolis frequency, f0 = 12M2 7.02× 10−5s−1 Planetary vorticity gradient, β 2.003× 10−11 m−1 s−1 Buoyancy frequency, N0 2pi/720s = 0.0087s−1 Vertical wavenumber, m 2pi/100m Pump wavenumber, l 2pi/125km Pump amplitude υ 1/5days Length scale, λ = (f0N20 /2f 2 0m 2β)1/3 19km Time scale, 1/βλ 30.5 days Eddy viscosity scale, λ3β 1.36× 102m2s−1 Non-dimensional pump amplitude, υ̂ = υ/(2βλ) 3.04 nondimesnional pump wavenumber, l̂ = λl 0.946 Table 2.1: Parameter values 21 Chapter 3 Dissipative Problem In this chapter we will explore the solutions to the PDEs (2.42) and (2.43) derived in the previous chapter. We will mostly focus on the regular eddy dissipation model by taking n = 1, which results in iA1t = −(1− iµ)A1yy + yA1 + υeilyA2, (3.1) −iA2t = −(1 + iµ)A2yy + yA2 + υe−ilyA1. (3.2) A set of energy equations can be derived for (3.1) and (3.2). The total energy density E can be defined as: E ≡ 12(u2 + v2) = 12 |LA|2 . (3.3) E is proportional to |A1|2+|A2|2, and thus En = |An|2 for n = 1, 2 can be interpreted as energy density for our equations. From (3.1) and (3.2), we can derive the following energy equations: E1t + J1y = −D1 + S, (3.4) E2t + J2y = −D2 + S, . (3.5) 22 with the meridional energy flux defined as Jn(y, t) ≡ i(−1)n ( A∗nAn y −AnA∗n y )− µH (An yA∗n +A∗n yAn) . (3.6) The energy loss due to dissipation is given by Dn(y, t) ≡ 2µH |A2n y| , (3.7) and finally the energy source due to PSI is: S = −iυ ( eilyA2A ∗ 1 − e−ilyA∗2A1 ) . (3.8) 3.1 An initial value problem The controlling parameters in (3.1) and (3.2) are the dissipation (µH), the amplitude of the pumpwave/forcing (υ), and the meridional wave number of the pumpwave (l). In the absence of the pumpwave and effect of dissipation, the equations reduce to the Schrödinger equation with a linear potential, which describes a particle falling under the influence of gravity. This is can be explained via a slowly varying wave approximation. Assuming that µH and υ are small (and thus can be ignored), the A1 equation takes the form: iA1t = −A1yy + yA1 (3.9) For the slowly varying wave approximation, we will assume A1 = eiφ, where φ = φ(y, t) is the phase. Defining the frequency ω̂ = −φt and wave number k̂ = φy, the dispersion relation for (3.9) is ω̂ = k̂2 + y. (3.10) Differentiating (3.10) with respect to y, and noting that ω̂y = −φty = −k̂t, we now have a PDE for k̂: k̂t + 2k̂k̂y = −1, (3.11) 23 which can be solved using the method of characteristics. The solution is given by k̂ = −t + k̂0 on the characteristic curves y − y0 = −t2 + 2k̂0t, where k̂0 and y0 are respectively the initial wave number and position of the wave packet. The form of the characteristic curves suggest that the wave packet will indeed follow a ballistic trajectory, with the wave-packet having velocity dy/dt = 2k = 2(t− k̂0); this suggest that for t > k̂0 the wave-packet will have negative velocity, and will be moving south. Eliminating t from the solution results in: k̂2 = k̂20 + (y0 − y), (3.12) which indicates that the wave number increases as the wave-packet travels south. With the above observations, we can now formulate an initial value problem for PSI. Consider a perturbation in the near-inertial wave field in the form of a Gaussian wave packet with wave-number k0 > 0 and : A1(y, 0) = e−(y−y0) 2/σ cos(k0y) = 1 2 ( e(y−y0) 2/σeik0y + e(y−y0) 2sigmae−ik0y ) (3.13) According to the analysis above, we expect the perturbation to split into two wave packets. The component with e−ik0y travels south, while the component with eik0y will first travel north up to y = k20 + y0 before turning south. In the presence of the pump-wave (i.e. υ > 0), the near-inertial wave-packets will grow due to PSI, and we would like to understand how the growth depends on the parameters. 3.1.1 Numerical scheme As with other wave-type problems on an infinite domain, We will demand that the amplitude of the near-intertial waves decay to 0 as y → ±∞; however there are some practical issues with truncating the problem to a finite domain [yL, yR] for the numerical calculations as the boundary conditions (BC) at infinity are now applied 24 on a finite interval. For the right boundary yR, the wave packet will be reflected before it reaches yR provided that the initial kinetic energy is not too high, and thus we can simply set A1(yR, t) = A2(yR, t) = 0. The situation at the left boundary yL is complicated by the presence of −Aiyy + yAi in (3.1) and (3.2); the slowly varying wave theory predicts that the wave number increases as √|y| as y → −∞ (see (3.12)). This presents us with a dilemma: on one hand the domain needs to be large in order for the wave amplitude to decay sufficiently to justify imposing a zero Dirichlet condition, but as the domain is extended by moving yR towards −∞ the grid needs to be refined to resolve the increasingly rapid oscillations. One way to resolve this issue is to impose an absorbing boundary condition (ABC) in conjunction with truncation of the domain. There are two main types of ABCs: PDE-based and material based (for a brief review, see Zheng (2007)). PDE- based ABC relies on the ability to factor the differential operator, which in turn can be used to derive a radiative boundary condition that only allows out-going waves at the boundary. On the other hand, the material based ABC incorporates a sponge- like layer near the endpoint to absorbed outgoing waves. One class of material based ABC is the perfectly matched layer (PML), which was used by Berenger (1994) for the study of Maxwell’s equations in electro-magnetics. It was also used sucessfully in both linear and non-linear Schrödinger’s equation with variable potentials (Zheng, 2007). As it is not clear how our equations can be factored to derive a radiative BC on the boundary, we will adopt the PML technique as outlined in Zheng (2007) for this numerical study. The main idea behind the method is to extend the computa- tional domain from [yL, yR] to [yLR, yR] with yLR < yL. The PDE with the PML 25 is identical to the original PDE in the original domain of interest, while the PDE is modified in the sponge layer [yLR, yL] by introducing an additional dissipative term that strongly damps out-going waves. The solutions are expected to decay exponentially, and hence a zero Dirichlet BC can be imposed at y = yLR. The primary function of the sponge layer is to ensure that the initial per- turbation is not reflected back into the boundary, and thus allowing us to consider cases where the dissipation is too small to damp out the perturbations before it reaches the boundary. With the introduction of a PML, the second derivatives are modified via: A1yy → − 11− i ζ(y) d dy ( 1 1− i ζ(y)A1y ) (3.14) A2yy → − 11 + i ζ(y) d dy ( 1 1 + i ζ(y) A2y ) . (3.15) The absorption function ζ(y) is chosen to be: ζ(y) = ζ0(y − yl) 4 if yLR < y < yL 0 if yL ≤ y < yR (3.16) Equations (3.1) and (3.2) together with the initial condition A1(y, 0) = e−(y+25) 2/5 cos( √ 26y), A2(y, 0) = 0 (3.17) are solved using a Crank-Nicolson finite difference scheme, which has the advantage of being unconditionally stable (Garcia, 1999). We will now present a sample solution, where the computational grid is uniform both in space and time, with time-step ∆t = 0.00005 and 4000 spatial grid- points over [−100, 50]. A sponge layer is located at [−100, 85] with the ζ0 = 0.01. The amplitudes of the near-inertial wave field, |A1| and |A2|, are shown in Figure 3.1. We can see that the wave packet splits up into two components that travel in 26 opposite directions. The northward propagating packet follows a ballistic trajectory while spreading out. At the same time we can see that as the wave packet moves through the domain it excites a standing wave that grows in amplitude over time. The standing wave is concentrated near the origin (the critical latitude), and the amplitude rapidly drops off to the south. This is not surprising as the effect of dissipation will be more potent as the wave-number increases to the south. ! " #$%&'( ) ' !*+ !,+ !-+ !.+ !/+ !0+ !1+ !2+ !)+ + )+ + 2 0 . , )+ )2 )0 ). ), 2+ " #3%&'( 2 ' !*+ !,+ !-+ !.+ !/+ !0+ !1+ !2+ !)+ + )+ + 2 0 . , )+ )2 )0 ). ), 2+ Figure 3.1: The results from a numerical computation of the PSI model. |A1| and |A2| are plotted in panels (a) and (b) respectively. The parameters used are υ = 0.25, l = 0, and µH = 0.01, with the spatial domain being [−85, 50] and a sponge layer located at [−100,−85]. Only the data from [−90, 10] is shown here. 3.2 Normal modes The standing wave solutions can be best understood via an analysis of normal modes. The standing wave solution in the IVP decays exponentially to the south of the origin, and therefore we do not need to include a sponge layer in our problem. The domain can be truncated at a point sufficiently far south, and a zero Dirichlet boundary condition can be imposed. Assuming An = Aneγt, where γ = γR + iγI , 27 (3.1) and (3.2) become: (1− iµ)A1yy = (y − iγ)A1 + υeilyA2, (3.18) (1 + iµ)A2yy = (y + iγ)A2 + υe−ilyA1. (3.19) This system is now an eigenvalue problem, and the boundary conditions are the same as the PDE: |An| → 0 as y → ±∞. 3.2.1 Numerical solution To investigate the effect of the three parameters υ, l, and µ on the eigenvalues, we will compute the eigenvalues numerically. The second derivatives in (3.18) and (3.19) are discretized using a centered final difference scheme, and the eigenvalue problem becomes a matrix eigenvalue problem; the matrix is large but sparse, which can be handled via the eigs function in MATLAB. eigs is an iterative numerical scheme that estimates the eigenvalues for matrices, and it is useful when only part of the spectra is of interest. We will compute the spectra over [yL, 50] for some yL, and impose zero Dirichlet boundary condition at both end-points. We will then take the limit as yL → −∞; if the eigenvalues converge to a limit, the limit can be interpreted as the eigenvalues of the problem over the infinite interval. We found that the eigenvalues and eigenfunctions only differ slightly with the zero Dirichlet BC is replaced by a sponge layer at the southern boundary, and therefore we will ignore the sponge layer in this section. Zonally propagating pump-wave: l = 0 Figure 3.2 shows the result for a series of calculations with υ = 1, and l = 0. The eigenvalue problem is solved over a domain of [−100, 50], and the real part of the 28 normalized eigenfunctions <(A1) over [−60, 10] are plotted on the left column of figure 3.2. µ = 1, 10−1, 10−2, 10−4, and 0 were used in the calculations; the plots are arranged in decreasing order of µ from top to bottom. One can see that in all cases, the eigenfuction is oscillatory to the south (when y < 0), while to the north of the turning point y = 0, the eigenfunction quickly decays to 0, indicating that the near-inertial waves are evanescent as expected. On the bottom row, the numerically computed eigenfunction (blue curve) is almost indistinguishable from the Airy function (red curve) in the absence of dissipation; in §3.2.2 we will show that the Airy function is the exact solution to the inviscid problem. When dissipation is increased to 10−4, the eigenfunction still resembles the airy function, but the amplitude to the south is significantly reduced. As dissipation is increased further, the wave becomes more and more localized to the turning latitude. The spectra for each case is plotted on the right column. In the inviscid limit µ = 0 (bottom panel) the eigenvalues are either purely real (growing or decaying modes), or purely imaginary (neutral modes). The fastest growing mode has a growth rate close to 1. As the dissipation increases, to µ = 10−2, the majority of the modes become decaying, with a handful of modes that remain unstable. When the dissipation is further increased to 1, all but one of the modes become stable, and the growth rate of the unstable mode is now significantly lower then the inviscid case. The eigenvalues becomes insensitive to yL for yL < −100, and therefore the computed spectra should be a good approximation to the spectra on the infinite interval. 29 −60 −50 −40 −30 −20 −10 0 10 −1 0 1 y −1 0 1−1 0 1 Re A1 −1 0 1−1 0 1 −5 0 5 −1 0 1 !I −1 0 1−1 0 1 !r −1 0 1−1 0 1 µ=0 µ=10−4 µ=10−2 µ=10−1 µ=1 Figure 3.2: Effect of dissipation on the spectrum: l = 0: The dissipation parameter used (from top to bottom) are µ = 1 , 10−1 , 10−2 , 10−4 , and 0. On the left column the real part of the eigenfuctions associated with the fastest growing modes are plotted in blue, while in the last two panel the airy function Ai(y) is also plotted in red. On the right column the spectra are plotted, with the abiscissa and ordinate being respectively the imaginary part and the real part (the growth rate) of the eigenvalues. υ = 0 for all calculations. More general case: l 6= 0 We will now turn to the case where l is no longer zero. In figure 3.3, the re- sults from a set of computations with υ = 1 and l = 0.1 is shown, with µ = 1, 10−1, 10−2, 10−3, 10−4, and 0. The computations are repeated over the domain [yL, 20], where yL was varied to determine the effect of truncation. The results with yL = −100 are plotted in blue, while the results with yL = −150 are plotted in red. Notice that only parts of the eigenfunctions over [−80, 10] are shown, and the scale for the frequency of the spectra plot is different from figure3.2. For large values of dissipation (top two rows in figure 3.3), the results from both values of yL agree, showing that neither the eigenfunction nor the physically 30 −80 −60 −40 −20 0 −1 0 1 y −1 0 1 −1 0 1 −1 0 1 Re A1 −1 0 1 −1 0 1 −2 −1 0 1 2 −1 0 1 !I −1 0 1 −1 0 1 −1 0 1 !r −1 0 1 −1 0 1 µ=10−4 µ=0 µ=10−3 µ=10−2 µ=10−1 µ=1 Figure 3.3: Effect of dissipation on the spectrum: l = 0.1: The dissipation parameter used (from top to bottom) are µ = 1, 10−1, 10−2, 10−3, 10−4, and 0. On the left column the real part of the eigenfuctions associated with the fastest growing modes are plotted; in blue are the results from using a domain of [−100, 20], while the curves in red are based on a domain of [−150, 20]. The corresponding spectra are respectively plotted with blue dots and red circes. The abiscissa and ordinate are respectively the imaginary part and the real part (the growth rate) of the eigenvalues. υ = 1 for all calculations 31 interesting portion of the spectra is affected by the change in yL. As in the l = 0 case, the eigenfunction decays quickly to zero due to dissipation, and there are a small number of discrete unstable modes. When dissipation is decreased to 10−2, the two computations agree for the fastest growing mode, but the spectra plot reveals that the domain length starts to affect the computation. The computation with yL = −100 suggests the presence of unstable modes that have non-zero frequency (i.e. γI 6= 0), while the computation with yL = −150 suggests that all unstable modes have zero frequency. As the dissipation is decreased further (bottom three rows), the picture is very different from the l = 0 case. The spectra plots show that most unstable modes now have non-zero frequency, and the frequency increases as the yL →∞. The plots of the eigenfunction reveal that the waves are no longer localized to the turning latitude; instead, the bulk of the wave is located at the center of the domain, and moves as the boundary moves southwards. We will refer these modes as ‘boundary modes’. To get a better understanding of the transition between the boundary and localized modes as µ changes, we performed a series of computations where yL and µ were varied systematically: the boundary was moved from yL = −80 to yL = −200, and for each series µ was varied from 10−4 to 1, with υ = 1 and l = 0.1. The results are presented in figure 3.4, where the frequency (γI) and growth rate (γR) of the most unstable mode are plotted as a function of dissipation. When the dissipation is small (µ ≈ 10−3), the growth rate (bottom panel) depends on yL, but the sensitivity is reduced as dissipation vanishes; the graph suggests that as we approach the inviscid limit, the growth rate will converge to a value that is independent of yL; however, while the frequency (top panel) is insensitive to dissipation, the frequency remains 32 10−4 10−3 10−2 10−1 100 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 µ !R Predicted yL=−80 yL=−100 yL=−120 yL=−140 yL=−160 yL=−180 yL=−200 10−4 10−3 10−2 10−1 100 0 0.2 0.4 0.6 0.8 1 !I −200 −180 −160 −140 −120 −100 −80 0.5 1 1.5 yL !I,R !I !R 0.1 yL 0.5 Figure 3.4: Effect of computational domain on the eigenvalues. The frequency (γI , top panel) and growth rate (γR, bottom panel) of the most unstable modes are plotted as a function of µ on a semi-log scale. The computation domain [yL, 20] ranges from yL = −80 to yL = −200, and the results from each series are marked by dots. The predicted growth rate from (3.74) is also plotted in the bottom panel as the green curve. Inset: γR and γI as a function of yl in the inviscid limit. υ = 1 and l = 0.1 for all calculations. 33 dependent on yL. The effect of yL on the eigenvalues for the inviscid problem (µ = 0) is shown in the inset of figure 3.4. The growth rate depends very weakly on yL, but the frequency increases approximately as l √ yL (plotted as the green dashed line). The eigenfunctions for the unstable modes in the inviscid problem resemble the boundary modes plotted in the bottom panel of figure 3.3. A question regarding these boundary modes is whether they are artifacts re- sulting from grid instabilities, but this can be discounted since the rapid oscillations in the eigenfunctions are well resolved, and a change in grid size does not affect the result. It is also not due to the zero Dirichlet boundary condition imposed at the southern boundary, as the eigenvalues and eigenfunction of the boundary modes are only affected minimally if the zero Dirichlet BC is replaced by a sponge layer or a radiative (Sommerfeld) type boundary condition that eliminates incoming waves at the boundary. Since the frequency of the boundary mode does not appear to converge, and the location of peak of the eigenfunctions changes as the yL is moved to the south, the conclusion is that these modes are un-physical. We will now return to figure 3.4. From the inviscid limit, the growth rate of these boundary modes rapidly decreases as dissipation is increased, with the mode from yL = −200 being most susceptible to dissipation. There is a sudden break in the graph when µ inceases to 10−2.5, and the numerically computed growth rates now converge to a single curve, and no longer depends on yL. The frequency also undergoes a similar transition: as the frequency increases pass the transition near µ ≈ 10−2.5, all the curves collapse to abruptly to 0. This marks the transition from the boundary modes to modes where both the eigenvalues and eigenfunctions are independent of the location of the boundary. These new modes are physically mean- ingful since the choice of the locaion of the boundary does not affect the structure 34 or growth rate of the mode. The eigenfunctions are similar to the localized modes in the top three panels of figure 3.3. A curious feature is that the effect of dissipation on these localized modes is not monotonic: while an increase in dissipation for larger values of µ stabilizes the system, the growth rate in fact increases with µ when dissipation is small (between 10−2.5 to 10−1.7); this result rather surprising, as one would normally expect dis- sipation to have a stabilizing effect. We will explore this behaviour via the energy equations in §3.2.3. We will now present a couple of plots in figure 3.5, which summarize the effect of the various parameters on the growth rate for υ = 4 (top panel) and υ = 6 (bottom panel). The two plots are qualitatively similar. For l = 0 (top series), γR ≈ υ when the dissipation is small, and decreases monotonically as µ increases. For l 6= 0, the non-monotonic effect of dissipation on growth rate observed earlier holds for a wide range of l. An increase in l reduces the growth rate, and for sufficiently large value of l, the normal modes become stable for small µ. For example, when l = 2.5 and υ = 4 (bottom series in the top panel of figure 3.5), the growth rate is negative for µ < 10−0.5. 35 10−2 10−1 100 101 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 µ γ R (a) υ = 4 10−2 10−1 100 101 −1 0 1 2 3 4 5 6 µ γ R A (b) υ = 6 Figure 3.5: Effect of the parameters on the growth rate. The growth rates are plotted as a function of µ for different values of pump wavenumber l (dots); from top to bottom, l = 0, 0.5, 1,1.5, 2, and 2.5. The dashed lines are the growth waves predicted by the short wave theory in §3.2.2(c.f. (3.74)). For all calculations, the eigenvalue problem was solve over the domain [−100, 20]. 36 3.2.2 Analytical results We will now develop some analytical results which provides support to our numeri- cally computed growth rates. The most important result is an analytical prediction for limit where 0 < l 1 and µ 1, as it confirms our numerical results showing that µ→ 0 stabilizes the system for a non-zero l. The problem with µ = l = 0 can be solvable exactly, which is then used as the basis of a perturbative solution with l and µ as a small parameter. A zonally propagating pump-wave in the inviscid limit We will first consider the inviscid limit (µ = 0), as well as a pump-wave that prop- agates zonally (l = 0), as in this regime the eigenvalue problem can be diagonalized and solved exactly:A1(y) A2(y) = υ −iγ −√υ2 − γ2 iγ + √ υ2 − γ2 υ c1Ai(y −√υ2 − γ2) c2Ai(y + √ υ2 − γ2) (3.20) The fundamental solutions are Ai(y± √ υ2 − γ2), where Ai(y) is the Airy function of the first kind. The argument of Ai(y) must be real for Ai(y) to remain bounded as y → ±∞, which implies that γ must be real and υ2 − γ2R ≥ 0. The latter condition implies that there is a continuous spectrum of normal modes, with growth rates satisfying the inequality: −υ ≤ γR ≤ υ . (3.21) The maximum growth rate is γRmax = υ, which agrees with the numerically ob- tained growth rates (c.f. figure 3.5 and the bottom panel of figure 3.3). The eigenfunctions associated with the fastest growing mode are A1(y) = Ai(y) and A2(y) = iAi(y). Using (3.6), It can be shown that the meridional energy flux is 37 identically 0. The interpretation is that the fastest growing near-inertial oscillations excited by a zonally propagating pump-wave will only transport energy zonally. Young et al. (2008) found the growth rate for a zonally propagating pump- wave at the critical latitude 28.9◦ to be (with the variables non-dimensionalized as the present study): γRmax = √ υ2 − (b+ σ)2 , (3.22) Curiously the factor b+σ in (3.22) is identical to the spatial shift we initially ignored in the IVP problem, yet it plays no role in the growth rate for our analysis. It is perhaps instructive to re-concile the two results. Discarding the β-effect term yAn, and re-inserting b+ σ into our equations, the eigenvalue problem becomes: A1yy = (−(b+ σ)− iγ)A1 + υ A2 , (3.23) A2yy = (−(b+ σ) + iγ)A2 + υ A1 , (3.24) Since y does not explicitly appear in the equations, we can seek solutions that are periodic in y: A1 = eil1y and A2 = e−il2y. To satisfy the resonant condition, the wave numbers must satisfy l1 + l2 = 0, and the eigenvalues are: iγ = i √ υ2 − [l21 − (b+ σ)]2 . (3.25) YTB considered the reduced 2-dimensional model in the absence of the β−effect, and seeked solutions that are uniform in y; this is equivalent to setting l1 in (3.25) to be 0 and we will recover their result (3.22). On the other hand, there is no reason to restrict ourselves to meridionally uniform solutions, and thus (3.25) is a more general result for the model of YTB. This equation suggests that maximum growth rate γmax = υ can be achieved when l1 = √ b+ σ, and therefore there is scale selection in the y-direction at the critical latitude. 38 This scale selection behaviour is consistent with our model. If we retain b+ σ in our governing equations, the corresponding eigenfuction is Ai (y − (b+ σ)) which has a turning point at y = b + σ; to the south of the turning point the Airy function is periodic with local wave-number √|y − (b+ σ)|, while to the north the wave becomes evanescent. Now if (b + σ) 1, then the critical latitude y = 0 lies to the south of the turning point; therefore the waves at the critical latitude will be periodic with wave-number √ (b+ σ), which is consistent with the analysis in the preceding paragraph. Short wave theory We can use the analytical solution to the µ = l = 0 problem in the previous section to make progress analytically for cases where the the parameters υ and l are non-zero but small. It is convenient to consider a change in variables: f(y) = A1e−ily/2 − iA2eily/2, ig(y) = A1e−ily/2 + iA2eily/2 (3.26) Our equations then become: f ′′ − lg′ − l 2 4 f + µg′′ + µlf ′ − µl 2 4 g − yf − (υ + γ)g = 0 , (3.27) g′′ + lf ′ − l 2 4 g − µf ′′ + µlg′ + µl 2 4 f − yg − (υ − γ)f = 0 . (3.28) In the previous section, the numerical results suggest that the unstable modes have frequency that vanishes (γI = 0), and therefore we will assume that γ, f and g in (3.27) and (3.28) are purely real. In addition we can easily verify that f(y) = Ai(y), g(y) = 0, and λ = υ are solutions to the above system. With a small parameter defined as = √ µ, as well as the assumption that l scales as 3/2 (i.e. l = l3/23/2), 39 (3.27) and (3.28) becomes f ′′ − l3/23/2g′ + 2g′′ − yf − (υ + γ)g = O(5/2) , (3.29) g′′ + l3/23/2f ′ − 2f ′′ − yg − (υ − γ)f = O(5/2) . (3.30) Inner solution We will seek an expansion of the form: f = Ai(y) + 1/2f1/2(y) + · · · , (3.31) g = 1/2g1/2(y) + g1(y) · · · , (3.32) γ = υ + γ1 + · · · . (3.33) Plugging (3.31) to (3.33) into (3.29) and (3.30) and gathering the O(1/2) terms, we have: f ′′1/2 − yf1/2 = 2υg1/2 , (3.34) g′′1/2 − yg1/2 = 0 . (3.35) (3.35) is satisfied by g1/2 = rAi(y), where r is a constant that is to be determined later. Since (Ai′)′′ = (Ai′′)′ = (yAi)′, a particular solution to (3.34) is f1/2 = 2υrAi ′(y) . (3.36) Moving on to O() terms in (3.30), we have: g′′1 − yg1 = γ1Ai, (3.37) which has the particular solution: g1 = −γ1Ai′(y) , (3.38) 40 To summarize, the inner solution is: f = Ai(y) + 1/22υrAi′(y) + · · · , (3.39) g = 1/2rAi(y)− γ1Ai′(y) + · · · , (3.40) The Outer Region As y → −∞, standard asymptotic formulae for the Airy function give (Abramowitz & Stegun, 1964): Ai(y) ∼ 1√ pi 4 √−y sin ( φ+ pi 4 ) and (3.41) Ai′(y) ∼ 4 √−y√ pi cos ( φ+ pi 4 ) , (3.42) where φ(y) = 23 |y|3/2 is the phase. Since the solution becomes increasingly oscil- latory, we anticipate that some of the derivatives we neglected in the inner region may become important. In addition, (3.42) suggests that the O() term in g(y) grows as y → −∞, which leads to a break in the order. Both of these observations point to the existence of an outer solution which we will need to match onto the inner solution (which we will relabel as fin and gin) given by (3.39) and (3.40). To proceed we will define Y = y as the outer variable; notice that in terms of the outer variable, f ∼ 1/4 [ 1 4 √−Y sin ( φ(Y ) 3/2 + pi 4 ) + 2υr 4 √−Y cos ( φ(Y ) 3/2 + pi 4 )] and (3.43) g ∼ 3/4 [ r 4 √−Y sin ( φ(Y ) 3/2 + pi 4 ) − γ1 4 √−Y cos ( φ(Y ) 3/2 + pi 4 )] , (3.44) 41 We have scaled the inner solutions to remove the constant factor of 1/ √ pi. Motivated by the form of the inner solution, we seek a WKBJ-type solution in the outer region: fout = 1/4 [ R1(Y ) sin ( φ(Y ) 3/2 ) +R2(Y ) cos ( φ(Y ) 3/2 )] and (3.45) gout = 3/4 [ S1(Y ) sin ( φ(Y ) 3/2 ) + S2(Y ) cos ( φ(Y ) 3/2 )] . (3.46) In terms of the outer variable, the scaled equations are: 2 d2fout dY 2 − l3/25/2 dgout dY + 4 d2gout dY 2 − Y f − (2υ + γ1)g = O(5/2) , (3.47) 2 d2gout dY 2 + l3/2 5/2dfout dY − 4d 2fout dY 2 − Y g + γ1f = O(5/2) , (3.48) We will now proceed with the asymptotic expansion by substituting (3.45) and (3.46) into (3.47) and (3.48). The leading order equation gives: dφ dY = ±√−Y (3.49) for both equations. We will take the negative root as its integral is consistent with the definition of the phase function φ(Y ) = 23 |Y |3/2. At the next order, the equations contain both sine and cosine terms, and their coefficients must vanish for the equations to hold. We can therefore extract four equations: 2 dR1 dY dφ dY +R1 d2φ dY 2 = 2υS2 (3.50) 2 dR2 dY dφ dY +R2 d2φ dY 2 = −2υS1 (3.51) 2 dS1 dY dφ dY + S1 d2φ dY 2 = −γ1R2 − l3/2 dφ dY R1 − dφ dY 2 R2 (3.52) 2 dS2 dY dφ dY + S2 d2φ dY 2 = +γ1R1 − l3/2 dφ dY R2 + dφ dY 2 R1 (3.53) These can be simplified by defining η = 4 √ 8υ √−Y , and R̂n, Ŝn = Rn√η, Sn√η. 42 With dφdY = − √−Y = −η/ 4√8υ, d2φ dY 2 = 4 √ 8υ/(2η), (3.50) to (3.53) simplifies to: R̂1η = 2υ 4 √ 8υ Ŝ2 (3.54) R̂2η = − 2υ4√8υ Ŝ1 (3.55) Ŝ1η = 1 4 √ 8υ [ −γ1R2 + η4√8υ l3/2R1 − η2√ 8υ R2 ] (3.56) Ŝ2η = 1 4 √ 8υ [ γ1R1 + η 4 √ 8υ l3/2R2 + η2√ 8υ R1 ] (3.57) Eliminating Ŝ1 and Ŝ2 from (3.54) to (3.57), and defining A = γ1 √ υ 2 and B = l3/2 2 4 √ υ 2 , we have a pair of coupled second order ODE: R̂1ηη − ( η2 4 +A ) R̂1 = BηR̂2 (3.58) R̂2ηη − ( η2 4 +A ) R̂2 = −BηR̂1 (3.59) It turns out that the pair can be uncoupled by defining F1 = R̂1 + iR̂2 and F2 = R̂1 − iR̂2, the differential equations for F1 and F2 are: F1ηη − ( 1 4 (η + 2iB)2 +B2 +A ) F1 = 0 , (3.60) F2ηη − ( 1 4 (η − 2iB)2 +B2 +A ) F2 = 0 , (3.61) both of which have parabolic cylinder functions as solutions. For the outer solutions fout(Y ) and gout(Y ) we demand that they decay to zero as Y → −∞, which is equivalent to |F1,2| → 0 as η →∞. The standard solutions to the parabolic cylinder equations are given by U(a, x) and V(a, x), and only U(a, x) decays as x → ∞ (Abramowitz & Stegun, 1964). Hence we have: F1(η) = c1U(B2 +A, η + 2iB) and F2(η) = c2U(B2 +A, η − 2iB) , (3.62) where c1 and c2 are constants that are determined by matching with the inner solution. 43 Matching To match up the inner and outer regions, the outer solution should take the form of (3.45) and (3.46) as Y → 0. To facilitate matching, we will expand (3.45) and (3.46) fin ∼ 1/4 [( 1 4 √−Y + 2υr 4 √−Y ) sin ( φ(Y ) 3/2 ) + ( 1 4 √−Y − 2υr 4 √−Y ) cos ( φ(Y ) 3/2 )] , (3.63) gin ∼ 3/4 [( r 4 √−Y + 2υ 4 √−Y ) sin ( φ(Y ) 3/2 ) + ( r 4 √−Y − 2υ 4 √−Y ) cos ( φ(Y ) 3/2 )] , (3.64) with the constant factors √ 2/2 being scaled out. Comparing the fin and fout from (3.45), we will need the coefficients of the sine and cosine terms to match as Y → 0. Remembering that R̂n = Rn√η and √−Y = η/ 4√8υ, it suggests that near Y = 0, the matching conditions are: R̂1(η) ∼ 4 √ 8υ + 2υrη and R̂2(η) ∼ 4 √ 8υ − 2υrη , =⇒ R̂1(0) = R̂2(0) = 4 √ 8υ and R̂′1(0) = −R̂′2(0) = 2υr . (3.65) Similarly, comparing gin and gout gives Ŝ1(η) ∼ r 4 √ 8υ − γ1η and Ŝ2(η) ∼ r 4 √ 8υ + γ1η . =⇒ Ŝ1(0) = Ŝ2(0) = r 4 √ 8υ and Ŝ′1(0) = −Ŝ′2(0) = −γ1 . (3.66) With the help of (3.54) to (3.57), we can see that the boundary conditions given by (3.66) are automatically satisfied as long as the boundary conditions in (3.65) are satisfied. The conditions R̂1(0) = R̂2(0) and R̂′1(0) = −R̂′2(0) from (3.65) in turn require that F1(0) = iF2(0) and F ′1(0) = −iF ′2(0) (3.67) 44 or equivalently, U(B2 +A, 2iB) −iU(B2 +A,−2iB) iU ′(B2 +A, 2iB) −U ′(B2 +A,−2iB) c1 c2 = 0 0 , (3.68) For the system (3.68) to have a non-trivial solution, the determinant of the matrix must vanish: U(B2 +A, 2iB)U ′(B2 +A,−2iB) + U(B2 +A,−2iB)U ′(B2 +A, 2iB) = 0 . (3.69) Now the Wronskian of U(a, x) and U(a,−x) is given by (U(a, x))′ U(a,−x)− U(a, x)(U(a,−x))′ = U ′(a, x)U(a,−x) + U(a, x)U ′(a,−x) , (3.70) where the prime here denotes the derivative with respect to the second argument; therefore (3.69) demands that the Wronskian of U(B2 + A, x) and U ′(B2 + A,−x) must vanish at x = 2iB. Since both U(B2 +A, x) and U(B2 +A,−x) are solutions to the parabolic cylinder equation d2U dx2 − ( 1 4 x2 +B2 +A ) U = 0 , (3.71) the Wronskian must be a constant due to Abel’s identity1. Therefore if the Wron- skian of U(B2 +A, x) and U ′(B2 +A,−x) were to vanish at x = 2iB, it must vanish at x = 0 as well. In other words, U(B2 +A, 0)U ′(B2 +A, 0) = 0 . (3.72) With the formulae for U and U ′ at the origin from Bender & Orszag (1999), as well as noting that that the Gamma function Γ(z) diverges when z is zero or a negative 1Abel’s identity states that for a linear second order differential equation y′′ + p(x)y′ + q(x)y = 0, the Wronskian is given by C exp (− ∫ p(x)dx). For the parabolic cylinder equa- tion p(x) = 0, and hence the Wronskian is a constant. 45 integer, we have: U(B2 +A, 0) = 0 ⇐⇒ 3 4 + 1 2 (B2 +A) = 0 ,−1 ,−2 ,−3 · · · , U ′(B2 +A, 0) = 0 ⇐⇒ 1 4 + 1 2 (B2 +A) = 0 ,−1 ,−2 ,−3 · · · , where N and Z denotes the set of natural numbers and the set of integers respectively. Together these two conditions imply B2+A = − (−12 + n), where n ∈ N. Reverting back to the original variables, we find that the correction to the eigenvalue is: γ1 = − ( 2n− 1 2 )√ 2 υ − l23/2 4 (3.73) The prediction for the eigenvalue is therefore: γ = υ + γ1 = υ − (2n− 1) √ µ 2υ − l 2 4µ (3.74) Therefore the fastest growing mode corresponds to n = 1, which in turn gives U ′(B2 +A, 0) = U ′(−1/2, 0) = 0. To complete our analysis, we will need to determine the constant r in the inner solution. Note that R̂1(Y ) = (F1 + F2)/2, and with the aid of (3.68) we can eliminate c1, leaving one undetermined coefficient for R̂1. In addition, we have not used the conditions R̂1(0) = 4 √ 8υ and R̂′1(0) = 2υr in (3.65) during our matching; after some algebra, these conditions will give: r = i 4 √ 8υ 2υ U ′(−1/2,−2iB) U(−1/2, 2iB) . (3.75) Notice that U(−1/2, x) = e−x2/4 satisfies (3.71) (Temme, 2010), and therefore r = − l3/2√ 2υ . (3.76) 46 3.2.3 Energetics To explain the unexpected results observed in the previous section, we will examine the energetics associated with the normal modes. If we substitute the normal mode solutions Aj(y)eγ−i(ω−a))t into the energy equations (3.4) and (3.5), we will have: 2γ|A1|2 = −J1y −D1 + S , (3.77) 2γ|A2|2 = −J2y −D2 + S . (3.78) Integrating over the entire domain, the equations become: γ = ( −D̂1 + Ŝ )/ 2|Â1|2 = ( −D̂2 + Ŝ )/ 2|Â2|2 , (3.79) where the hatted variables are the integrated values over the entire domain. In addition we have assumed that Jn vanishes as y → ±∞ since |An| → 0. The energy equation for the normal mode problem simply states that the growth rate of the instability is determined by the difference between energy source due to PSI activity (Ŝ) and energy sink due to dissipation (D̂1). Figure 3.6 shows the results for υ = 1 and l = 0, where total energy source Ŝ and total dissipation D̂1 are calculated using the most unstable eigenfuctions, and are normalized by the wave- energy |Â1|2. Not surprisingly, the dissipative term D̂1 (green dots) increases as the dissipation increases; at the same time, the energy source Ŝ (blue dots) also decreases significantly (20 % compared to the inviscid case), indicating that dissipation also disrupts the energy transfer. At higher dissipation a significant portion of the energy (≈ 50%) generated via PSI is dissipated away. The results for l = 0.1 is plotted in figure 3.7. For larger dissipation, the situation is very similar to the l = 0 case: The total dissipation D̂1 increases signif- icantly while the energy source Ŝ drops. On the other hand, as µ → 0 the system 47 10−4 10−3 10−2 10−1 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 µ Ŝ/2 |̂A1| 2 D̂/2 |̂A1| 2 γR Prediction Figure 3.6: Energetics of PSI: l = 0. Ŝ (blue dots) and D̂1 (green dots) from (3.79) are calculated and normalized by using the eigenfunctions, and are plotted here as a function of µ. The numerically computed growth rate (red dots) and the asymptotic prediction from (3.74) (dashed line) are also plotted for comparison. υ = 1 and yL = −200. 48 10−3 10−2 10−1 100 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 µ A B C Ŝ/2Â1 D̂/2Â1 γR Prediction Figure 3.7: Energetics of PSI: l = 0.1. Ŝ (blue dots) and D̂1 (green dots) from (3.79) are calculated and normalized by using the eigenfunctions, and are plotted here as a function of µ. The numerically computed growth rate (red dots) and the asymptotic prediction from (3.74) (dashed line) are also plotted for comparison. The vertical line separates the boundary modes and the localized modes. υ = 1 and yL = −1500. behaves very differently compared to the l = 0 case: when µ drops below 10−2, the energy source Ŝ drops off rapidly while the energy dissipation D̂1 increases. These two effects combine to stabilize the system when the dissipation is sufficiently small. To further understand the stability problem in the small dissipation limit, we will now examine the energetics more closely. If we integrate (3.77) from y to ∞, ∫ ∞ y S(s)ds = ∫ ∞ y D1(s)ds+ 2γ ∫ ∞ y |A1(s)|2ds− J1(y) , (3.80) 49 (3.80) states that total energy generated via PSI to the north of y is either dissi- pated, converted into PSI wave, or transported out of the region between y and ∞. (which is represented by J1(y). We will consider the three modes that are close to the stability boundary (γR = 0) with dissipation µ1 = 0.002238, µ2 = 0.002512, and µ3 = 0.002818; these modes correspond to the growth rates marked by A, B, and C respectively in figure 3.7. The four terms in (3.80) are computed using the eigenfunctions (normalized such that max(|A1|) = 1) and are plotted in 3.8 together with the eigenfunction. From the plot of |A1|2 in panel a, we can see clearly that the waves become more and more de-localized as µ decreases; the eigenfunctions are very similar near the turning latitude, but the eigenfunction for µ3 (blue curve) decays significantly faster than the eigenfunction for µ1 (red curve). A plot of ∫∞ y |A1(s)|2ds (panel b) illustrates the delocalization more precisely: the energy density integral for µ3 becomes constant when y ≈ −120, while for µ1 it converges only when y < −280. Earlier we saw that the dissipation increases as µ decreases, and the plot of energy flux J1 in panel c provides a possible explanation: as µ decreases, the merio- dional energy flux becomes more negative, meaning that more energy is transported southwards by the PSI waves. As the waves become increasingly oscillatory to the south, dissipation is more effective since it increases with the wave-number. Thus in this regime the energy generated near the critical latitude is transported to the south where it is dissipated away, and the end result is that the global dissipation increases despite a decrease in the dissipation parameter. The other phenomenon we have observed earlier is that the overall energy source via PSI Ŝ also decreases with µ. ∫∞ y S(s)ds is plotted in panel (e) of figure 3.8, which demonstrates how the source term changes over the domain. For µ3 50 (blue curve), it first increases rapidly as we move south and away from the turning point, indicating there is significant energy input near the critical latitude. For µ2 (green curve), the integral increases much slower near the critical latitude compared to the previous case, and therefore energy input via PSI is greatly reduced. When µ further decreases to µ3, the integral becomes negative near the turning latitude, showing that the energy is in fact being drained from the PSI waves into the pump wave. The integral continues to decrease as we head south until y = −50, where a minimum is reached and the integral starts to increase; thus the source term is positive to the south of y = −50, and energy is transfered from the pump-wave into the PSI waves. The net transfer of energy Ŝ is small as the drain of energy to the north is balanced by the gain of energy to the south. 51 −400 −300 −200 −100 0 0 5 10 15 20 25 ∫ ∞ y| A 1 |2 d y (b) −400 −300 −200 −100 0 −1.5 −1 −0.5 0 0.5 J 1 (c) −400 −300 −200 −100 0 0 2 4 6 8 ∫ ∞ yD 1 d y y (d) −400 −300 −200 −100 0 −2 −1 0 1 2 3 ∫ ∞ yS d y y (e) −100 −90 −80 −70 −60 −50 −40 −30 −20 −10 0 0 0.2 0.4 0.6 0.8 1 |A 1 |2 y (a)µ = 0.002238 µ = 0.002512 µ = 0.002818 Figure 3.8: Energetics of PSI near the stability boundary. The eigenfunctions for µ = 0.002238 (red), µ = 0.002512 (green), and µ = 0.002818 (blue) are used to compute the various terms in the energy equation (3.80). (a) |A1|2. (b) ∫∞ y |A1(s)|2ds. (c) J1. (d) ∫∞ y D1(s)ds. (e) ∫∞ y S(s)ds. υ = 1 and yL = −1500. Please note that the range for the horizontal axis in (a) is −100 < y < 5 but the range for (b) to (e) is −400 < y < 5. 52 3.3 More initial value problems A couple of Initial value problems will now be presented to further illustrate some of the properties of our PSI model discussed in §3.2. 3.3.1 Transient growth Figure 3.9 shows the result of an initial value problem where υ = 2, l = 0.5, µ = 0.028. The computational domain is [−200, 50], with a PML located between [−200, 185]. The normal mode analysis shows that the system is asymptotically stable, as the eigenvalue calculation shows the maximum growth rate to be γR = −0.3267. The initial perturbation is in the form of a Gaussian wave packet: A1(y, 0) = e−(y+25) 2/5 cos( √ 26y), A2(y, 0) = 0 (3.81) In the top panel of figure 3.9, |A1| is plotted for 0 ≤ t ≤ 2.5 and −80 ≤ y ≤ 5. As in the example shown in §3.1, the wave packet splits into two components that travel in opposite directions; between the wave packets, the wave amplitude experiences a growth due to PSI. This growth continues exponentially, and the wave amplitude increases from O(1) to O(105) by t = 20 (see the second panel); the growth comes to an end at t ≈ 20, and the wave amplitude subsequently decays exponentially at a rate that is consistent with the eigenvalue analysis. In figure 3.10 the wave amplitude |A1| at various times are plotted, which illustrates the transient growth more clearly. The numerical result suggests that a very significant transient growth is possible, even though the system is asymptotically stable. We will be examining transient growth in much greater detail in §4.3.2. 53 3.3.2 Boundary modes The second example is one where a boundary mode eventually dominates the nu- merical solution. In figure 3.11, the amplitude of the wave |A1| is plotted at various times for υ = 2, l = 0.1 and µ = 10−4; the computational domain is similar to the previous example. Panel (a) shows the initial perturbation given by a Gaussian wave packet centred at y = −75: A1(y, 0) = e−(y+75) 2/5 cos( √ 75y), A2(y, 0) = 0 . (3.82) The initial dynamics is similar to the previous example, where an exponentially growing mode is excited between the wave packets (panels (b) to panel (f)); the difference however is that the mode eventually gives way to an exponentially growing boundary mode that dominates the solution (panels (g) to (l)). Note that the wave profile appears to be solid due to the rapid oscillations, which are well resolved by the grid used in the numerical computation. As in the normal mode anlysis in §3.2.1, the type of boundary conditions (zero dirichlet, sponge layer or radiative) had minimal effect on the growth rate of the unstable boundary modes in the IVP. 54 Figure 3.9: An initial value problem with transient growth The results from a nu- merical computation of the PSI model. |A1| is plotted in the top two panel, with the top panel showing the initial growth from t = 0 to 2.5. The second panel shows the entire computation from t = 0 to 40. Note that the color scale is different for the two plots. A similar set of plots for |A2| are produced in the bottom two panels. The parameters used are υ = 2, l = 0.5, and µ = 0.028, with the spatial domain being [−185, 50] and a sponge layer located at [−200,−85]. Only the data from [−80, 5] is shown here. 55 −80 −60 −40 −20 0 20 0 0.5 1 t=0(a) −80 −60 −40 −20 0 20 0 0.1 0.2 0.3 0.4 t=0.75(b) −80 −60 −40 −20 0 20 0 0.2 0.4 t=1.5(c) −80 −60 −40 −20 0 20 0 0.5 1 t=2.25(d) −80 −60 −40 −20 0 20 0 5 10 t=4(e) −80 −60 −40 −20 0 20 0 50 100 150 200 t=6(f) −80 −60 −40 −20 0 20 0 500 1000 1500 t=7.5(g) −80 −60 −40 −20 0 20 0 1 2 3 x 104 t=10(h) −80 −60 −40 −20 0 20 0 2 4 6 8 x 105 t=17.5(i) −80 −60 −40 −20 0 20 0 1 2 3 4 x 105 t=22.5(j) −80 −60 −40 −20 0 20 0 5 10 x 104 t=27.5(k) −80 −60 −40 −20 0 20 0 2000 4000 6000 8000 t=35(l) Figure 3.10: An initial value problem with stable normal modes. |A1| at various times are plotted. The parameters used are υ = 2, l = 0.1, and µ = 0.028, with the spatial domain being [−185, 50] and a sponge layer located at [−200,−85]. Only the data from [−80, 30] is shown here. 56 −200 −150 −100 −50 0 0 0.5 1 t=0(a) −200 −150 −100 −50 0 0 0.2 0.4 0.6 t=0.75(b) −200 −150 −100 −50 0 0 0.2 0.4 0.6 0.8 t=1.5(c) −200 −150 −100 −50 0 0 1 2 3 t=2.25(d) −200 −150 −100 −50 0 0 200 400 600 t=5(e) −200 −150 −100 −50 0 0 2 4 6 8 x 108 t=12.5(f) −200 −150 −100 −50 0 0 2 4 6 8 x 1012 t=17.5(g) −200 −150 −100 −50 0 0 1 2 3 4 x 1016 t=22.5(h) −200 −150 −100 −50 0 0 1 2 3 x 1018 t=25(i) −200 −150 −100 −50 0 0 1 2 x 1020 t=27.5(j) −200 −150 −100 −50 0 0 0.5 1 1.5 2 x 1022 t=30(k) −200 −150 −100 −50 0 0 0.5 1 1.5 2 x 1026 t=35(l) Figure 3.11: An initial value problem with boundary mode. |A1| at various times are plotted here. The parameters used are υ = 2, l = 0.1, and µ = 10−4, with the spatial domain being [−185, 50] and a sponge layer located at [−200,−85]. Only the data from [−200, 10] is shown. 57 3.4 Summary In this chapter we have presented an initial value problem in which a perturbation in form of a Gaussian wave packet results in exponential growth, and a normal mode analysis confirms the existence of unstable modes. The stability problem was investigated numerically, as well as analytically in the limit of small dissipation and small pump wavenumber (i.e. l 1 and µ 1). When l is non-zero, the unexpected result is a vanishing dissipation has a stabilizing effect, which is confirmed by an initial value problem; however the IVP also showed that a significant transient growth is possible even when the PSI-waves are asymptotically stable. In the next chapter we will explore the transient growth, as well as the stability of the PSI-model in the inviscid limit. 58 Chapter 4 PSI in the Inviscid Limit The goal of this chapter is to further understand the PSI model in the inviscid limit. There are two questions we would like to address: 1. What is the shape of the spectrum for the inviscid problem? In the previous chapter, as we decreased the dissipation parameter (µ), the spectra eventually became contaminated by unstable boundary modes that have little physical meaning. While the analytical and numerical results suggests that the discrete unstable modes will eventually become stable as µ→ 0, it is not clear if there are other physical modes present. 2. What is the nature of the large transient growth observed in some initial value problems? We are interested in understanding the mathematical explanation behind the transient growth in an asymptotically stable system. 59 4.1 Stability of the inviscid problem To address the first question, we will solve the inviscid problem analytically via a Fourier transform. Setting µ in (3.18) and (3.19) to zero, the inviscid equations are: A1yy = (y − iγ)A1 + υeilyA2, (4.1) A2yy = (y + iγ)A2 + υe−ilyA1. (4.2) Applying Fourier transform F{·} to the equations, and defining F (q) = F{A1(y)e−ily/2} and G(q) = F{A2(y)eily/2}, we have dF dq = i ( q2 + 14 l 2 + ql − iγ)F + iυG , (4.3) dG dq = iυF + i ( q2 + 14 l 2 − ql + iγ)G . (4.4) Note that under Fourier transform, y is transformed into i ddq . These equations can be further simplified by a change of variable: F̂ (q) = ei “ 1 3 q 3+ 1 4 l 2q ” F (q) and Ĝ(q) = ei “ 1 3 q 3+ 1 4 l 2q ” G(q) . (4.5) (4.3) and (4.4) become: dF̂ dq = i(ql − iγ)F̂ + iυĜ , (4.6) dĜ dq = iυF̂ − i(ql − iγ)Ĝ . (4.7) We can eliminate one of the equations to form a second order ODE. For example, using (4.7) to eliminate Ĝ in (4.6), we have: d2F̂ dq2 + [ (lq − iγ)2 + υ2 − il] F̂ = 0 . (4.8) Defining z ≡ √ il 2 ( q − iγl ) and H(z) = F̂ (q) the equation can be manipulated into the standard parabolic cylinder equation: d2H dq2 − [ 1 4 z2 + ( − iυ 2 2l − 1 2 )] H = 0 . (4.9) 60 In this case the solutions are U(a, z) and U(a,−z), where a = − iυ22l − 12 ; their Wronskian is given by (Temme, 2010): W = √ 2pi Γ(12 + a) = √ 2pi Γ(− iυ22l ) , (4.10) which will not vanish since the gamma function is finite along the imaginary axis. Thus these two solutions to (4.14) are linearly independent, and form a satisfactory basis.To carryout the inverse Fourier transform, H must remain bounded as q → ±∞; however we can demonstrate that this cannot be true if γR 6= 0. We will first consider q → ∞. In this case arg z ≈ pi/4, and the asymptotic formula for |z| 1 is given by (Temme, 2010): U(a, z) ∼ e− 14 z2z−a−1/2 +O(1/z2) . (4.11) Now since z2 = il 2 ( q2 − 2 iγ l q − γ 2 l2 ) = il 2 ( q2 − γ 2 l2 ) + γq , (4.12) e−z2/4 will be bounded as q → ∞ if and only if γR ≥ 0. On the other hand, as q → −∞, arg z ≈ −3pi/4, and the appropriate asymptotic formula is U(a, z) ∼ e− 14 z2z−a− 12 − i √ 2pi Γ(12 − a) eipiae 1 4 z2za− 1 2 +O(1/z2) . (4.13) In this case the first term contains e−z2/4, which will increase without bounds unless γR ≤ 0. Combining the two constraints, we arrive at the conclusion that γR = 0 for U(a, z) to be bounded. For the other solution U(a,−z), a similar analysis will lead to the same conclusion. Therefore for if H is to remain bounded as q → ±∞, γR = 0. What the analysis above suggests is that for the inviscid problem, there are no unstable normal modes. On the other hand, both U(a, z) and U(a,−z) remains bounded whenever γ is purely imaginary. Thus a continuous spectrum of eigenvalues exists, and the corresponding eigenfunction can be obtained by reversing the Fourier transform. 61 4.2 Non-normal operators and transient growth Traditionally in fluid mechanics, the hydrodynamical stability is considered through normal mode analysis, where infinitesimally small perturbations can grow exponen- tially into finite amplitude disturbances for unstable flows. The growth rates for the normal modes are given by the eigenvalues, while the spatial structure is de- scribed by the eigenfunctions. Normal mode analysis have been quite successful for problems such as the Rayleigh-Bèrnard instability and Taylor-Couette instability, as experimental and analytical prediction for instability agree (Trefethen et al., 1993). However for shear flows such as plane Couette and Poiseuille flows, the onset of turbulence occurs at Reynolds number much lower than the threshold predicted by theory (Trefethen et al., 1993; Butler & Farrell, 1992). For example, Couette flows are theoretically stable for all Reynolds numbers, but transition to turbulence was observed in laboratory at Reynolds number of ≈ 360 (Tillmark & Alfredsson, 1992). For plane Poiseuille flow, the critical Reynolds number unstable normal modes for the Orr-Sommerfeld equation is Re > 5772.22 (Orszag, 1971), while transition to full turbulence was observed for Re as low as 1000 (Carlson et al., 1982). Since the 1990’s, the limitation of the normal mode analysis was gradually recognized. It turns out that for a problem where the operator is ‘normal’ (i.e. the eigenfunctions are orthogonal), the behaviour of the system is completely determined by the spectra. The linear operators in the Rayleigh-Bèrnard and Taylor-Couette stability problems fall into this category, and therefore the normal mode analysis succeeded in predicting the onset of instability (Trefethen et al., 1993). In contrast, the Orr-Sommerfeld operator for shear instabilities is in fact ‘non-normal’, in the sense that the eigenfunctions are not orthogonal. A consequence of non-orthogonal eigenfunctions is that while the long term behavior of a system is dictated by the 62 spectra, these systems can have significant transient growth, which the normal mode analysis fails to capture. 4.2.1 A simple example The following example was used by Farrell & Ioannou (1996) to illustrate some of the features of non-normal operators. Consider a simple 2× 2 matrix: A = −1 − cot θ 0 −2 , (4.14) A matrix M is defined to be normal if MM † = M †M , where † denotes the Hermitian transpose1. Here we have AA† −A†A = cot2 θ cot θ cot2 θ 0 , (4.15) and therefore A is normal only if θ = pi/2, and becomes increasingly non-normal as θ → 0. It turns out that the non-normality can play an important role in the transient behaviour of a dynamical system. Farrell & Ioannou (1996) considered the following dynamical system: dx dt = Ax . (4.16) The matrix A has eigenvalues r = −1 and r = −2 with eigenvectors v1 = (1, 0) and v2 = (cos θ, sin θ), and the general solution is given by: x = c1e−t 1 0 + c2e−2t cos θ sin θ . (4.17) Consequently all solutions are asymptotically stable (i.e. x → 0 as t → ∞ for all initial conditions x0), with decay rate −1. Farrell & Ioannou (1996) continued on 1A Hermitian transpose is also known as the conjugate transpose 63 to demonstrate that an asymptotic stability analysis does not capture the dynamics of the system as θ → 0. In this limit, v2 → v1, and the eigenvectors become almost parallel to each other. The ramification is that if the initial condition is almost orthogonal to these eigenvectors, c1 and c2 have to be large and have opposite signs in order for the initial condition to be satisfied. For example, the solution with initial condition x0 = (sin θ,− cos θ), which is orthogonal to v2, will be given by c1 = sin θ + cos θ cot θ = csc θ and c2 = − cot θ; clearly in this case as θ → 0, c2 decreases without bound, and c1 → −c2. The solution to the initial value problem is x = e−t csc θ 0 + e−2t sin θ − csc θ − cos θ . (4.18) Now as t increases, the second term will quickly vanish, while the first term can have a magnitude much greater than unity if θ is small. As |x0| = 1, this indicates that the amplitude experiences a significant growth before decaying. To further illustrate this, notice that the solution to (4.16) with initial condition x(0) = x0 can be written as: x = eAtx0 , (4.19) If we restrict ourselves to initial conditions that satisfies |x0| = 1, then max |x| = max |x0|=1 |eAtx0| = ||eAt|| , (4.20) where || · || denotes the induced matrix 2-norm. ||eAt|| is therefore a measure of the maximum magnitude of x. In figure 4.1, ||eAt|| for θ = pi/2,pi/4,pi/10, and pi/100 are plotted as a function of time. For all four cases, the curves all become a straight line with slope -1 in accordance with the asymptotic stability analysis. The main difference is that for θ = pi/2, which corresponds to A being a normal matrix, the solution decays with a decay rate of -1 for all t. On the other hand, as θ decreases 64 0 0.5 1 1.5 2 2.5 3 3.5 4 10−1 100 101 t ||e A t || θ = pi/2 θ = pi/4 θ = pi/10 θ = pi/100 Slope = −2 Slop e = 12 [−3 + csc(pi/100) ] Figure 4.1: Transient growth for a simple dynamical system. ||eAt|| is plotted for four values of θ: θ = pi/2,pi/4,pi/10, and pi/100, where A is defined in (4.14). Also plotted as dash line is the initial growth rate give by 12(−3+csc θ), and the asymptotic decay rate -1 for θ = pi/2. Note that the vertical axis is in logarithmic scale. to pi/4 and pi/10, A becomes increasingly non-normal, and the decay in ||eAt|| is delayed. For pi/100, the difference is even more remarkable, as ||eAt|| increased almost an order of magnitude before decaying. The initial growth rate near t = 0 is given by the maximum eigenvalue of (A+A†)/2, which in this case is 12(−3 + csc θ) (Farrell & Ioannou, 1996). This suggests that for positive initial growth the criterion is θ < pi/9.244. From figure 4.1, we can see that for θ = pi/10 < pi/9.244, ||eAt|| indeed increases slightly near t = 0. For θ = pi/100, the curve for ||eAt|| agrees well with the predicted initial growth rate of 12(−3 + cscpi/100). 65 4.2.2 Non-normal operators and fluid dynamics This simple example in the previous section illustrates that while the dynamics of a normal system can be characterized by the eigenvalues, they do not capture the transient behaviour for non-normal systems. This fact was demonstrated clearly in a study by Reddy et al. (1993) where the Orr-Somerfeld operator in the Poiseuille problem was analyzed; for a Reynolds number in the stable regime, the energy increased by a factor of 51 before the perturbation decayed away. Similarly, Butler & Farrell (1992) investigated the optimal excitation of perturbations in a plane Couette flow, and found energy amplifications of O(1000) despite Couette flow being stable. The implication of these findings is that transient growth can amplify small perturbations and trigger the transition to turbulence without the need to normal mode instability. Non-normal operators also received attention in the geophysical fluid dynamics community in recent years, for example in the studies of baroclinic instability, error amplification in weather forecast, El Niño, and the stability of thermohaline circulation (Farrell & Ioannou, 1996; Trefethen & Embree, 2005, Ch. 23, and references therein). 4.3 Pseudospectra and PSI 4.3.1 Pseudospectra Properties of a non-normal operator/matrix can be characterized by the pseudospec- tra. For ≥ 0, the -pseudospectra of matrix A, σ(A), is defined as the set of z ∈ C such that z is an eigenvalue of A+E for some perturbation matrix ||E|| ≤ (Reddy et al., 1993). For normal operators, a perturbation to the matrix A of size ||E|| ≤ will perturb the spectra of A, σ(A), by a distance no more than . On the other 66 hand, for a non-normal matrix the eigenvalues can be perturbed much larger than the perturbation to the matrix (Reddy et al., 1993). The pseudospectra has to be computed numerically for most cases, and it provides a convenient way to visualize the degree of non-normality of an matrix operator. One important information that can be extracted from the pseudospectra is the lower bound to the maximum of ||eAt||. Suppose that α(A) is the supremum of the real part of σ(A), then the following inequality holds: sup t≥0 ||eAt|| ≥ α(A)/ . (4.21) Since (4.2) holds for all ≥ 0, a good estimate for the lower bound is the Kreiss constant, which is defined as K(A) = sup ≥0 α(A)/ , (4.22) which leads to sup t≥0 ||eAt|| ≥ K(A) . (4.23) For details, please see Reddy et al. (1993) or Chapters 14 ad 15 in Trefethen & Embree (2005). 4.3.2 Transient growth We will now consider our stability problem in the context of pseudospectra. Re- writing (3.1) and (3.2), A1t = −i [ −(1− iµ)A1yy + yA1 + υeilyA2 ] , (4.24) A2t = −i [ (1 + iµ)A2yy − yA2 − υe−ilyA1 ] . (4.25) (4.24) and (4.25) can be written as a dynamical system: d ~A dt = −iM ~A , (4.26) 67 where M is the matrix resulting from a spatial discretization. In this section we are going to discretized the spatial derivative using a spectral collocation (or pseudo- spectral) method using Chebyshev polynomial as interpolant (see Weideman & Reddy (2000) or Trefethen (2000) for more details). While the resulting matrix is dense compared to the finite difference scheme used in chapter 3, its smaller size allow for the use of eig in MATLAB to compute the full set of eigenvalues; in addition it results in much better convergence in the calculation of the pseudospec- tra. We have compared the spectra for the same set of parameters using the same method, and they agree over the region of interest. In §3.3.1, we have seen that for υ = 2, l = 0.5, and µ = 0.028, the PSI waves grew significantly before decaying as predicted by the normal mode analysis, and we would like to investigate this transient growth with the aid of the pseudospectra as- sociated with the spatial operator. The pseudospectra for the corresponding matrix M is calculated using the eigtool package2 for MATLAB. Due to the extra factor of −i in (4.26), we will have to consider the imaginary part of the spectrum and pseudospectra for information pertaining to growth. In figure 4.2, we can see that all eigenvalues (black dots) lie below the real axis, and therefore the normal modes are stable, confirming the long term behaviour exhibited by the IVP in §3.3. On the other hand, the pseudospectra protrudes significantly above the real-axis, indi- cating that the discretized matrix M is highly non-normal. By considering σ(M) for = 10−3, we can see that a perturbation of O(10−3) results in an O(1) shift in the eigenvalues. Since the operator is non-normal, the transient growth we have observed in the IVP should not be surprising. We will now take a closer look at the transient growth in the corresponding 2Thomas G. Wright. EigTool. http://www.comlab.ox.ac.uk/pseudospectra/ eigtool/, 2002. 68 γI γ R −10 −8 −6 −4 −2 0 2 4 6 8 10 −2 −1.5 −1 −0.5 0 0.5 1 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 Figure 4.2: Pseudospectra for a stable system. The -pseudospectra for the dis- cretized spatial operator M for (4.24) and (4.25) is plotted (from top to bottom) as level curves for = 10−2, 10−3, . . . 10−14. The eigenvalues are plotted as black dots. IVP by considering the total energy contained in the PSI waves: Etot(t) = ∫ ∞ −∞ (|A1|+ |A2|) dy (4.27) In figure 4.3, Etot (normalized so that the Etot(0) = 1) is plotted as a function of time (blue curve). We can see that the energy increased to almost 1013 times the initial energy before the decay sets in at a rate predicted by the eigenvalues (dashed black line). The magnitude of growth is consistent with the lower bound given by (4.21): from figure 4.2, the -pseudospectra for = 10−7 have a maximum imaginary part of 0.187, and therefore the lower bound for the maximum growth of ||eAt|| will 69 be 1.87×106, which is equivalent to an increase in energy by a factor of ≈ 3.5×1012. This estimation is plotted in figure 4.3 as the blue dotted line, which shows that the observed amplification is in line with the prediction given by the pseudo spectra. In figure 4.3, Etot for µ = 0.033 (green curve) and µ = 0.038 (red curve) are also plotted for comparison. The dissipation parameter µ appears to only have a minor effect on the transient growth rate; the onset of transient growth is delayed slightly for µ = 0.038, but in all three cases the transient growth is exponential with comparable growth rates. The transient growth continues pass t > 15, before the growth/decay predicted by the normal mode analysis sets in. The growth rate for Etot in the transient phase is estimated to be ≈ 2.6, whereas the asymptotic growth rates are -0.65, 0.015, and 0.50 respectively for µ = 0.028, 0.033, and 0.038. This indicates that the transient growth plays a significant role in the evolution of the PSI-waves, particularly in cases where the normal modes are neutral or stable; even when unstable normal modes are present, the transient growth can still occur at a significantly shorter time scale compared to the asymptotic growth. We will present a plot to illustrate the dependence of transient growth on the initial condition. In figure 4.4, we have plotted Etot for five different initial configurations (see table 4.1). We have not explored in depth the optimal excitation that leads to the maximum transient growth, but from the figure we can notice that the wave-number of the perturbation affects the transient growth significantly. The physical significance of our result is that transient growth, rather than unstable normal modes, appears to be the dominant force in the excitation of PSI- waves in certain regimes; even in cases where the PSI-waves are asymptotically stable, the transient growth may be sufficiently large for a transition into a non- 70 0 5 10 15 20 25 30 35 40 100 105 1010 1015 1020 t E t o t µ= 0.028 µ= 0.033 µ= 0.038 Figure 4.3: Transient growth in a time dependent problem. The total energy for the IVPs considered in §?? are plotted as a function of time. The black dashed lines on the right are the asymptotic growth/ decay predicted by the eigenvalues, while the blue dotted line is the minimum growth predicted by the pseduospectra. linear regime. 71 0 2 4 6 8 10 12 14 16 18 20 10−2 100 102 104 106 108 1010 1012 1014 t E t ot IC1 IC2 IC3 IC4 IC5 Figure 4.4: Initial condition and transient growth. The total energy Etot for different initial conditions are plotted as a function of time. υ = 2, l = 0.5, and µ = 0.028. Table 4.1: Initial conditions for figure 4.4 A1(y, 0) A2(y, 0) IC1 e−(y+25)2/5 cos( √ 5y) 0 IC2 e−(y+25)2/5 cos( √ 5y) ie−(y+25)2/5 cos( √ 5y) IC3 e−(y+25)2/5 cos( √ 40y) 0 IC4 e−(y+25)2/5 cos( √ 15y) 0 IC5 e−(y+25)2/5 cos( √ 25y) 0 72 4.3.3 Numerical difficulties The non-normality in operators/matrices also results in complications during nu- merical computations of eigenvalues, since the operator is sensitive to errors; the round off errors act as a perturbation to the matrix, and if the operator is non- normal, the perturbation is sufficient to shift the eigenvalues. An example of the round-off errors can be detected in figure 4.2; near the center of the plot where the three branches meet, the eigenvalues are not symmetric about the imaginary axis, and they are sensitive to slight changes in the domain length. Reddy et al. (1993) considered the Airy operator i Re d2 dy2 + y , (4.28) as a model for the Orr-Sommerfeld operator in shear-flow instabilities. The spectra of the Airy operator consists of three branches in a ‘Y’-shape, which resembles the spectrum in figure 4.2. Reddy et al. found the operator to be non-normal, and the eigenvalues near the intersection were difficult to compute as they were sensitive to perturbations; the sensitivity increased with as the Reynolds number Re increased. Notice that the Airy operator is in fact part of the spatial operator in our PSI-model, and Re → ∞ in Reddy et al. corresponds to µ → 0 in our model. Thus in hind sight, we should not be surprised by the non-normality of our eigenvalue problem in the inviscid limit. The effect of round-off error is much more evident in figure 4.5, where the eigenvalues for discretized operator with υ = 6, l = 1.5 and µ = 0.1 are computed using two different level of precision. In the right panel, the eigenvalues are calcu- lated using double precision floating point, and we can see that the eigenvalues line up along the boundary of the pseduospectra with = 10−14. In the left panel, the 73 eigenvalues for the same matrix were computed with single precision; the eigenvalues near the center now lie along the pseudospectra for = 10−6. At the same time, the eigenvalues in the lower left and right regions of the spectra are very similar for the two computations, showing that a change in numerical precision did not affect the accuracy significantly for these parts of the spectra. When the non-normality becomes severe, the round-off error prevents the eigenvalues near the center from being calculated accurately. The growth rate pre- dicted by our asymptotic theory in §3.2.2 (c.f. (3.74)) for υ = 6, l = 1.5 and µ = 0.1 is 0.346. In panel (b) of figure 4.5, the predicted eigenvalue of 0.346i is right in the middle of the area enclosed by the eigenvalues perturbed due to round-off error. Due to the severe non-normality, the true eigenvalue cannot be resolved even by the eigs command in MATLAB. This effect of non-normality is the most prominent in the transition region between the boundary modes in the inviscid limit to the localized modes. The example considered above is marked in panel (b) of figure 3.5 by A. As the largest γR reported by the eigenvalue calculation near A is the perturbed eigenvalue due to the round-off error rather than the true eigenvalue, the growth rate is potentially overestimated. This explains why the numerically calculated growth rate agrees with the prediction for larger values of µ, but they break off suddenly from the prediction when µ becomes too small. By computing the pseudospectra near the breaking point for decreasing values of µ, we have confirmed that the round off error starts affecting the eigenvalue with the largest growth rate precisely at the breaking point in the plot. The non-normality of the matrix also affected the use of iterative schemes, such as the eigs function in MATLAB, for the computation of eigenvalues. In 74 −20 −15 −10 −5 0 5 10 15 20 −2 −1 0 1 2 3 4 (a) Single Precision −20 −15 −10 −5 0 5 10 15 20 −2 −1 0 1 2 3 4 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 (b) Double Precision Figure 4.5: Effect of non-normality on the calculation of eigenvalues The eigenvalues computed for υ = 6, l = 1.5 and µ = 0.1 for (a) single precision floating point, and (b) double precision floating point. particular, when the eigenvectors become almost linearly dependent, the convergence in such iterative schemes can be significantly slowed (Trefethen & Embree, 2005); this has been a problem for our calculations, especially in the invisicid limit. An additional complication is that even when the iterations converge, the estimated eigenvalues (Ritz values) are not guaranteed to be close to the actual eigenvalues (Trefethen & Embree, 2005). Indeed in the inviscid limit we encountered difficulty in using eigs to compute the eigenvalues as the convergence was slow. 4.4 Summary In this chapter we considered the stability of the inviscid PSI problem, and concluded that there are no unstable normal modes. In addition, as µ→ 0, the spatial operator in our model becomes non-normal. A consequence of the non-normality is that a significant transient growth can occur with suitable perturbations, and the transient 75 growth can in fact dominate the dynamics of PSI. Therefore in the inviscid limit, the PSI model may best be understood via its transient dynamics, rather than the traditional modal instability. 76 Chapter 5 Discussion 5.1 Horizontal dissipation In this thesis, we have formulated a linear model for PSI, and explored its behaviour in Chapter 3. Solving the PDE numerically as an initial value problem revealed the presence of exponentially growing normal modes, which suggest that the PSI problem can be formulated in a normal mode stability analysis. The main result is that for a meridionally propagating pumpwave (i.e. l 6= 0), instability can only occur when the horizontal eddy viscosity µH is sufficiently large. Through the analysis of the energy equation, the stability in the inviscid limit can be attributed to two factors. A reduction in µH increases the meridional energy flux towards the equator, which prevents a build up of energy in the PSI waves near the critical latitude. At the same time the energy transfer from the pump wave to the PSI waves is also less effective. The β-effect and the dissipation together limit the horizontal extent of PSI. As the PSI-waves are near-inertial at the critical latitude of 28.9◦, the decrease in Coriolis force to the south implies that the horizontal wave number must increase, 77 and therefore the wave energy is dissipated away more rapidly to the south. This provides an alternative mechanism to the localization of PSI seen in MacKinnon & Winters (2005), where the localization was attributed to a vertical divergence in energy flux preventing a build up of energy to the south of the critical latitude. Depending on the horizontal scale, the horizontal eddy diffusivity measured from dye-release experiments ranges from 10−2 m2s−1 for a horizontal scale of 0.1km to 103m2s−1 at a scale of 1000km (Ledwell et al., 1998). The intermediate values of 100 − 101m2s−1 correspond to µ ∼ 10−2 − 10−1 since the eddy diffusivity is scaled by λ3β = 1.36 × 102m2s−1. For this range of µ, our model predicts that the PSI activity will be localized to between 15-30 units, or 300-600km, to the south of the critical latitude (c.f. figure 3.3 for example). For comparison, the Hawaiian Islands, which is a strong site for generation of mode-1 internal tides, lies between 19.5◦N to 22◦N, or approximately between 750 to 1000 km south of the critical latitude. Thus the meridional extent of PSI imposed by horizontal eddy diffusion can be significant compared to other length scales, such as the source location for the pump-wave. For an near-intertial wave, we can demonstrate that the horizontal eddy diffusivity is significant compared to vertical eddy diffusion. For a plane wave, the horizontal and vertical eddy diffusion can respectively be characterized by µH(k2+l2) and µVm2. We wish to consider the ratio between the two quantities. Using the dispersion relation for an internal wave with the near-inertial approximation ω = f0, f = f0 + βy, and N f µH(k2 + l2) µVm2 = µH µV ω2 − f2 N2 − f2 ≈ µH µV −2f0βy N2 . (5.1) Taking µV = 10−5m2s−1, µH = 10m2s−1 and the parameters from table 2.1, the 78 last term in (5.1) is: µH(k2 + l2) µVm2 = −y 27km (5.2) What this tells us is that near the origin, the vertical dissipation is much larger than the horizontal dissipation, which is expected since near the critical latitude, the vertical structure of near-inertial waves is much finer than the horizontal structure; however as the wave packet moves south, the horizontal wave number increases, and at 27km south of the critical latitude the horizontal dissipation is the same order of magnitude as the vertical dissipation. This justifies the include of the horizontal dissipation in our model. In summary, while horizontal dissipation was included originally as a cure for numerical issues (boundary modes), it is clear that it is in fact dynamically significant for an near-inertial wave on the β-plane. The normal mode problem of our PSI model is also interesting from a math- ematical perspective. We were able to calculate the exact solution to the eigenvalue problem with µ = l = 0, where a continuous spectrum of unstable modes are possi- ble when υ > 0. On the other hand, while a continuous spectrum of eigenfunctions also exists for the inviscid problem with µ = 0 but l 6= 0, the eigenfunctions are all stable. This suggests that µ = l = 0 is a singular limit as we cannot recover the solutions by taking l→ 0. The inviscid solutions appear to be singular as well since the growth rate of the discrete normal modes in the dissipative problem tends to −∞ as µ→ 0. 79 5.2 Transient growth Another interesting mathematical property of our model is the non-normality of the spatial operator. As we approach the inviscid limit, the eigenfunctions for the operator become parallel to each other, which caused difficulties in our numerical computations for eigenvalues. A closer look at our spatial differential operator re- vealed connections to the Orr-Sommerfeld operator for the stability problem in plane Poiseuille and Couette flows, and the Airy operator, which are all non-normal. One of the ways the non-normality manifests itself in a time dependent prob- lem is through transient growth. Through a series of initial value problems, we have demonstrated that transient growth is indeed possible, despite the system being asymptotically stable. Even in cases where the system is asymptotically unstable, the transient growth rate is in fact much higher than the growth rate for the unstable normal modes. The amplitude and wavelength of the pump-wave considered in YTB corre- sponds to υ ≈ 3 and l ≈ 1 with our scaling. With 10−2 < µ < 10−1, the PSI is in the regime where the growth rate increases with dissipation (c.f. panel (a) of figure 3.3). For µ = 0.05, the normal modes are all stable, yet the pseudo-spectra indicates the possibility of large amplitude transient growth (see figure 5.1). The implication is that the transient growth we observed may be oceanographically relevant; transient growth may be the mechanism via which a transition into non-linear regime takes place, rather than via normal mode instability. 80 γI γ R −8 −6 −4 −2 0 2 4 6 8 −1 −0.5 0 0.5 1 −8 −7 −6 −5 −4 −3 −2 −1 Figure 5.1: Pseudospectra for a stable system. The -pseudospectra for υ = 3, l = 1, and µ = 0.05 are plotted as level curves for = 10−1, 10−3, . . . 10−8. The eigenvalues are plotted as black dots. The operator is discretized using a spectral collocation method using Chebyshev polynomials as interpolant over the domain [−75, 75] 5.3 Outstanding issues In this thesis we have established that dissipation plays an important role in the stability of PSI. However a number of key questions remain: 1. What are the implications of a realistically stratified ocean? Through-out this thesis we have assumed a linearly stratified ocean, and the vertical structure of the PSI wave and that of the pump wave are assumed to be independent. Further more, the pump-wave is assumed to be a plane wave. 81 A more realistic profile for the ocean is the Gill’s profile, where the ocean is modeled by a mixed layer with uniform density near the surface, and the density increases smoothly below the mixed layer (Gill, 1984). Mathematically, the buoyancy frequency is given by: N(z) = 0 if 0 ≤ z ≤ hmixN0/(z − z0) if z ≥ hmix (5.3) The internal-tides associated with this profile can be solved analytically, and these can be used in place of the planar pump-wave in our model. YTB investigated PSI with the Gill’s profile on the f -plane, and found intensification of PSI at the base of the mixed layer. On the other hand, the measurements of Alford et al. (2007) found PSI activity extending much deeper into the ocean. By generalizing our model to a realistically stratified ocean, we can perhaps better reproduce the observations. 2. How is the vertical scale of the PSI-waves determined? Another aspect of PSI that demands clarification is the vertical scale selection. The numerical results of MacKinnon & Winters (2005) and the field data of Alford et al. (2007) both suggest that PSI favours a particular vertical scale. In the current model the vertical wavenumber m is a scaling parameter, and therefore it is not clear immediately whether vertical scale selection occurs. 3. What is the optimal excitation for transient growth? For parameters that are oceanographically relevant, our analysis earlier in- dicated the potential for significant transient growth. We have also demon- strated that the transient growth depends strongly on the initial perturbation. 82 Therefore to fully understand PSI in the ocean, we must understand the rela- tionship between the perturbation and the subsequent transient growth. 4. When does non-linearity become dynamically important? Both YTB and the current study is based on the linear model of PSI, but It is clear in our model that perturbations are either transiently or asymptotically unstable, and small perturbations can be amplified to the point where a con- stant amplitude for the pump-wave is no longer valid. In addition the main motivation for the study of PSI is to understand the energy transfer out of internal tides, which requires a model with genuine three-mode coupling that is similar to (1.7) to (1.9). Our linear model essentially takes the form: i ∂A1 ∂t + propagation terms = ΥA2 , (5.4) i ∂A∗2 ∂t + propagation terms = ΥA∗1 , (5.5) which is similar to the three mode model of Bretherton (1964). In the deriva- tion of our PSI model (chapter 2), the evolution equation for the pump-wave can be derived if we allow the background variables to vary on the slow time scale. We expect the third equation to take the form of: i ∂Υ ∂t + propagation terms = A1A2 . (5.6) The fully non-linear equations will allow us to explore a wider range of prob- lems and resolve some outstanding issues. One important issue is the dis- crepancy between the catastrophic loss of energy from the internal tide in the numerical model of MacKinnon & Winters (2005), and the more moderate loss measured by Alford et al. (2007); while the results of YTB and our study sug- gests that the time-scale over which the instability occurs is oceanographically 83 relevant, one must study non-linear effects such as saturation of the instability and the feedback on the pump-wave to fully understand the PSI in the ocean. 84 Bibliography Abramowitz, Milton & Stegun, Irene A. 1964 Handbook of mathematical functions with formulas, graphs, and mathematical tables, National Bureau of Standards Applied Mathematics Series, vol. 55. For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. Alford, M. H, Mackinnon, J. A, Zhao, Zhongxiang, Pinkel, Rob, Kly- mak, Jody & Peacock, Thomas 2007 Internal waves across the pacific. Geo- physical Research Letters 34 (24), 6. Bender, Carl M. & Orszag, Steven A. 1999 Advanced mathematical methods for scientists and engineers. I . New York: Springer-Verlag, asymptotic methods and perturbation theory, Reprint of the 1978 original. Berenger, J.P 1994 A perfectly matched layer for the absorption of electromagnetic-waves. Journal of Computational Physics 114 (2), 185–200. Bretherton, F. P 1964 Resonant interactions between waves. the case of discrete oscillations. Journal of Fluid Mech 20, 457–479. Butler, KM & Farrell, B 1992 Three-dimensional optimal perturbations in viscous shear flow. Physics of Fluids A 4 (8), 1637–1650. Carlson, D. R., Widnall, S. E. & Peeters, M. F. 1982 A flow-visualization study of transition in plane Poiseuille flow. Journal of Fluid Mechanics 121, 487– 505. Carter, G & Gregg, M 2006 Persistent near-diurnal internal waves observed above a site of m 2 barotropic-to-baroclinic conversion. Journal of Physical Oceanography 36, 1136–1146. Dickey, JO, Bender, PL, Faller, JE, Newhall, XX, Ricklefs, RL, Ries, JG, Shelus, PJ, Veillet, C, Whipple, AL & Wiant, JR 1994 Lunar laser ranging- a continuing legacy of the apollo program. Science 265 (5171), 482–490. 85 Egbert, GD 1997 Tidal data inversion: Interpolation and inference. Progress in Oceanography 40 (1-4), 53–80. Egbert, GD & Ray, RD 2000 Significant dissipation of tidal energy in the deep ocean inferred from satellite altimeter data. Nature 405 (6788), 775–778. Egbert, GD & Ray, R 2001 Estimates of m 2 tidal energy dissipation from topex/poseidon altimeter data. Journal of Geophysical Research 106 (C10), 22475–22502. Farrell, B & Ioannou, P 1996 Generalized stability theory. part i: Autonomous operators. Journal of the Atmospheric Sciences 53 (14), 2025–2040. Garcia, Alejandro L. 1999 Numerical Methods for Physics (2nd Edition). Upper Saddle River, NJ, USA: Prentice-Hall, Inc. Garrett, C 2003 Internal tides and ocean mixing. Science 301 (5641), 1858–1859. Garrett, Chris & Kunze, Eric 2007 Internal tide generation in the deep ocean. Annual Review of Fluid Mechanics 39, 57–87. Garrett, C & Munk, W 1972 Space-time scales of internal waves. Geophysical & Astrophysical Fluid Dynamics 3 (1), 225. Garrett, C & Munk, W 1975 Space-time scales of internal waves: A progress report. Journal of Geophysical Researchophys. Res. 80 (3), 291–297. Gill, A. E. 1984 On the behavior of internal waves in the wakes of storms. Journal of Physical Oceanography 14 (7), 1129–1151. Hasselmann, K 1966 Feynman diagrams and interaction rules of wave-wave scat- tering processes. Reviews of Geophysics 4 (1), 1–32. Hasselmann, K 1967 Nonlinear interactions treated by the methods of theoret- ical physics (with application to the generation of waves by wind). Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 299 (1456), 77–103. Hibiya, T, Nagasawa, M & Niwa, Y 2002 Nonlinear energy transfer within the oceanic internal wave spectrum at mid and high latitudes. Journal of Geophysical Research-Oceans 107 (C11), 3207. Hibiya, T, Niwa, Y & Fujiwara, K 1998 Numerical experiments of nonlinear energy transfer within the oceanic internal wave spectrum. Journal of Geophysical Research 103 (C9), 18715–18722. 86 Jeffreys, H 1921 Tidal friction in shallow seas. Philosophical Transactions of the Royal Society of London. Series A 221, 239–264. Laurent, L St. & Garrett, C 2002 The role of internal tides in mixing the deep ocean. Journal of Physical Oceanography 32, 2882–2898. Ledwell, JR, Watson, AJ & Law, CS 1993 Evidence for slow mixing across the pycnocline from and open-ocean tracer-release experiment. Nature 364 (6439), 701–703. Ledwell, J.R, Watson, A.J & Law, C.S 1998 Mixing of a tracer in the pycno- cline. Journal of Geophysical Research-Oceans 103 (C10), 21499–21529. MacKinnon, JA & Winters, KB 2005 Subtropical catastrophe: Significant loss of low-mode tidal energy at 28.9 degrees. Geophysical Research Letters 32 (15), L15605. McComas, C & Bretherton, F 1977 Resonant interaction of oceanic internal waves. Journal of Geophysical Research-Oceans 82 (9), 1397–1412. Miller, G 1966 The flux of tidal energy out of deep oceans. Journal of Geophysical Research 71 (10), 2485–2489. Müller, P, Holloway, G, Henyey, F & Pomphrey, N 1986 Nonlinear inter- actions among internal gravity waves. Reviews of Geophysics 24 (3), 493–536. Munk, W. & Wunsch, C. 1998 Abyssal recipes ii: energetics of tidal and wind mixing. Deep-Sea Res Pt I 45 (12), 1977–2010. Olbers, D & Pomphrey, N 1981 Disqualifying two candidates for the energy balance of the oceanic internal wave field. Journal of Physical Oceanography 11, 1423–1425. Orszag, S. A. 1971 Accurate solution of the Orr-Sommerfeld stability equation. Journal of Fluid Mechanics 50, 689–703. Phillips, O. M 1960 On the dynamics of unsteady gravity waves of finite amplitude. i. the elementary interactions. J Fluid Mech 9, 193–217. Phillips, O. M 1961 On the dynamics of unsteady gravity waves of finite amplitude. ii. local properties of a random wave field. J Fluid Mech 11, 143–155. Rainville, L & Pinkel, R 2006 Baroclinic energy flux at the hawaiian ridge: Observations from the r/p flip. Journal of Physical Oceanography 36 (6), 1104– 1122. 87 Ray, R & Mitchum, G 1997 Surface manifestation of internal tides in the deep ocean: Observations from altimetry and island gauges. Progress in Oceanography 40 (1-4), 135–162. Reddy, SC, Schmid, PJ & Henningson, DS 1993 Pseudospectra of the orr- sommerfeld operator. SIAM Journal on Applied Mathematics 53 (1), 15–47. Simmons, Harper L 2008 Spectral modification and geographic redistribution of the semi-diurnal internal tide. Ocean Model 21 (3-4), 126–138. Taylor, G 1920 Tidal friction in the irish sea. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Phys- ical Character 220, 1–33. Temme, N M 2010 Digital library of mathematical functions. National Institute of Standards and Technology. Thorpe, S.A. 2005 The Turbulent Ocean. Cambridge University Press. Tillmark, N. & Alfredsson, P. H. 1992 Experiments on transition in plane Couette flow. Journal of Fluid Mechanics 235, 89–102. Trefethen, L, Trefethen, A, Reddy, S & Driscoll, T 1993 Hydrodynamic stability without eigenvalues. Science 261 (5121), 578–584. Trefethen, Lloyd N. 2000 Spectral methods in MatLab. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics. Trefethen, Lloyd N. & Embree, Mark 2005 Spectra and pseudospectra. The behavior of nonnormal matrices and operators. Weideman, J. A. & Reddy, S. C. 2000 A matlab differentiation matrix suite. ACM Trans. Math. Softw. 26 (4), 465–519. Young, W. R, Tsang, Y. K & Balmforth, N. J 2008 Near-inertial parametric subharmonic instability. J Fluid Mech 607, 25–49. Zheng, C 2007 A perfectly matched layer approach to the nonlinear schrödinger wave equations. Journal of Computational Physics 227 (1), 537–556. 88
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Parametric subharmonic instability and the β-effect
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Parametric subharmonic instability and the β-effect Chan, Ian 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Parametric subharmonic instability and the β-effect |
Creator |
Chan, Ian |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | Parametric subharmonic instability (PSI) is a nonlinear interaction between a resonant triad of waves, in which energy is transferred from low wavenumber, high frequency modes to high wavenumber, low frequency modes. In the ocean, PSI is thought to be one of the mechanisms responsible for transferring energy from M₂ internal tides (internal gravity waves with diurnal tidal frequency) to near-inertial waves (internal gravity waves with frequency equal to the local Coriolis frequency) near the latitude of 28.9 degrees. Due to their small vertical scale, near-inertial waves generate large vertical shear and are much more efficient than M₂ internal tides at generating shear instability needed for vertical mixing, which is required to maintain ocean stratification. The earlier estimate of the time-scale for the instability is an order of magnitude larger than the time-scale observed in a recent numerical simulation (MacKinnon and Winters) (MW). An analytical model was developed by (Young et al. 2008) (YTB), and their findings agreed with the MW estimation; however as YTB assumed a constant Coriolis force, the model cannot explain the intensificaiton of PSI near 28.9 degrees as observed in the model of MW; in addition, the near-ineartial waves can propagate a significant distance away from the latitude of 28.9 degrees. This thesis extends the YTB model by allowing for a linearly varying Coriolis parameter (β-effect) as well as eddy diffusion. A linear stability analysis shows that the near-inertial wave field is unstable to perturbations. We show that the β-effect results in a shortening in wave length as the near-inertial waves propagate south; horizontal eddy diffusion is therefore enhanced to the south, and limits the meridional extent of PSI. The horizontal diffusion also affects the growth rate of the instability. A surprising result is that as the horizontal diffusion vanishes, the system becomes stable; this can be demonstrated both analytically and numerically. Mathematically, the β-effect renders the spatial differential operator nonnormal, which is characterized with the aid of pseudo-spectra. Our results suggest the possibility of large amplitude transient growth in near-inertial waves in regimes that are asymptotically stable to perturbations. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0071275 |
URI | http://hdl.handle.net/2429/28200 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_fall_chan_ian.pdf [ 2.29MB ]
- Metadata
- JSON: 24-1.0071275.json
- JSON-LD: 24-1.0071275-ld.json
- RDF/XML (Pretty): 24-1.0071275-rdf.xml
- RDF/JSON: 24-1.0071275-rdf.json
- Turtle: 24-1.0071275-turtle.txt
- N-Triples: 24-1.0071275-rdf-ntriples.txt
- Original Record: 24-1.0071275-source.json
- Full Text
- 24-1.0071275-fulltext.txt
- Citation
- 24-1.0071275.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071275/manifest