UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Parametric subharmonic instability and the β-effect Chan, Ian 2010-12-31

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


ubc_2010_fall_chan_ian.pdf [ 2.29MB ]
JSON: 1.0071275.json
JSON-LD: 1.0071275+ld.json
RDF/XML (Pretty): 1.0071275.xml
RDF/JSON: 1.0071275+rdf.json
Turtle: 1.0071275+rdf-turtle.txt
N-Triples: 1.0071275+rdf-ntriples.txt
Original Record: 1.0071275 +original-record.json
Full Text

Full Text

Parametric Subharmonic Instability and the  -e ectbyIan ChanB.Sc., University of British Columbia, 2007A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE STUDIES(Mathematics)The University of British Columbia(Vancouver)August 2010c Ian Chan, 2010AbstractParametric subharmonic instability (PSI) is a nonlinear interaction between a res-onant triad of waves, in which energy is transferred from low wavenumber, highfrequency modes to high wavenumber, low frequency modes. In the ocean, PSI isthought to be one of the mechanisms responsible for transferring energy from M2internal tides (internal gravity waves with diurnal tidal frequency) to near-inertialwaves (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 wavesgenerate large vertical shear and are much more e cient than M2 internal tides atgenerating shear instability needed for vertical mixing, which is required to maintainocean strati cation.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  ndings agreed with the MW estimation; however as YTBassumed a constant Coriolis force, the model cannot explain the intensi caiton ofPSI near 28.9 as observed in the model of MW; in addition, the near-ineartial wavescan propagate a signi cant distance away from the latitude of 28.9 .This thesis extends the YTB model by allowing for a linearly varying Coriolisparameter ( -e ect) as well as eddy di usion. A linear stability analysis showsthat the near-inertial wave  eld is unstable to perturbations. We show that the -e ect results in a shortening in wave length as the near-inertial waves propagatesouth; horizontal eddy di usion is therefore enhanced to the south, and limits themeridional extent of PSI. The horizontal di usion also a ects the growth rate ofthe instability. A surprising result is that as the horizontal di usion vanishes, thesystem becomes stable; this can be demonstrated both analytically and numerically.Mathematically, the  -e ect renders the spatial di erential operator nonnor-mal, which is characterized with the aid of pseudo-spectra. Our results suggest thepossibility of large amplitude transient growth in near-inertial waves in regimes thatare asymptotically stable to perturbations.iiContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 IGWs and the dynamics of the ocean . . . . . . . . . . . . . . . . . . 31.2 Parametric subharmonic instability . . . . . . . . . . . . . . . . . . . 81.3 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Mathematic Description of PSI . . . . . . . . . . . . . . . . . . . . . 122.1 Derivation of the PSI model . . . . . . . . . . . . . . . . . . . . . . . 132.1.1 The pumpwave train . . . . . . . . . . . . . . . . . . . . . . . 172.1.2 The simpli ed amplitude equation . . . . . . . . . . . . . . . 182.2 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 203 Dissipative Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1 An initial value problem . . . . . . . . . . . . . . . . . . . . . . . . . 23iii3.1.1 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . 243.2 Normal modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.2.1 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . 283.2.2 Analytical results . . . . . . . . . . . . . . . . . . . . . . . . . 373.2.3 Energetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.3 More initial value problems . . . . . . . . . . . . . . . . . . . . . . . 533.3.1 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . 533.3.2 Boundary modes . . . . . . . . . . . . . . . . . . . . . . . . . 543.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584 PSI in the Inviscid Limit . . . . . . . . . . . . . . . . . . . . . . . . . 594.1 Stability of the inviscid problem . . . . . . . . . . . . . . . . . . . . 604.2 Non-normal operators and transient growth . . . . . . . . . . . . . . 624.2.1 A simple example . . . . . . . . . . . . . . . . . . . . . . . . 634.2.2 Non-normal operators and  uid dynamics . . . . . . . . . . . 664.3 Pseudospectra and PSI . . . . . . . . . . . . . . . . . . . . . . . . . . 664.3.1 Pseudospectra . . . . . . . . . . . . . . . . . . . . . . . . . . 664.3.2 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . 674.3.3 Numerical di culties . . . . . . . . . . . . . . . . . . . . . . . 734.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 775.1 Horizontal dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . 775.2 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.3 Outstanding issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84ivList of Figures3.1 An initial value problem . . . . . . . . . . . . . . . . . . . . . . . . . 273.2 E ect of dissipation on the spectrum: l = 0 . . . . . . . . . . . . . . 303.3 E ect of dissipation on the spectrum: l = 0:1 . . . . . . . . . . . . . 313.4 E ect of computational domain on the eigenvalues . . . . . . . . . . 333.5 E ect of the parameters on the growth rate. . . . . . . . . . . . . . . 363.6 Energetics of PSI: l = 0 . . . . . . . . . . . . . . . . . . . . . . . . . 483.7 Energetics of PSI: l = 0:1 . . . . . . . . . . . . . . . . . . . . . . . . 493.8 Energetics of PSI near the stability boundary . . . . . . . . . . . . . 523.9 An initial value problem with transient growth . . . . . . . . . . . . 553.10 An initial value problem with stable normal modes . . . . . . . . . . 563.11 An initial value problem with boundary modes . . . . . . . . . . . . 574.1 Transient growth for a simple dynamical system . . . . . . . . . . . . 654.2 Pseudospectra for a stable system . . . . . . . . . . . . . . . . . . . . 694.3 Transient growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.4 Initial condition and transient growth . . . . . . . . . . . . . . . . . 724.5 E ect of non-normality on calculation of eigenvalues . . . . . . . . . 755.1 Pseudospectra for a stable system . . . . . . . . . . . . . . . . . . . . 81vAcknowledgementsFirst and foremost, I would like to thank my supervisor Neil Balmforth for being aninspiring mentor over the last three years. It has been a pleasure working with himand I have learnt a great deal from him, which will no doubt be invaluable to mycareer. I would also like to thank Kraig Winters and Bill Young for their commentsfor this manuscript.Navigating myself through the degree would’ve been much more di cultwithout the help from the faculty members at UBC. I would like to thank DouwSteyn, Philip Loewen, Eric Cytrynbaum, Mark Maclean, George Bluman, andRichard Anstee for their constant guidance, which helped me grow immensely inmy research and teaching.I would also like to thank the sta of the Department of Mathematics andthe Institude of Applied Mathematics, especially Lee Yupitan for answering thecountless administrative questions I had and Marek Labecki for his help on solvingcomputer 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 encouragementsand love.The research is partially supported by the NSERC PGS-M scholarship.viChapter 1IntroductionThe 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 existdue to gravity: if the water surface is perturbed upwards, gravity will act on thewater 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 willnow act on the water surrounding the depression to  ll in the depression. Thissequence will occur repeatedly, which results in the familiar waves we observe onthe ocean surface.Gravity waves are not con ned to surface of the ocean, as they also existunderneath the surface. Due to di erences in temperature and salinity, the densityin the ocean typically increases with depth, and therefore when a parcel of water isperturbed downwards, the denser ambient  uid provides a positive buoyancy thatacts 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 di er from their coun-1terparts on the surface. Firstly, the sharp density contrast between air and waterresults in a much higher frequency of oscillation for the surface waves; the typicaltime 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 ais al afrequency N:N =s g 0d dz; (1.1)where g is the gravitational acceleration,  0 is the characteristic density, and z is thevertical coordinate. Below the ‘mixed layer’ on the surface of the ocean, which isuniform in density and extends up to 50m, the buoyancy frequency decreases from10 2s 1 to about 10 4s 1 in the deep ocean (about 4000m) (Thorpe, 2005); thesefrequencies correspond to time-scales ranging from minutes to hours.The IGWs also have a much larger spatial scale compared to surface gravitywaves. For example the internal waves associated with the diurnal tide can span theentire depth of the ocean (103m), while having horizontal scales measured in tens ofkilometres, which is in stark contrast to surface waves that typically have heights of0:1 10m and horizontal wavelengths of perhaps a tens of meters. A consequenceof the large horizontal scale of IGWs is that IGWs are susceptible to the Coriolise ect, which is re ected in the dispersion relation:k2 +l2m2 =!2 f2N2 !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 theCoriolis parameter, which depends on the earth’s rate of rotation ( = 7:27  1Throughout this thesis we will adhere to the convention in oceanography by orientingthe positive x axis towards east, positive y axis towards north, and positive z upwards210 5s 1) and latitude ( ):f = 2  sin : (1.3)Since the wavenumbers are real, and N  f in general, an implication of (1.2) isthat an IGWs can only exist if N  ! f. In addition one can use (1.2) to showthe direction of the group velocity makes an angle of  with the horizontal plane,where  satis es:tan =rk2 +l2m2 =r!2 f2N2 !2: (1.4)Since the energy propagates along the direction of the group velocity, this shows thatin general IGWs transports energy both in the horizontal and vertical direction,which is in contrast to the surface waves where energy is propagated along thesurface. Two limiting cases can be considered:1. !!N: in this case the angle of propagation tends to 90 , which correspondsto vertical oscillations with frequency N.2. !!f (The near-inertial limit): in this limit the angle of propagation tendsto 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 equationto 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 thehorizontal scale (k2 +l2) 1=2.1.1 IGWs and the dynamics of the oceanAs warm water has a lower density than cold water, the constant heating at thesurface will in principle result in a warm layer near the surface, and the rest of theocean will consist of cold dense water underneath. In reality the density pro le is3much smoother due to mixing, and one can  nd circulations that span the entiredepth of the ocean (the thermo-haline circulation for example ). Molecular di usionitself cannot account for the mixing, as the di usivity of heat and salt are in the orderof 10 6 10 7m2s 1, while the di usivity needed to maintain the strati cation isestimated to be 10 4m2s 1 (Munk & Wunsch, 1998). It is believed that the breakingof internal waves is responsible for the turbulence, and the diapycnal (across densitysurfaces) mixing that follows (Laurent & Garrett, 2002).Munk & Wunsch (1998) estimated that 2:1TW of energy is needed to main-tain the strati cation in the abyssal ocean via mixing; as an additional 2:6TW ofenergy is dissipated through turbulence mixing in the shallow oceans, a total energysource of 4.7TW needed. Wind induced mixing, either directly near the surface, orindirectly 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 oftides due to the sun and the moon; amongst the modes, the semi-diurnal componentof the lunar tide (denoted by M2) is the dominate source, accounting for 2.5TW.The e ect of the moon’s gravity  eld on the earth is not uniform everywhereon the surface, and the minute di erence is su cient to cause a deformation inthe ocean: the result is that the ocean bulges out along the earth-moon axis, bothtowards and away from the moon. As the earth rotates around its axis, a point onthe surface over the ocean will experience a rise and a subsequent fall in the sea-levelas it passes through a bulge. As most locations will pass through both bulges, theresult is the semi-diurnal tide that occurs approximately twice daily, or every 12hours and 25 minutes.Since the earth’s rotation rate (7:292 10 5s 1) is much higher than the rateat which the moon orbits around the earth (2:662 10 6s 1), we can view the earth4as a sphere moving through a viscous  lm of  uid held in place by the gravity of themoon. The friction between the earth and the ocean causes the earth’s rotationalrate 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 togain rotational momentum and increase its distance from the earth. The rotationalenergy is partially dissipated in the ocean, and partially transferred into the orbitalenergy of the moon, causing it to move further away from the earth. The modernvalue 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 oceanThe subject of tidal dissipation can be dated back to Sir G.I. Taylor (1920) and SirHarold Je reys (1921). Taylor estimated the the frictional dissipation over the Irishsea alone to be 0.051 TW, while Je reys estimated the global dissipation over shallowseas to be 1.1 TW; later Miller (1966) pegged the value at 1.6TW, while recentvalues inferred by satellite data is for the M2 mode is 1.8TW (Egbert, 1997; Munk& Wunsch, 1998). The discrepancy between the estimated dissipation over shallowseas and the total dissipation calculated from astronomical data can be attributed tothe 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 deepocean has been thought to be insigni cant until recently (Garrett, 2003). Part ofthe reason for the revision is due to the availability of satellite altimetry data, whichsuggests that dissipation over the deep ocean is responsible for converting 1TW ofenergy (Egbert & Ray, 2000, 2001).5The mechanism behind the tidal conversion is simple: as the tide movesaround the earth, the water in the ocean is moved back and forth; in the presence ofsubmarine ridges, the density surfaces are perturbed and internal gravity waves arelaunched in a fashion similar to lee waves generated by mountains. These internalgravity waves (called internal tides) then propagate from the bottom up into theopen ocean, carrying energy upwards that ultimately leads to overturning and mix-ing (Garrett, 2003). These internal tides can be observed with the aid of satellitealtimetry: for example, Ray & Mitchum (1997) reported internal tides propagatingaway up to 500km away from the Hawaiian Islands, primary as mode-1 baroclinictide. The baroclinic tides often have the same frequency as the underlying barotropictide due to the weak background tidal current (Garrett & Kunze, 2007), while thehorizontal structure is around 20-100km (Laurent & Garrett, 2002). The Richard-son number for these internal tides exceeds unity, suggesting that the internal tidecannot lead to overturning directly; it is believed that energy is cascaded down vianon-linear wave-wave interaction into smaller scale internal waves, which are moree ective at causing over-turnings (Laurent & Garrett, 2002)..Non-linear wave interactions in internal gravity waves are thought to beresponsible for shaping and maintaining the energy spectra of internal waves in theopen ocean into the so-called Garrett-Munk (GM) spectrum (M uller et al., 1986;Thorpe, 2005). Garrett & Munk (1972, 1975) proposed a model for the observedenergy spectra of internal waves that is in remarkable agreement with experimentaldata from di erent sites; this suggests that there is a single universal mechanismresponsible for cascading energy from large scale motion to small scale mixing.As abyssal mixing is an essential component to the dynamics of large scale oceancirculations that in turns play a pivotal role in our climate, the GM spectrum and6ocean 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-numbersatisfy 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 fromthe  rst two modes into the third. Bretherton (1964) considered the interactionbetween a triad of discrete modes, where !1 + !2 = !3 and ~k1 + ~k2 = ~k3, andprovided a heuristic argument for the evolution equation of the amplitudes ai:i!1@a1@t = ^Da 2a3; (1.7)i!2@a2@t = ^Da 1a3; (1.8)i!3@a3@t = ^Da1a2; (1.9)where ^D is a real coe cient that describes the interaction, and generally depends onthe wave-numbers. Hasselmann (1966, 1967) studied weak non-linear wave interac-tion for a continuous  eld of waves that is assumed to be Gaussian and statisticallystationary (‘the random phase hypothesis’). McComas & Bretherton (1977) used asimilar method to investigate the resonance interaction in the GM spectrum; theynumerically computed the energy transfer in the spectrum, and concluded that theprocess is dominated by three classes of interactions: induced di usion, elastic scat-tering, and parametric subharmonic instability. For Induced di usion, wave action(wave energy divided by frequency) di uses from a high-frequency, high-wavenumber7mode to one that has similar frequency and wavenumber. Elastic scattering occursbetween 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 twohigh frequency modes leaving the third low frequency mode largely unchanged. Thelast class of interaction is parametric subharmonic instability (PSI), where a largescale wave decays into two small scale waves with half the frequency. An importantobservation is that among the three pathways PSI is the only mechanism whereenergy is transferred into waves with a much di erent wave number. For a thoroughdiscussion, please see McComas & Bretherton (1977) or M uller et al. (1986).1.2 Parametric subharmonic instabilityParametric subharmonic instability allows a wave to spontaneously excite daughterwaves 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 angularfrequency M2 = 2 =12:42hr = 1:405  10 4s 1, the daughter waves will havefrequency ! = M2=2 = 7:026 10 5s 1. As mentioned earlier, internal gravity wavescan exist only if the frequency is greater than the local Coriolis frequency; thus forthe 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 tothe critical latitude 28:9 . At the critical latitude the daughter waves have frequencythat is equal to the local Coriolis frequency, and therefore they are near-inertialoscillations (NIO).While PSI have received some attention over the 1970’s, the interest wanedas it was believed that the energy transfer via PSI occurs over a time-scale of 100days for the mode-1 internal tide, which is too slow compared to other processes to8remove energy e ciently (Olbers & Pomphrey, 1981); however the interest in PSIhave 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 internalwave  eld was perturbed by injecting energy into the  rst vertical mode at the M2freqency. For the experiment at 49 N, the spectra did not di er signi cantly over10 inertial period; on the other hand, the energy was quickly transfered to the highwave-number modes in the f <!< 2f range within 10 inertial periods at 28 N.MacKinnon & Winters (2005) (henceforth MW) solved the full 3-D equationof motions numerically with a linearly strati ed background; a mode-1 internal tidewas excited at the southern boundary, which propagated northward in the presenceof background noise. After the initial spin-up, instabilities started to develop nearthe critical latitude. The instability had much  ner vertical structure compared tothe mode-1 internal tide, and a plot of temporal spectra plot reveals a high spectralcontent over the M2=2 frequency. At the same time the energy  ux associatedwith internal tide decreased by 62.5% as it crossed the critical latitude. In additionthe time scale of the energy transfer was estimated to be 10 days, which was anorder of magnitude smaller than the estimation of Olbers & Pomphrey (1981). MWattributed the disparity to the fact that the estimation of Olbers & Pomphrey wasbased on the random phase assumption, but from the numerical simulation theinteraction were clearly much more coherent. The e ciency of PSI led the authorsto describe this as the ‘subtropical catastrophe’. Simmons (2008) ran a global oceanmodel with realistic geometry, and observed energy transfer into the M2=2 modesnear the critical latitude as well. A longitudinal cross section near a internal tidegeneration site showed that PSI intensi ed near the surface, possibly due to thenon-linear strati cation.9Recent  eld measurements also pointed to the presence of PSI activity inthe ocean. For example, Rainville & Pinkel (2006) compared the vertical pro lesat a site near Hawaii with strong internal tide activity (near eld) to a site that isfurther away, and found a strong link between the semi-diurnal energy  ux at thenear eld site and the diurnal energy  ux at the far eld site; this suggested that asthe internal tide radiates away from the generation site it transfers energy down tothe diurnal mode via PSI. Carter & Gregg (2006) came to a similar conclusion basedon another set of data from the Hawaii area. Through a bicoherence analysis theywere able to demonstrate that the strength of the M2=2 diurnal mode was correlatedto the strength of the M2 internal tide. Alford et al. (2007) estimated the reductionin energy  ux as the internal tide crosses the critical latitude to be 10-20%, whichis signi cantly smaller than the 62% in MW; they suggested that the presence ofhigher modes de-tuned the resonance and leads to a reduction in e ciency.Young et al. (2008) (henceforth YTB) constructed an analytical model basedon a multiple time scale approach, where the time scale over which evolution of thePSI instability takes place was slower than the inertial time scale f 1. With thisassumption the Boussinesq equations were simpli ed into a PDE that describesthe evolution of the PSI wave. For simplicity the latitude was restricted to thecritical latitude 28:9 , so that the daughter waves are near-inertial. The stability ofa plane wave internal tide in a linearly strati ed ocean was considered, and usingrealistic parameters the growth rate of the unstable near-inertial mode agreed withMW in order of magnitude. YTB also considered the stability problem in a morerealistic density strati cation; using Gill’s model for the ocean strati cation, themost unstable normal modes for the near-inertial  eld were found to be concentratedjust under the mixed layer.101.3 GoalsThe goal of this thesis is to expand on the model of PSI derived by YTB; in particularwe will allow the Coriolis force to vary with latitude, and assess its e ect on PSI.The rational behind this modi cation is that over the time-scale for PSI (10 days),there is signi cant propagation of the near-inertial waves, and thus the constantCoriolis parameter assumption may not be justi ed. In addition, we would like toinvestigate the e ect of dissipation to determine if it is crucial to the instability,in particular to the understanding to vertical scale selection that appeared in thenumerical simulations.11Chapter 2Mathematic Description of PSIA familiar example of resonance occurs when a simple harmonic oscillator is forcedat its natural frequency. Here we will  rst consider the following equation, whichdescribes an oscillator forced by a weak (  1) periodic driver with frequency 2:y00 +!2y =  y e2it: (2.1)The presence of y (the complex conjugate of y) in the forcing term indicates thatthe forcing is moderated by feedback from the oscillator. Since the natural frequencyfor the oscillation is  !, the homogenous solutions are ei!t and e i!t. These twosolutions will result in forcing terms that are proportional to ei(2 !)t and ei(2+!)trespectively on the RHS of (2.1). As long as !6= 1, the forcing will be at a frequencythat di ers from the natural frequency, and hence no resonance will occur.On the other hand, when ! = 1 both the homogenous solution and thecorresponding forcing will contain eit, and in this case resonance will occur. Onecan analyse (2.1) through a multiple time scale approach. As the coupling is weakdue to the small parameter  , we would expect the amplitude of y to grow over a slowtime scale T =  t. If we seek an ansatz to (2.1) of the form y(t;T) = A(T)ei(t  =4),12we will have:2i A0(T)ei(t  =4) +O( 2) =  A(T)ei(t+ =4); (2.2)Equating the O( ) terms leads to a di erential equation for A(T):dAdT =12A (2.3)The above equations can be considered as a toy model for PSI, as the naturalfrequency of the oscillator is half the frequency of the driver, and therefore theresonant condition !1 +!2 = !3 (where !1 = !2 = !=2) is satis ed. In addition themodel predicts that the amplitude of y(t) will grow as e t, indicating that the driveris transferring energy into the oscillator, and the amplitude grows exponentially overthe 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 modi cations. The  rst is the inclusion of the latitudinaldependence 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 notconsidered by Young et al. As these two modi cations can be easily included in therigorous 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 approachby assuming a priori that the PSI modes will have frequency M2=2, and derive thegoverning equations.2.1 Derivation of the PSI modelThe starting point is the hydrostatic Boussinesq equations linearized around thebackground  ow  elds (U;V;W;B;P), which denote the velocities, buoyancy, and13pressure respectively. Since we are interested in the PSI in a mode-1 internal tidewith frequency 2f0 = M2, the background  elds will be proportional to e 2if0t andhave 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 ofthe background  elds are time independent, and do not diminish as energy is beingextracted from the internal tide. It can be viewed as a linearized model for the triadresonance, where the amplitudes of the daughter waves are small, and therefore thefeedback 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 +Vuy +Wuz +uUx +vUy +wUz fv +px + ( r2)nu = 0; (2.4)vt +Uvx +Vvy +Wvz +uVx +vVy +wVz +fu+py + ( r2)nv = 0; (2.5) b+pz = 0; (2.6)ux +vy +wz = 0; (2.7)bt +Ubx +Vby +Wbz +uBx +vBy +wBz +wN20 = 0; (2.8)where r2 = @2x +@2y +@2z. The last terms in (2.4) and (2.5) describe the dissipativeprocess, and the exponent n = 0;1;2 corresponds to the linear drag, eddy di usion,and hyperdissipation model respectively. In the ocean, the presence of strati cationimplies that the turbulence is not necessarily isotropic, and horizontal mixing occursmuch more readily than mixing vertically due to stable vertical strati cation. Forexample, the measured eddy horizontal di usivity in a dye-release experiment inopen ocean ranged from 10 2 to 103m2s 1 depending on the scale (Ledwell et al.,1998), while the vertical eddy di usivity is of the order of 10 5m2s 1 (Ledwell et al.,1993, 1998); therefore in the eddy di usion model, the dissipation term is split into14its vertical and horizontal components: ( r2)nu =  H( r2H)nu+ V ( @zz)nu; (2.9)wherer2H = @2x+@2y, and  H and  V are the horizontal and vertical eddy di usivityrespectively. We will ignore the di usion of buoyancy b.We can further simplify the equations using the assumptions stated in YTB.Since the background  ow has a much larger vertical scale, any z-derivatives of thebackground variables can be ignored in (2.4) and (2.5). YTB assumed the buoyancyperturbations 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 thecomplex variablesU = u+iv;  = x+iy; (2.10)the horizontal derivatives can then be written as:@ = 12(@x i@y); @  = 12(@x +i@y): (2.11)We can then simplify (2.4) to (2.8) to:Ut +UUx +VUy +WUz + (Ux +iVx)u+ (Uy +iVy)v+ifU + 2p  + H( r2H)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 speci cally, we will assume that the15solutions are of the form:U = LA(x;y;z;t)e if0t; (2.16)where L = @zf20N 2@z. The other near-inertial  elds 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 theassumption that the near-inertial period f 10 is much larger than the time scale onwhich b evolves, and the variation of the amplitude can be ignore. b can then beobtained by integrating wN2 with respect to time. We will also linearize the Coriolisparameter 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  rst 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 thathave 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) canthen be discarded. Among the coupling terms, only (Ux + iVx)u + (Uy + iVy)v areresonant:(Ux +iVx)u+ (Uy +iVy)v = 12U(Ux +Vy +iVx iUy) + 12U (Ux Vy +iVx +iUy)= (U +iV)  LA e if0t +NRT; (2.22)16We will denote (U +iV)  by 12 . Ignoring the NRT, the amplitude equation (2.12)then becomes:LAt + i yLA+ 12if0r2HA+ H( r2H)nLA+ V ( @zz)nLA+ 12 LA = 0: (2.23)2.1.1 The pumpwave trainAs in Young et al, we will consider a pumpwave that is periodic in the horizontaldirections with horizontal wave numbers (k;l), and frequency ! = 2f0 + 2 , where2  2f0 is de ned as the de-tuning frequency. Assuming that the pressure  eldassociated with the pumpwave is given by:P = acos(kx+ly !t)p1(z) (2.24)where p1(z) is the eigenfunction of the  rst vertical mode. We can then retrieve thecorresponding horizontal velocities through the linear Boussinesq equations. De n-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 de nition of  , and retaining only the parts that resultsin resonance: = 2(U +iV)  = ia(k +il)22(! f0) ei(kx+ly 2 t)p1(z): (2.27)We can ignore the argument of (k + il)2 since it just represents a constant phaseshift of a. De ning = a(k2 +l2)2(! f0)  a(k2 +l2)2f0 (2.28)17as 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 tobe linearly strati ed and have a constant buoyancy frequency: N(z) = N0. In thiscase the  rst vertical mode will be:p1(z) = cos  Hz ; (2.30)where H is the depth of the ocean. The dispersion relationship for internal gravitywaves then results in a constraint on the horizontal wavenumbers k and l:k2 +l2 = !2 f20N20  !2 2H2  3 2f20N20H2 ; (2.31)where the buoyancy frequency N0 is assumed to be much larger than f0.2.1.2 The simpli ed amplitude equationWe now look for solutions of the formA cos(mz)e i theik1xA1(y;t) + eik2xA 2(y;t)i; (2.32)where m is the local vertical wave number, and the horizontal wave numbers satisfythe resonant condition:k1 +k2 = k: (2.33)If we assume that the vertical scale of the instabilities, m 1, is small compared toscale over which the buoyancy frequency and the vertical structure of the pump wavevary, 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 f0mN0; (2.34)18and the di erential operator L becomes:LA= ~m2A; (2.35)Substituting (2.29) and (2.32) into the amplitude equation (2.23) and separating thetwo harmonic components, we obtain a pair of coupled partial di erential equationsfor the amplitudes Aj(y;t). With the horizontal length scale ( ) and time scale ( )de ned as = f02 ~m2  1=3;  = 1  ; (2.36)the non-dimensionalized equations are:iA1~t = A1~y~y + ~y ~b ~a ~  A1 + ~ ei~l~yA 2 i~ H  1A1 + ( 1)nA1~y2n ;(2.37) iA2~t = A2~y~y + ~y ~b+ ~a ~  A2 + ~ e i~l~yA 1 +i~ H  2A2 + ( 1)nA2~y2n :(2.38)In the above equations, the ~t =   t and ~y = y= are the non-dimensionalize in-dependent variables, and ~kj =  kj and ~l =  l are the non-dimensionalized wavenumbers. The non-dimensionalized pump wave amplitude ~ , de-tuning frequency ~ and dissipation parameter ~ are given by:~ =  2  ; ~ =    ; ~ =  H 2n+1 : (2.39)The parameters ~a and ~b are de ned as:~a =  2(k21 k22)2 ;~b =  2(k21 +k22)2 : (2.40)The e ective wave number is de ned by: i = ~k2ni 1 +  Vm2n Hk2ni (2.41)19(2.37) and (2.38) can be further simpli ed by removing some parameters, as theire ects on the system can be understood easily. For example ~b+ ~ can be removedvia a shift in spatial coordinates ~y ~b ~ ! ~y together with a suitable constantphase shift in A1;2, and therefore the e ect of these two parameters is limit to ashift of the origin. If we assume that k1 = k2, the dissipative term  i~ U i can alsobe removed by a change of variable: Ai!e ~ H itAi ; this suggests that the e ectof dissipation in the zonal and vertical direction is to suppress the amplitude of thePSI waves without altering the meridional structure, and hence we can ignore it.Finally we can also remove the parameter ~a via (A1, A2)! ei~atA1, e i~atA2 , andtherefore the e ect of ~a is limited to a frequency shift.With the above simpli cations, 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 estimationWe will close this chapter by providing a summary of the oceanographically relevantvalues for the physical parameters used in our scaling. In table 2.1, the pumpwavenumber l and amplitude  are similar to the ones used in YTB. The buoyancyfrequency of 0.0087s 1 is a value that is typically found in regions just underneaththe mixed layer. These values, together with a vertical wave number of 100m, givea time scale of 30.5 days and a length scale of 19km; this length scale correspondsto a change of 0.171 in latitude.20M2 tidal frequency 1:405 10 4s 1Coriolis frequency, f0 = 12M2 7:02 10 5s 1Planetary vorticity gradient,  2:003 10 11 m 1 s 1Buoyancy frequency, N0 2 =720s = 0:0087s 1Vertical wavenumber, m 2 =100mPump wavenumber, l 2 =125kmPump amplitude  1=5daysLength scale,  = (f0N20=2f20m2 )1=3 19kmTime scale, 1=  30.5 daysEddy viscosity scale,  3 1:36 102m2s 1Non-dimensional pump amplitude, ^ =  =(2  ) 3:04nondimesnional pump wavenumber, ^l =  l 0:946Table 2.1: Parameter values21Chapter 3Dissipative ProblemIn this chapter we will explore the solutions to the PDEs (2.42) and (2.43) derivedin the previous chapter. We will mostly focus on the regular eddy dissipation modelby taking n = 1, which results iniA1t = (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 densityE can be de ned as:E 12(u2 +v2) = 12jLAj2: (3.3)E is proportional tojA1j2+jA2j2, and thusEn =jAnj2 forn = 1; 2 can be interpretedas energy density for our equations. From (3.1) and (3.2), we can derive the followingenergy equations:E1t +J1y = D1 +S; (3.4)E2t +J2y = D2 +S;: (3.5)22with the meridional energy  ux de ned asJn(y;t) i( 1)n A nAny AnA ny   H AnyA n +A nyAn : (3.6)The energy loss due to dissipation is given byDn(y;t) 2 HjA2nyj; (3.7)and  nally the energy source due to PSI is:S = i  eilyA2A 1 e ilyA 2A1 : (3.8)3.1 An initial value problemThe controlling parameters in (3.1) and (3.2) are the dissipation ( H), the amplitudeof the pumpwave/forcing ( ), and the meridional wave number of the pumpwave (l).In the absence of the pumpwave and e ect of dissipation, the equations reduce tothe Schr odinger equation with a linear potential, which describes a particle fallingunder the in uence of gravity. This is can be explained via a slowly varying waveapproximation. Assuming that  H and  are small (and thus can be ignored), theA1 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. De ning the frequency ^! =  t and wave number ^k =  y, thedispersion relation for (3.9) is^! = ^k2 +y: (3.10)Di erentiating (3.10) with respect to y, and noting that ^!y =  ty = ^kt, we nowhave a PDE for ^k:^kt + 2^k^ky = 1; (3.11)23which can be solved using the method of characteristics. The solution is given by^k = t + ^k0 on the characteristic curves y y0 = t2 + 2^k0t, where ^k0 and y0 arerespectively the initial wave number and position of the wave packet. The form ofthe characteristic curves suggest that the wave packet will indeed follow a ballistictrajectory, with the wave-packet having velocity dy=dt = 2k = 2(t ^k0); this suggestthat for t> ^k0 the wave-packet will have negative velocity, and will be moving south.Eliminating t from the solution results in:^k2 = ^k20 + (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 problemfor PSI. Consider a perturbation in the near-inertial wave  eld in the form of aGaussian wave packet with wave-number k0 > 0 and :A1(y;0) = e (y y0)2= cos(k0y) = 12 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 wavepackets. The component with e ik0y travels south, while the component with eik0ywill  rst travel north up to y = k20 +y0 before turning south. In the presence of thepump-wave (i.e.  > 0), the near-inertial wave-packets will grow due to PSI, andwe would like to understand how the growth depends on the parameters.3.1.1 Numerical schemeAs with other wave-type problems on an in nite domain, We will demand that theamplitude of the near-intertial waves decay to 0 as y ! 1; however there aresome practical issues with truncating the problem to a  nite domain [yL;yR] for thenumerical calculations as the boundary conditions (BC) at in nity are now applied24on a  nite interval. For the right boundary yR, the wave packet will be re ectedbefore it reaches yR provided that the initial kinetic energy is not too high, and thuswe 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 thewave number increases as pjyj as y! 1 (see (3.12)). This presents us with adilemma: on one hand the domain needs to be large in order for the wave amplitudeto decay su ciently to justify imposing a zero Dirichlet condition, but as the domainis extended by moving yR towards  1 the grid needs to be re ned to resolve theincreasingly 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 ofABCs: PDE-based and material based (for a brief review, see Zheng (2007)). PDE-based ABC relies on the ability to factor the di erential operator, which in turn canbe used to derive a radiative boundary condition that only allows out-going waves atthe 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 basedABC is the perfectly matched layer (PML), which was used by Berenger (1994) forthe study of Maxwell’s equations in electro-magnetics. It was also used sucessfully inboth linear and non-linear Schr odinger’s equation with variable potentials (Zheng,2007).As it is not clear how our equations can be factored to derive a radiative BCon the boundary, we will adopt the PML technique as outlined in Zheng (2007) forthis 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 PML25is identical to the original PDE in the original domain of interest, while the PDEis modi ed in the sponge layer [yLR;yL] by introducing an additional dissipativeterm that strongly damps out-going waves. The solutions are expected to decayexponentially, 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 re ected back into the boundary, and thus allowing us to considercases where the dissipation is too small to damp out the perturbations before itreaches the boundary. With the introduction of a PML, the second derivatives aremodi ed via:A1yy! 11 i (y) ddy 11 i (y)A1y (3.14)A2yy! 11 +i (y) ddy 11 +i (y)A2y : (3.15)The absorption function  (y) is chosen to be: (y) =8><>: 0(y yl)4 if yLR <y<yL0 if yL y<yR(3.16)Equations (3.1) and (3.2) together with the initial conditionA1(y;0) = e (y+25)2=5 cos(p26y); A2(y;0) = 0 (3.17)are solved using a Crank-Nicolson  nite di erence scheme, which has the advantageof being unconditionally stable (Garcia, 1999).We will now present a sample solution, where the computational grid isuniform 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  eld, jA1j and jA2j, are shown in Figure3.1. We can see that the wave packet splits up into two components that travel in26opposite directions. The northward propagating packet follows a ballistic trajectorywhile spreading out. At the same time we can see that as the wave packet movesthrough 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 theamplitude rapidly drops o to the south. This is not surprising as the e ect ofdissipation will be more potent as the wave-number increases to the south.!"#$%&’()’!*+ !,+ !-+ !.+ !/+ !0+ !1+ !2+ !)+ + )++20.,)+)2)0).),2+"#3%&’(2’!*+ !,+ !-+ !.+ !/+ !0+ !1+ !2+ !)+ + )++20.,)+)2)0).),2+Figure 3.1: The results from a numerical computation of the PSI model. jA1j andjA2jare 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 layerlocated at [ 100; 85]. Only the data from [ 90;10] is shown here.3.2 Normal modesThe 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 theorigin, and therefore we do not need to include a sponge layer in our problem.The domain can be truncated at a point su ciently far south, and a zero Dirichletboundary 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 thesame as the PDE: jAnj!0 as y! Numerical solutionTo investigate the e ect 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  nal di erence scheme, and the eigenvalueproblem becomes a matrix eigenvalue problem; the matrix is large but sparse, whichcan be handled via the eigs function in MATLAB. eigs is an iterative numericalscheme that estimates the eigenvalues for matrices, and it is useful when only partof the spectra is of interest.We will compute the spectra over [yL; 50] for some yL, and impose zeroDirichlet boundary condition at both end-points. We will then take the limit asyL! 1; if the eigenvalues converge to a limit, the limit can be interpreted as theeigenvalues of the problem over the in nite interval. We found that the eigenvaluesand eigenfunctions only di er slightly with the zero Dirichlet BC is replaced by asponge layer at the southern boundary, and therefore we will ignore the sponge layerin this section.Zonally propagating pump-wave: l = 0Figure 3.2 shows the result for a series of calculations with  = 1, and l = 0. Theeigenvalue problem is solved over a domain of [ 100; 50], and the real part of the28normalized eigenfunctions <(A1) over [ 60;10] are plotted on the left column of gure 3.2.  = 1, 10 1, 10 2, 10 4, and 0 were used in the calculations; the plotsare arranged in decreasing order of  from top to bottom. One can see that inall cases, the eigenfuction is oscillatory to the south (when y < 0), while to thenorth of the turning point y = 0, the eigenfunction quickly decays to 0, indicatingthat the near-inertial waves are evanescent as expected. On the bottom row, thenumerically computed eigenfunction (blue curve) is almost indistinguishable fromthe Airy function (red curve) in the absence of dissipation; inx3.2.2 we will show thatthe Airy function is the exact solution to the inviscid problem. When dissipationis increased to 10 4, the eigenfunction still resembles the airy function, but theamplitude to the south is signi cantly 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 decayingmodes), or purely imaginary (neutral modes). The fastest growing mode has agrowth rate close to 1. As the dissipation increases, to  = 10 2, the majority ofthe modes become decaying, with a handful of modes that remain unstable. Whenthe dissipation is further increased to 1, all but one of the modes become stable, andthe growth rate of the unstable mode is now signi cantly lower then the inviscidcase. The eigenvalues becomes insensitive to yL for yL < 100, and therefore thecomputed spectra should be a good approximation to the spectra on the in niteinterval.29−60 −50 −40 −30 −20 −10 0 10−101y−101−101Re A1−101−101−5 0 5−101!I−101−101!r−101−101µ=0µ=10−4µ=10−2µ=10−1µ=1Figure 3.2: E ect of dissipation on the spectrum: l = 0: The dissipation parameterused (from top to bottom) are  = 1;10 1;10 2;10 4 , and 0. On the left columnthe real part of the eigenfuctions associated with the fastest growing modes areplotted in blue, while in the last two panel the airy function Ai(y) is also plotted inred. On the right column the spectra are plotted, with the abiscissa and ordinatebeing respectively the imaginary part and the real part (the growth rate) of theeigenvalues.  = 0 for all calculations.More general case: l6= 0We will now turn to the case where l is no longer zero. In  gure 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 e ect of truncation. The results withyL = 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 scalefor the frequency of the spectra plot is di erent from  gure3.2.For large values of dissipation (top two rows in  gure 3.3), the results fromboth values of yL agree, showing that neither the eigenfunction nor the physically30−80 −60 −40 −20 0−101y−101−101−101Re A1−101−101  −2 −1 0 1 2−101!I−101−101−101!r−101−101  µ=10−4µ=0µ=10−3µ=10−2µ=10−1µ=1Figure 3.3: E ect of dissipation on the spectrum: l = 0:1: The dissipation parameterused (from top to bottom) are  = 1; 10 1; 10 2; 10 3; 10 4, and 0. On the leftcolumn the real part of the eigenfuctions associated with the fastest growing modesare plotted; in blue are the results from using a domain of [ 100;20], while thecurves in red are based on a domain of [ 150;20]. The corresponding spectra arerespectively plotted with blue dots and red circes. The abiscissa and ordinate arerespectively the imaginary part and the real part (the growth rate) of the eigenvalues. = 1 for all calculations31interesting portion of the spectra is a ected by the change in yL. As in the l = 0case, the eigenfunction decays quickly to zero due to dissipation, and there are asmall number of discrete unstable modes. When dissipation is decreased to 10 2, thetwo computations agree for the fastest growing mode, but the spectra plot revealsthat the domain length starts to a ect the computation. The computation withyL =  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 unstablemodes have zero frequency.As the dissipation is decreased further (bottom three rows), the picture isvery di erent from the l = 0 case. The spectra plots show that most unstable modesnow have non-zero frequency, and the frequency increases as the yL!1. The plotsof the eigenfunction reveal that the waves are no longer localized to the turninglatitude; instead, the bulk of the wave is located at the center of the domain, andmoves as the boundary moves southwards. We will refer these modes as ‘boundarymodes’.To get a better understanding of the transition between the boundary andlocalized 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 resultsare presented in  gure 3.4, where the frequency ( I) and growth rate ( R) of the mostunstable 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 isreduced as dissipation vanishes; the graph suggests that as we approach the inviscidlimit, 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 remains3210−410−310−210−11000.µ!R  PredictedyL=−80yL=−100yL=−120yL=−140yL=−160yL=−180yL=−20010−410−310−210−110000.!I−200 −180 −160 −140 −120 −100 −800.511.5yL!I,R  !I!R0.1 yL0.5Figure 3.4: E ect of computational domain on the eigenvalues. The frequency ( I,top panel) and growth rate ( R, bottom panel) of the most unstable modes areplotted 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 markedby dots. The predicted growth rate from (3.74) is also plotted in the bottom panelas the green curve. Inset:  R and  I as a function of yl in the inviscid limit.  = 1and l = 0:1 for all calculations.33dependent on yL. The e ect of yL on the eigenvalues for the inviscid problem ( = 0)is shown in the inset of  gure 3.4. The growth rate depends very weakly on yL, butthe frequency increases approximately as lpyL (plotted as the green dashed line).The eigenfunctions for the unstable modes in the inviscid problem resemble theboundary modes plotted in the bottom panel of  gure 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 oscillationsin the eigenfunctions are well resolved, and a change in grid size does not a ect theresult. It is also not due to the zero Dirichlet boundary condition imposed at thesouthern boundary, as the eigenvalues and eigenfunction of the boundary modes areonly a ected minimally if the zero Dirichlet BC is replaced by a sponge layer ora radiative (Sommerfeld) type boundary condition that eliminates incoming wavesat the boundary. Since the frequency of the boundary mode does not appear toconverge, and the location of peak of the eigenfunctions changes as the yL is movedto the south, the conclusion is that these modes are un-physical.We will now return to  gure 3.4. From the inviscid limit, the growth rate ofthese boundary modes rapidly decreases as dissipation is increased, with the modefrom yL = 200 being most susceptible to dissipation. There is a sudden break inthe graph when  inceases to 10 2:5, and the numerically computed growth ratesnow converge to a single curve, and no longer depends on yL. The frequency alsoundergoes 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 fromthe boundary modes to modes where both the eigenvalues and eigenfunctions areindependent of the location of the boundary. These new modes are physically mean-ingful since the choice of the locaion of the boundary does not a ect the structure34or growth rate of the mode. The eigenfunctions are similar to the localized modesin the top three panels of  gure 3.3.A curious feature is that the e ect of dissipation on these localized modes isnot monotonic: while an increase in dissipation for larger values of  stabilizes thesystem, the growth rate in fact increases with  when dissipation is small (between10 2:5 to 10 1:7); this result rather surprising, as one would normally expect dis-sipation to have a stabilizing e ect. We will explore this behaviour via the energyequations in x3.2.3.We will now present a couple of plots in  gure 3.5, which summarize the e ectof the various parameters on the growth rate for = 4 (top panel) and = 6 (bottompanel). The two plots are qualitatively similar. For l = 0 (top series),  R  whenthe dissipation is small, and decreases monotonically as  increases. For l6= 0, thenon-monotonic e ect of dissipation on growth rate observed earlier holds for a widerange of l. An increase in l reduces the growth rate, and for su ciently large valueof l, the normal modes become stable for small  . For example, when l = 2:5 and = 4 (bottom series in the top panel of  gure 3.5), the growth rate is negative for < 10 0:5.3510−210−1100101−1−0.500.511.522.533.54µγR(a)  = 410−210−1100101−10123456µγRA(b)  = 6Figure 3.5: E ect of the parameters on the growth rate. The growth rates areplotted as a function of  for di erent values of pump wavenumber l (dots); fromtop to bottom, l = 0, 0.5, 1,1.5, 2, and 2.5. The dashed lines are the growth wavespredicted by the short wave theory in x3.2.2(c.f. (3.74)). For all calculations, theeigenvalue problem was solve over the domain [ 100; 20].363.2.2 Analytical resultsWe will now develop some analytical results which provides support to our numeri-cally computed growth rates. The most important result is an analytical predictionfor limit where 0 <l 1 and   1, as it con rms our numerical results showingthat  !0 stabilizes the system for a non-zero l. The problem with  = l = 0 canbe solvable exactly, which is then used as the basis of a perturbative solution withl and  as a small parameter.A zonally propagating pump-wave in the inviscid limitWe will  rst 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 diagonalizedand solved exactly:264A1(y)A2(y)375 =264   i  p 2  2i +p 2  2  375264c1Ai(y p 2  2)c2Ai(y +p 2  2)375 (3.20)The fundamental solutions are Ai(y p 2  2), where Ai(y) is the Airy function ofthe  rst kind. The argument of Ai(y) must be real for Ai(y) to remain bounded asy! 1, which implies that  must be real and  2  2R 0. The latter conditionimplies that there is a continuous spectrum of normal modes, with growth ratessatisfying the inequality:    R  : (3.21)The maximum growth rate is  Rmax =  , which agrees with the numerically ob-tained growth rates (c.f.  gure 3.5 and the bottom panel of  gure 3.3). Theeigenfunctions associated with the fastest growing mode are A1(y) = Ai(y) andA2(y) = iAi(y). Using (3.6), It can be shown that the meridional energy  ux is37identically 0. The interpretation is that the fastest growing near-inertial oscillationsexcited 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 asthe present study): Rmax =p 2 (b+ )2; (3.22)Curiously the factor b+ in (3.22) is identical to the spatial shift we initially ignoredin the IVP problem, yet it plays no role in the growth rate for our analysis. It isperhaps instructive to re-concile the two results.Discarding the  -e ect 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 areperiodic in y: A1 = eil1y and A2 = e il2y. To satisfy the resonant condition, thewave numbers must satisfy l1 +l2 = 0, and the eigenvalues are:i = iq 2  l21 (b+ ) 2: (3.25)YTB considered the reduced 2-dimensional model in the absence of the   e ect,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 noreason to restrict ourselves to meridionally uniform solutions, and thus (3.25) is amore general result for the model of YTB. This equation suggests that maximumgrowth rate  max =  can be achieved when l1 = pb+ , and therefore there isscale selection in the y-direction at the critical latitude.38This scale selection behaviour is consistent with our model. If we retainb+ 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 Airyfunction is periodic with local wave-number pjy (b+ )j, while to the north thewave becomes evanescent. Now if (b +  ) 1, then the critical latitude y = 0 liesto the south of the turning point; therefore the waves at the critical latitude will beperiodic with wave-number p(b+ ), which is consistent with the analysis in thepreceding paragraph.Short wave theoryWe can use the analytical solution to the  = l = 0 problem in the previous sectionto make progress analytically for cases where the the parameters  and l are non-zerobut 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:f00 lg0 l24 f + g00 + lf0  l24 g yf ( + )g = 0; (3.27)g00 +lf0 l24 g  f00 + lg0 + l24 f yg (   )f = 0: (3.28)In the previous section, the numerical results suggest that the unstable modes havefrequency 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 parameterde ned as  =p , as well as the assumption that l scales as  3=2 (i.e. l = l3=2 3=2),39(3.27) and (3.28) becomesf00 l3=2 3=2g0 + 2g00 yf ( + )g = O( 5=2); (3.29)g00 +l3=2 3=2f0  2f00 yg (   )f = O( 5=2): (3.30)Inner solutionWe 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, wehave:f001=2 yf1=2 = 2 g1=2; (3.34)g001=2 yg1=2 = 0: (3.35)(3.35) is satis ed by g1=2 = rAi(y), where r is a constant that is to be determinedlater. Since (Ai0)00 = (Ai00)0 = (yAi)0, a particular solution to (3.34) isf1=2 = 2 rAi0(y): (3.36)Moving on to O( ) terms in (3.30), we have:g001  yg1 =  1Ai; (3.37)which has the particular solution:g1 =  1Ai0(y); (3.38)40To summarize, the inner solution is:f = Ai(y) + 1=22 rAi0(y) +   ; (3.39)g =  1=2rAi(y)   1Ai0(y) +   ; (3.40)The Outer RegionAs y! 1, standard asymptotic formulae for the Airy function give (Abramowitz& Stegun, 1964):Ai(y) 1p 4p y sin  +  4 and (3.41)Ai0(y) 4p yp cos  +  4 ; (3.42)where  (y) = 23jyj3=2 is the phase. Since the solution becomes increasingly oscil-latory, we anticipate that some of the derivatives we neglected in the inner regionmay become important. In addition, (3.42) suggests that the O( ) term in g(y)grows as y! 1, which leads to a break in the order. Both of these observationspoint to the existence of an outer solution which we will need to match onto theinner solution (which we will relabel as fin and gin) given by (3.39) and (3.40). Toproceed we will de ne Y =  y as the outer variable; notice that in terms of the outervariable,f  1=4 14p Y sin  (Y) 3=2 + 4 + 2 r 4p Y cos  (Y) 3=2 + 4  and (3.43)g  3=4 r4p Y sin  (Y) 3=2 + 4   1 4p Y cos  (Y) 3=2 + 4  ; (3.44)41We have scaled the inner solutions to remove the constant factor of 1=p . Motivatedby 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: 2d2foutdY2  l3=2 5=2dgoutdY + 4d2goutdY2  Y f (2 +  1)g = O( 5=2); (3.47) 2d2goutdY2 +l3=2 5=2dfoutdY   4d2foutdY2  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 = p Y (3.49)for both equations. We will take the negative root as its integral is consistentwith the de nition of the phase function  (Y) = 23jYj3=2. At the next order, theequations contain both sine and cosine terms, and their coe cients must vanish forthe equations to hold. We can therefore extract four equations:2dR1dY d dY +R1 d2 dY2 = 2 S2 (3.50)2dR2dY d dY +R2 d2 dY2 = 2 S1 (3.51)2dS1dY d dY +S1 d2 dY2 =  1R2 l3=2d dY R1 d dY2R2 (3.52)2dS2dY d dY +S2 d2 dY2 = + 1R1 l3=2d dY R2 +d dY2R1 (3.53)These can be simpli ed by de ning  = 4p8 p Y, and ^Rn; ^Sn = Rnp ; Snp .42With d dY = p Y =  = 4p8 , d2 dY2 = 4p8 =(2 ), (3.50) to (3.53) simpli es to:^R1 = 2 4p8 ^S2 (3.54)^R2 = 2 4p8 ^S1 (3.55)^S1 = 14p8    1R2 +  4p8 l3=2R1  2p8 R2 (3.56)^S2 = 14p8   1R1 +  4p8 l3=2R2 +  2p8 R1 (3.57)Eliminating ^S1 and ^S2 from (3.54) to (3.57), and de ning A =  1p 2 and B =l3=22 4p 2, we have a pair of coupled second order ODE:^R1     24 +A ^R1 = B ^R2 (3.58)^R2     24 +A ^R2 = B ^R1 (3.59)It turns out that the pair can be uncoupled by de ning F1 = ^R1 + i^R2 and F2 =^R1 i^R2, the di erential equations for F1 and F2 are:F1    14( + 2iB)2 +B2 +A F1 = 0; (3.60)F2    14(  2iB)2 +B2 +A F2 = 0; (3.61)both of which have parabolic cylinder functions as solutions. For the outer solutionsfout(Y) and gout(Y) we demand that they decay to zero as Y ! 1, which isequivalent tojF1;2j!0 as  !1. The standard solutions to the parabolic cylinderequations are given by U(a;x) and V(a;x), and only U(a;x) decays as x ! 1(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 innersolution.43MatchingTo match up the inner and outer regions, the outer solution should take the formof (3.45) and (3.46) as Y ! 0. To facilitate matching, we will expand (3.45) and(3.46)fin  1=4  14p Y + 2 r4p Y sin  (Y) 3=2 + 14p Y  2 r4p Y cos  (Y) 3=2  ;(3.63)gin  3=4  r4p Y + 2 4p Y sin  (Y) 3=2 + r4p Y  2 4p Y cos  (Y) 3=2  ;(3.64)with the constant factors p2=2 being scaled out.Comparing the fin and fout from (3.45), we will need the coe cients of thesine and cosine terms to match as Y ! 0. Remembering that ^Rn = Rnp andp Y =  = 4p8 , it suggests that near Y = 0, the matching conditions are:^R1( ) 4p8 + 2 r and ^R2( ) 4p8  2 r ;=) ^R1(0) = ^R2(0) = 4p8 and ^R01(0) = ^R02(0) = 2 r: (3.65)Similarly, comparing gin and gout gives^S1( ) r 4p8   1 and ^S2( ) r 4p8 + 1 :=) ^S1(0) = ^S2(0) = r 4p8 and ^S01(0) = ^S02(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 satis ed as long as the boundary conditions in (3.65) aresatis ed. The conditions ^R1(0) = ^R2(0) and ^R01(0) =  ^R02(0) from (3.65) in turnrequire thatF1(0) = iF2(0) and F01(0) = iF02(0) (3.67)44or equivalently,264U(B2 +A;2iB)  iU(B2 +A; 2iB)iU0(B2 +A;2iB)  U0(B2 +A; 2iB)375264c1c2375 =26400375; (3.68)For the system (3.68) to have a non-trivial solution, the determinant of the matrixmust vanish:U(B2 +A;2iB)U0(B2 +A; 2iB) +U(B2 +A; 2iB)U0(B2 +A;2iB) = 0: (3.69)Now the Wronskian of U(a;x) and U(a; x) is given by(U(a;x))0U(a; x) U(a;x)(U(a; x))0 =U0(a;x)U(a; x) +U(a;x)U0(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 U0(B2 + A; x)must vanish at x = 2iB. Since both U(B2 +A;x) and U(B2 +A; x) are solutionsto the parabolic cylinder equationd2Udx2   14x2 +B2 +A U = 0; (3.71)the Wronskian must be a constant due to Abel’s identity1. Therefore if the Wron-skian ofU(B2 +A;x) andU0(B2 +A; x) were to vanish at x = 2iB, it must vanishat x = 0 as well. In other words,U(B2 +A;0)U0(B2 +A;0) = 0: (3.72)With the formulae forU andU0 at the origin from Bender & Orszag (1999), as wellas noting that that the Gamma function  (z) diverges when z is zero or a negative1Abel’s identity states that for a linear second order di erential equation y00 +p(x)y0 +q(x)y = 0, the Wronskian is given by C exp  R p(x)dx . For the parabolic cylinder equa-tion p(x) = 0, and hence the Wronskian is a constant.45integer, we have:U(B2 +A;0) = 0 () 34 + 12(B2 +A) = 0; 1; 2; 3   ;U0(B2 +A;0) = 0 () 14 + 12(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 n2N. Revertingback to the original variables, we  nd that the correction to the eigenvalue is: 1 =  2n 12 r2  l23=24 (3.73)The prediction for the eigenvalue is therefore: =  +  1 =   (2n 1)r 2  l24 (3.74)Therefore the fastest growing mode corresponds to n = 1, which in turn givesU0(B2 +A;0) =U0( 1=2;0) = 0.To complete our analysis, we will need to determine the constant r in theinner solution. Note that ^R1(Y) = (F1 + F2)=2, and with the aid of (3.68) we caneliminate c1, leaving one undetermined coe cient for ^R1. In addition, we have notused the conditions ^R1(0) = 4p8 and ^R01(0) = 2 r in (3.65) during our matching;after some algebra, these conditions will give:r = i4p8 2 U0( 1=2; 2iB)U( 1=2;2iB) : (3.75)Notice that U( 1=2;x) = e x2=4 satis es (3.71) (Temme, 2010), and thereforer = l3=2p2 : (3.76)463.2.3 EnergeticsTo explain the unexpected results observed in the previous section, we will examinethe energetics associated with the normal modes. If we substitute the normal modesolutions Aj(y)e  i(! a))t into the energy equations (3.4) and (3.5), we will have:2 jA1j2 = J1y D1 +S; (3.77)2 jA2j2 = J2y D2 +S: (3.78)Integrating over the entire domain, the equations become: =  ^D1 + ^S .2j^A1j2 =  ^D2 + ^S .2j^A2j2 ; (3.79)where the hatted variables are the integrated values over the entire domain. Inaddition we have assumed that Jn vanishes as y !  1 since jAnj ! 0. Theenergy equation for the normal mode problem simply states that the growth rateof the instability is determined by the di erence between energy source due to PSIactivity ( ^S) and energy sink due to dissipation ( ^D1). Figure 3.6 shows the resultsfor  = 1 and l = 0, where total energy source ^S and total dissipation ^D1 arecalculated using the most unstable eigenfuctions, and are normalized by the wave-energyj^A1j2. Not surprisingly, the dissipative term ^D1 (green dots) increases as thedissipation increases; at the same time, the energy source ^S (blue dots) also decreasessigni cantly (20 % compared to the inviscid case), indicating that dissipation alsodisrupts the energy transfer. At higher dissipation a signi cant portion of the energy( 50%) generated via PSI is dissipated away.The results for l = 0:1 is plotted in  gure 3.7. For larger dissipation, thesituation is very similar to the l = 0 case: The total dissipation ^D1 increases signif-icantly while the energy source ^S drops. On the other hand, as  ! 0 the system4710−410−310−210−110000.µˆS/2ˆ|A1|2ˆD/2ˆ|A1|2γRPredictionFigure 3.6: Energetics of PSI: l = 0. ^S (blue dots) and ^D1 (green dots) from(3.79) are calculated and normalized by using the eigenfunctions, and are plottedhere as a function of  . The numerically computed growth rate (red dots) andthe asymptotic prediction from (3.74) (dashed line) are also plotted for comparison. = 1 and yL = 200.4810−310−210−1100−0.4−µABCˆS/2ˆA1ˆD/2ˆA1γRPredictionFigure 3.7: Energetics of PSI: l = 0:1. ^S (blue dots) and ^D1 (green dots) from (3.79)are calculated and normalized by using the eigenfunctions, and are plotted here as afunction of  . The numerically computed growth rate (red dots) and the asymptoticprediction from (3.74) (dashed line) are also plotted for comparison. The verticalline separates the boundary modes and the localized modes.  = 1 and yL = 1500.behaves very di erently compared to the l = 0 case: when  drops below 10 2, theenergy source ^S drops o rapidly while the energy dissipation ^D1 increases. Thesetwo e ects combine to stabilize the system when the dissipation is su ciently 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 to1, Z1yS(s)ds =Z 1yD1(s)ds+ 2 Z 1yjA1(s)j2ds 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 and1.(which is represented by J1(y). We will consider the three modes that are close tothe 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  gure 3.7. The four terms in (3.80) are computed using theeigenfunctions (normalized such that max(jA1j) = 1) and are plotted in 3.8 togetherwith the eigenfunction.From the plot of jA1j2 in panel a, we can see clearly that the waves becomemore and more de-localized as  decreases; the eigenfunctions are very similar nearthe turning latitude, but the eigenfunction for  3 (blue curve) decays signi cantlyfaster than the eigenfunction for  1 (red curve). A plot of R1y jA1(s)j2ds (panelb) illustrates the delocalization more precisely: the energy density integral for  3becomes 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 ofenergy  ux J1 in panel c provides a possible explanation: as  decreases, the merio-dional energy  ux becomes more negative, meaning that more energy is transportedsouthwards by the PSI waves. As the waves become increasingly oscillatory to thesouth, dissipation is more e ective since it increases with the wave-number. Thusin this regime the energy generated near the critical latitude is transported to thesouth where it is dissipated away, and the end result is that the global dissipationincreases despite a decrease in the dissipation parameter.The other phenomenon we have observed earlier is that the overall energysource via PSI ^S also decreases with  . R1y S(s)ds is plotted in panel (e) of  gure3.8, which demonstrates how the source term changes over the domain. For  350(blue curve), it  rst increases rapidly as we move south and away from the turningpoint, indicating there is signi cant energy input near the critical latitude. For  2(green curve), the integral increases much slower near the critical latitude comparedto 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 pumpwave. The integral continues to decrease as we head south until y =  50, wherea minimum is reached and the integral starts to increase; thus the source term ispositive to the south of y = 50, and energy is transfered from the pump-wave intothe PSI waves. The net transfer of energy ^S is small as the drain of energy to thenorth is balanced by the gain of energy to the south.51−400 −300 −200 −100 00510152025integraltext∞y|A1|2dy(b)−400 −300 −200 −100 0−1.5−1−0.500.5J1(c)−400 −300 −200 −100 002468integraltext∞yD1dyy(d)−400 −300 −200 −100 0−2−10123integraltext∞ySdyy(e)−100 −90 −80 −70 −60 −50 −40 −30 −20 −10|A1|2y(a)µ =0.002238µ =0.002512µ =0.002818Figure 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 computethe various terms in the energy equation (3.80). (a) jA1j2. (b) R1y jA1(s)j2ds. (c)J1. (d)R1y D1(s)ds. (e) R1y S(s)ds.  = 1 and yL =  1500. Please note that therange for the horizontal axis in (a) is  100 <y < 5 but the range for (b) to (e) is 400 <y< 5.523.3 More initial value problemsA couple of Initial value problems will now be presented to further illustrate someof the properties of our PSI model discussed in x3. Transient growthFigure 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 asymptoticallystable, 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(p26y); A2(y;0) = 0 (3.81)In the top panel of  gure 3.9,jA1jis plotted for 0 t 2:5 and 80 y 5. As inthe example shown in x3.1, the wave packet splits into two components that travelin opposite directions; between the wave packets, the wave amplitude experiencesa growth due to PSI. This growth continues exponentially, and the wave amplitudeincreases from O(1) to O(105) by t = 20 (see the second panel); the growth comesto an end at t  20, and the wave amplitude subsequently decays exponentiallyat a rate that is consistent with the eigenvalue analysis. In  gure 3.10 the waveamplitude jA1j at various times are plotted, which illustrates the transient growthmore clearly.The numerical result suggests that a very signi cant transient growth ispossible, even though the system is asymptotically stable. We will be examiningtransient growth in much greater detail in x4.3.2.533.3.2 Boundary modesThe second example is one where a boundary mode eventually dominates the nu-merical solution. In  gure 3.11, the amplitude of the wavejA1jis plotted at varioustimes for  = 2, l = 0:1 and  = 10 4; the computational domain is similar tothe previous example. Panel (a) shows the initial perturbation given by a Gaussianwave packet centred at y = 75:A1(y;0) = e (y+75)2=5 cos(p75y); A2(y;0) = 0: (3.82)The initial dynamics is similar to the previous example, where an exponentiallygrowing mode is excited between the wave packets (panels (b) to panel (f)); thedi erence however is that the mode eventually gives way to an exponentially growingboundary mode that dominates the solution (panels (g) to (l)). Note that the wavepro le appears to be solid due to the rapid oscillations, which are well resolvedby the grid used in the numerical computation. As in the normal mode anlysis inx3.2.1, the type of boundary conditions (zero dirichlet, sponge layer or radiative)had minimal e ect on the growth rate of the unstable boundary modes in the IVP.54Figure 3.9: An initial value problem with transient growth The results from a nu-merical computation of the PSI model. jA1j is plotted in the top two panel, withthe top panel showing the initial growth from t = 0 to 2.5. The second panel showsthe entire computation from t = 0 to 40. Note that the color scale is di erent forthe two plots. A similar set of plots forjA2jare produced in the bottom two panels.The parameters used are  = 2, l = 0:5, and  = 0:028, with the spatial domainbeing [ 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 2000.51t=0(a)−80 −60 −40 −20 0 2000. t=0.75(b)−80 −60 −40 −20 0 2000.20.4t=1.5(c)−80 −60 −40 −20 0 2000.51t=2.25(d)−80 −60 −40 −20 0 200510t=4(e)−80 −60 −40 −20 0 20050100150200t=6(f)−80 −60 −40 −20 0 20050010001500t=7.5(g)−80 −60 −40 −20 0 200123x 104t=10(h)−80 −60 −40 −20 0 2002468x 105t=17.5(i)−80 −60 −40 −20 0 2001234x 105t=22.5(j)−80 −60 −40 −20 0 200510x 104t=27.5(k)−80 −60 −40 −20 0 2002000400060008000t=35(l)Figure 3.10: An initial value problem with stable normal modes. jA1j at varioustimes are plotted. The parameters used are  = 2, l = 0:1, and  = 0:028, with thespatial domain being [ 185;50] and a sponge layer located at [ 200; 85]. Only thedata from [ 80;30] is shown here.56−200 −150 −100 −50 000.51t=0(a)−200 −150 −100 −50−200 −150 −100 −50−200 −150 −100 −50 00123t=2.25(d)−200 −150 −100 −50 00200400600t=5(e)−200 −150 −100 −50 002468x 108t=12.5(f)−200 −150 −100 −50 002468x 1012t=17.5(g)−200 −150 −100 −50 001234x 1016t=22.5(h)−200 −150 −100 −50 00123x 1018t=25(i)−200 −150 −100 −50 0012x 1020t=27.5(j)−200 −150 −100 −50 000.511.52x 1022t=30(k)−200 −150 −100 −50 000.511.52x 1026t=35(l)Figure 3.11: An initial value problem with boundary mode. jA1j at various timesare plotted here. The parameters used are  = 2, l = 0:1, and  = 10 4, with thespatial domain being [ 185;50] and a sponge layer located at [ 200; 85]. Only thedata from [ 200;10] is shown.573.4 SummaryIn this chapter we have presented an initial value problem in which a perturbationin form of a Gaussian wave packet results in exponential growth, and a normalmode analysis con rms the existence of unstable modes. The stability problem wasinvestigated numerically, as well as analytically in the limit of small dissipation andsmall pump wavenumber (i.e. l 1 and   1). When l is non-zero, the unexpectedresult is a vanishing dissipation has a stabilizing e ect, which is con rmed by aninitial value problem; however the IVP also showed that a signi cant transientgrowth is possible even when the PSI-waves are asymptotically stable. In the nextchapter we will explore the transient growth, as well as the stability of the PSI-modelin the inviscid limit.58Chapter 4PSI in the Inviscid LimitThe 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 previouschapter, as we decreased the dissipation parameter ( ), the spectra eventuallybecame contaminated by unstable boundary modes that have little physicalmeaning. While the analytical and numerical results suggests that the discreteunstable modes will eventually become stable as  !0, it is not clear if thereare other physical modes present.2. What is the nature of the large transient growth observed in some initial valueproblems? We are interested in understanding the mathematical explanationbehind the transient growth in an asymptotically stable system.594.1 Stability of the inviscid problemTo address the  rst question, we will solve the inviscid problem analytically via aFourier 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 transformFf gto the equations, and de ningF(q) =FfA1(y)e ily=2gand G(q) =FfA2(y)eily=2g, we havedFdq = i q2 + 14l2 +ql i  F + i G; (4.3)dGdq = i F + i q2 + 14l2 ql + i  G: (4.4)Note that under Fourier transform, y is transformed into i ddq. These equations canbe further simpli ed by a change of variable:^F(q) = ei 13q3+14l2q F(q) and ^G(q) = ei 13q3+14l2q G(q): (4.5)(4.3) and (4.4) become:d ^Fdq = i(ql i )^F + i ^G; (4.6)d ^Gdq = i ^F i(ql i ) ^G: (4.7)We can eliminate one of the equations to form a second order ODE. For example,using (4.7) to eliminate ^G in (4.6), we have:d2 ^Fdq2 + (lq i )2 + 2 il ^F = 0: (4.8)De ning z qil2 q i l and H(z) = ^F(q) the equation can be manipulated intothe standard parabolic cylinder equation:d2Hdq2   14z2 +  i 22l  12  H = 0: (4.9)60In this case the solutions are U(a;z) and U(a; z), where a =  i 22l  12; theirWronskian is given by (Temme, 2010):W =p2  (12 +a) =p2  ( i 22l ) ; (4.10)which will not vanish since the gamma function is  nite along the imaginary axis.Thus these two solutions to (4.14) are linearly independent, and form a satisfactorybasis.To carryout the inverse Fourier transform, H must remain bounded as q ! 1; however we can demonstrate that this cannot be true if  R6= 0.We will  rst consider q!1. In this case argz  =4, and the asymptoticformula for jzj 1 is given by (Temme, 2010):U(a;z) e 14z2z a 1=2 +O(1=z2): (4.11)Now sincez2 = il2 q2 2i l q  2l2 = il2 q2  2l2 + q; (4.12)e z2=4 will be bounded as q !1 if and only if  R  0. On the other hand, asq! 1, argz  3 =4, and the appropriate asymptotic formula isU(a;z) e 14z2z a 12  ip2  (12 a)ei ae14z2za 12 +O(1=z2): (4.13)In this case the  rst 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 forU(a;z) to be bounded. For the other solution U(a; z), a similar analysis will leadto the same conclusion.Therefore for if H is to remain bounded as q ! 1,  R = 0. What theanalysis above suggests is that for the inviscid problem, there are no unstable normalmodes. 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 thecorresponding eigenfunction can be obtained by reversing the Fourier transform.614.2 Non-normal operators and transient growthTraditionally in  uid mechanics, the hydrodynamical stability is considered throughnormal mode analysis, where in nitesimally small perturbations can grow exponen-tially into  nite amplitude disturbances for unstable  ows. The growth rates forthe normal modes are given by the eigenvalues, while the spatial structure is de-scribed by the eigenfunctions. Normal mode analysis have been quite successful forproblems such as the Rayleigh-B ernard instability and Taylor-Couette instability, asexperimental and analytical prediction for instability agree (Trefethen et al., 1993).However for shear  ows such as plane Couette and Poiseuille  ows, the onset ofturbulence occurs at Reynolds number much lower than the threshold predicted bytheory (Trefethen et al., 1993; Butler & Farrell, 1992). For example, Couette  owsare theoretically stable for all Reynolds numbers, but transition to turbulence wasobserved in laboratory at Reynolds number of 360 (Tillmark & Alfredsson, 1992).For plane Poiseuille  ow, the critical Reynolds number unstable normal modes forthe Orr-Sommerfeld equation is Re > 5772:22 (Orszag, 1971), while transition tofull 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 graduallyrecognized. It turns out that for a problem where the operator is ‘normal’ (i.e. theeigenfunctions are orthogonal), the behaviour of the system is completely determinedby the spectra. The linear operators in the Rayleigh-B ernard and Taylor-Couettestability problems fall into this category, and therefore the normal mode analysissucceeded 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 thesense that the eigenfunctions are not orthogonal. A consequence of non-orthogonaleigenfunctions is that while the long term behavior of a system is dictated by the62spectra, these systems can have signi cant transient growth, which the normal modeanalysis fails to capture.4.2.1 A simple exampleThe following example was used by Farrell & Ioannou (1996) to illustrate some ofthe features of non-normal operators. Consider a simple 2 2 matrix:A =264 1  cot 0  2375; (4.14)A matrix M is de ned to be normal if MMy = MyM, wherey denotes the Hermitiantranspose1. Here we haveAAy AyA =264cot2 cot cot2 0375; (4.15)and therefore A is normal only if  =  =2, and becomes increasingly non-normalas  ! 0. It turns out that the non-normality can play an important role in thetransient behaviour of a dynamical system. Farrell & Ioannou (1996) considered thefollowing dynamical system:dxdt = Ax: (4.16)The matrix A has eigenvalues r = 1 and r = 2 with eigenvectors v1 = (1;0) andv2 = (cos ;sin ), and the general solution is given by:x = c1e t26410375+c2e 2t264cos sin 375: (4.17)Consequently all solutions are asymptotically stable (i.e. x! 0 as t!1 for allinitial conditions x0), with decay rate  1. Farrell & Ioannou (1996) continued on1A Hermitian transpose is also known as the conjugate transpose63to demonstrate that an asymptotic stability analysis does not capture the dynamicsof the system as  !0. In this limit, v2!v1, and the eigenvectors become almostparallel to each other. The rami cation is that if the initial condition is almostorthogonal to these eigenvectors, c1 and c2 have to be large and have opposite signsin order for the initial condition to be satis ed. For example, the solution withinitial condition x0 = (sin ; cos ), which is orthogonal to v2, will be given byc1 = sin + cos cot = csc and c2 =  cot ; clearly in this case as  ! 0, c2decreases without bound, and c1 ! c2. The solution to the initial value problemisx = e t264csc 0375+e 2t264sin  csc  cos 375: (4.18)Now as t increases, the second term will quickly vanish, while the  rst term canhave a magnitude much greater than unity if  is small. As jx0j= 1, this indicatesthat the amplitude experiences a signi cant growth before decaying. To furtherillustrate this, notice that the solution to (4.16) with initial condition x(0) = x0 canbe written as:x = eAtx0; (4.19)If we restrict ourselves to initial conditions that satis es jx0j= 1, thenmaxjxj= maxjx0j=1jeAtx0j=jjeAtjj; (4.20)where jj jj denotes the induced matrix 2-norm. jjeAtjj is therefore a measure of themaximum magnitude of x. In  gure 4.1, jjeAtjj for  =  =2, =4, =10, and  =100are plotted as a function of time. For all four cases, the curves all become a straightline with slope -1 in accordance with the asymptotic stability analysis. The maindi erence is that for  =  =2, which corresponds to A being a normal matrix, thesolution decays with a decay rate of -1 for all t. On the other hand, as  decreases640 0.5 1 1.5 2 2.5 3 3.5 410−1100101t||eAt||θ = pi/2 θ = pi/4θ = pi/10θ = pi/100Slope = −2Slope =12[−3 + csc(pi/100) ]Figure 4.1: Transient growth for a simple dynamical system. jjeAtjjis plotted for fourvalues of  :  =  =2, =4, =10, and  =100, where A is de ned in (4.14). Also plottedas dash line is the initial growth rate give by 12( 3+csc ), and the asymptotic decayrate -1 for  =  =2. Note that the vertical axis is in logarithmic scale.to  =4 and  =10, A becomes increasingly non-normal, and the decay in jjeAtjj isdelayed. For  =100, the di erence is even more remarkable, as jjeAtjj increasedalmost an order of magnitude before decaying. The initial growth rate near t = 0 isgiven by the maximum eigenvalue of (A+Ay)=2, which in this case is 12( 3 + csc )(Farrell & Ioannou, 1996). This suggests that for positive initial growth the criterionis  <  =9:244. From  gure 4.1, we can see that for  =  =10 <  =9:244, jjeAtjjindeed increases slightly near t = 0. For  =  =100, the curve for jjeAtjj agrees wellwith the predicted initial growth rate of 12( 3 + csc =100).654.2.2 Non-normal operators and  uid dynamicsThis simple example in the previous section illustrates that while the dynamics ofa normal system can be characterized by the eigenvalues, they do not capture thetransient behaviour for non-normal systems. This fact was demonstrated clearly ina study by Reddy et al. (1993) where the Orr-Somerfeld operator in the Poiseuilleproblem was analyzed; for a Reynolds number in the stable regime, the energyincreased by a factor of 51 before the perturbation decayed away. Similarly, Butler& Farrell (1992) investigated the optimal excitation of perturbations in a planeCouette  ow, and found energy ampli cations of O(1000) despite Couette  ow beingstable. The implication of these  ndings is that transient growth can amplify smallperturbations and trigger the transition to turbulence without the need to normalmode instability. Non-normal operators also received attention in the geophysical uid dynamics community in recent years, for example in the studies of baroclinicinstability, error ampli cation in weather forecast, El Ni~no, and the stability ofthermohaline circulation (Farrell & Ioannou, 1996; Trefethen & Embree, 2005, Ch.23, and references therein).4.3 Pseudospectra and PSI4.3.1 PseudospectraProperties of a non-normal operator/matrix can be characterized by the pseudospec-tra. For   0, the  -pseudospectra of matrix A,   (A), is de ned as the set of z2Csuch that z is an eigenvalue of A+E for some perturbation matrixjjEjj  (Reddyet al., 1993). For normal operators, a perturbation to the matrix A of sizejjEjj  will perturb the spectra of A,  (A), by a distance no more than  . On the other66hand, for a non-normal matrix the eigenvalues can be perturbed much larger thanthe perturbation to the matrix (Reddy et al., 1993). The pseudospectra has to becomputed numerically for most cases, and it provides a convenient way to visualizethe degree of non-normality of an matrix operator.One important information that can be extracted from the pseudospectra isthe lower bound to the maximum ofjjeAtjj. Suppose that   (A) is the supremum ofthe real part of   (A), then the following inequality holds:supt 0jjeAtjj   (A)= : (4.21)Since (4.2) holds for all   0, a good estimate for the lower bound is the Kreissconstant, which is de ned asK(A) = sup  0  (A)= ; (4.22)which leads tosupt 0jjeAtjj 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 growthWe will now consider our stability problem in the context of pseudospectra. Re-writing (3.1) and (3.2),A1t = ih (1 i )A1yy +yA1 + eilyA2i; (4.24)A2t = ih(1 +i )A2yy yA2  e ilyA1i: (4.25)(4.24) and (4.25) can be written as a dynamical system:d~Adt = iM~A; (4.26)67where M is the matrix resulting from a spatial discretization. In this section we aregoing 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 matrixis dense compared to the  nite di erence scheme used in chapter 3, its smaller sizeallow for the use of eig in MATLAB to compute the full set of eigenvalues; inaddition 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 samemethod, and they agree over the region of interest.Inx3.3.1, we have seen that for  = 2, l = 0:5, and  = 0:028, the PSI wavesgrew signi cantly before decaying as predicted by the normal mode analysis, and wewould like to investigate this transient growth with the aid of the pseudospectra as-sociated with the spatial operator. The pseudospectra for the corresponding matrixM is calculated using the eigtool package2 for MATLAB. Due to the extra factorof  i in (4.26), we will have to consider the imaginary part of the spectrum andpseudospectra for information pertaining to growth. In  gure 4.2, we can see thatall eigenvalues (black dots) lie below the real axis, and therefore the normal modesare stable, con rming the long term behaviour exhibited by the IVP in x3.3. Onthe other hand, the pseudospectra protrudes signi cantly 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) shiftin the eigenvalues. Since the operator is non-normal, the transient growth we haveobserved in the IVP should not be surprising.We will now take a closer look at the transient growth in the corresponding2Thomas 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.500.51−14−13−12−11−10−9−8−7−6−5−4−3−2Figure 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) aslevel 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) =Z 1 1(jA1j+jA2j)dy (4.27)In  gure 4.3, Etot (normalized so that the Etot(0) = 1) is plotted as a function oftime (blue curve). We can see that the energy increased to almost 1013 times theinitial energy before the decay sets in at a rate predicted by the eigenvalues (dashedblack line). The magnitude of growth is consistent with the lower bound given by(4.21): from  gure 4.2, the  -pseudospectra for  = 10 7 have a maximum imaginarypart of 0.187, and therefore the lower bound for the maximum growth ofjjeAtjjwill69be 1:87 106, which is equivalent to an increase in energy by a factor of 3:5 1012.This estimation is plotted in  gure 4.3 as the blue dotted line, which shows that theobserved ampli cation is in line with the prediction given by the pseudo spectra.In  gure 4.3, Etot for  = 0:033 (green curve) and  = 0:038 (red curve) arealso plotted for comparison. The dissipation parameter  appears to only have aminor e ect on the transient growth rate; the onset of transient growth is delayedslightly for  = 0:038, but in all three cases the transient growth is exponential withcomparable growth rates. The transient growth continues pass t > 15, before thegrowth/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 asigni cant role in the evolution of the PSI-waves, particularly in cases where thenormal modes are neutral or stable; even when unstable normal modes are present,the transient growth can still occur at a signi cantly shorter time scale comparedto the asymptotic growth.We will present a plot to illustrate the dependence of transient growth onthe initial condition. In  gure 4.4, we have plotted Etot for  ve di erent initialcon gurations (see table 4.1). We have not explored in depth the optimal excitationthat leads to the maximum transient growth, but from the  gure we can notice thatthe wave-number of the perturbation a ects the transient growth signi cantly.The physical signi cance of our result is that transient growth, rather thanunstable 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 asymptoticallystable, the transient growth may be su ciently large for a transition into a non-700 5 10 15 20 25 30 35 40100105101010151020tEtotµ=0.028µ=0.033µ=0.038Figure 4.3: Transient growth in a time dependent problem. The total energy for theIVPs considered inx?? are plotted as a function of time. The black dashed lines onthe right are the asymptotic growth/ decay predicted by the eigenvalues, while theblue dotted line is the minimum growth predicted by the pseduospectra.linear regime.710 2 4 6 8 10 12 14 16 18 2010−2100102104106108101010121014tEtot  IC1IC2IC3IC4IC5Figure 4.4: Initial condition and transient growth. The total energy Etot for di erentinitial conditions are plotted as a function of time.  = 2, l = 0:5, and  = 0:028.Table 4.1: Initial conditions for  gure 4.4A1(y;0) A2(y;0)IC1 e (y+25)2=5 cos(p5y) 0IC2 e (y+25)2=5 cos(p5y) ie (y+25)2=5 cos(p5y)IC3 e (y+25)2=5 cos(p40y) 0IC4 e (y+25)2=5 cos(p15y) 0IC5 e (y+25)2=5 cos(p25y) 0724.3.3 Numerical di cultiesThe non-normality in operators/matrices also results in complications during nu-merical computations of eigenvalues, since the operator is sensitive to errors; theround o errors act as a perturbation to the matrix, and if the operator is non-normal, the perturbation is su cient to shift the eigenvalues. An example of theround-o errors can be detected in  gure 4.2; near the center of the plot where thethree 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 operatoriRed2dy2 +y; (4.28)as a model for the Orr-Sommerfeld operator in shear- ow instabilities. The spectraof the Airy operator consists of three branches in a ‘Y’-shape, which resembles thespectrum in  gure 4.2. Reddy et al. found the operator to be non-normal, and theeigenvalues near the intersection were di cult to compute as they were sensitive toperturbations; 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!1 in Reddy et al. corresponds to  ! 0 in our model. Thus in hindsight, we should not be surprised by the non-normality of our eigenvalue problemin the inviscid limit.The e ect of round-o error is much more evident in  gure 4.5, where theeigenvalues for discretized operator with  = 6, l = 1:5 and  = 0:1 are computedusing two di erent level of precision. In the right panel, the eigenvalues are calcu-lated using double precision  oating point, and we can see that the eigenvalues lineup along the boundary of the pseduospectra with  = 10 14. In the left panel, the73eigenvalues for the same matrix were computed with single precision; the eigenvaluesnear 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 forthe two computations, showing that a change in numerical precision did not a ectthe accuracy signi cantly for these parts of the spectra.When the non-normality becomes severe, the round-o error prevents theeigenvalues near the center from being calculated accurately. The growth rate pre-dicted by our asymptotic theory inx3.2.2 (c.f. (3.74)) for  = 6, l = 1:5 and  = 0:1is 0.346. In panel (b) of  gure 4.5, the predicted eigenvalue of 0.346i is right inthe middle of the area enclosed by the eigenvalues perturbed due to round-o error.Due to the severe non-normality, the true eigenvalue cannot be resolved even by theeigs command in MATLAB.This e ect of non-normality is the most prominent in the transition regionbetween the boundary modes in the inviscid limit to the localized modes. Theexample considered above is marked in panel (b) of  gure 3.5 by A. As the largest R reported by the eigenvalue calculation near A is the perturbed eigenvalue due tothe round-o error rather than the true eigenvalue, the growth rate is potentiallyoverestimated. This explains why the numerically calculated growth rate agreeswith the prediction for larger values of  , but they break o suddenly from theprediction when  becomes too small. By computing the pseudospectra near thebreaking point for decreasing values of  , we have con rmed that the round o errorstarts a ecting the eigenvalue with the largest growth rate precisely at the breakingpoint in the plot.The non-normality of the matrix also a ected the use of iterative schemes,such as the eigs function in MATLAB, for the computation of eigenvalues. In74−20 −15 −10 −5 0 5 10 15 20−2−101234(a) Single Precision−20 −15 −10 −5 0 5 10 15 20−2−101234−14−13−12−11−10−9−8−7−6−5−4−3−2−1(b) Double PrecisionFigure 4.5: E ect of non-normality on the calculation of eigenvalues The eigenvaluescomputed for  = 6, l = 1:5 and  = 0:1 for (a) single precision  oating point, and(b) double precision  oating point.particular, when the eigenvectors become almost linearly dependent, the convergencein such iterative schemes can be signi cantly slowed (Trefethen & Embree, 2005);this has been a problem for our calculations, especially in the invisicid limit. Anadditional complication is that even when the iterations converge, the estimatedeigenvalues (Ritz values) are not guaranteed to be close to the actual eigenvalues(Trefethen & Embree, 2005). Indeed in the inviscid limit we encountered di cultyin using eigs to compute the eigenvalues as the convergence was slow.4.4 SummaryIn this chapter we considered the stability of the inviscid PSI problem, and concludedthat there are no unstable normal modes. In addition, as  !0, the spatial operatorin our model becomes non-normal. A consequence of the non-normality is that asigni cant transient growth can occur with suitable perturbations, and the transient75growth 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 thetraditional modal instability.76Chapter 5Discussion5.1 Horizontal dissipationIn this thesis, we have formulated a linear model for PSI, and explored its behaviourin Chapter 3. Solving the PDE numerically as an initial value problem revealedthe presence of exponentially growing normal modes, which suggest that the PSIproblem can be formulated in a normal mode stability analysis. The main result isthat for a meridionally propagating pumpwave (i.e. l6= 0), instability can only occurwhen the horizontal eddy viscosity  H is su ciently large. Through the analysisof the energy equation, the stability in the inviscid limit can be attributed to twofactors. A reduction in  H increases the meridional energy  ux towards the equator,which prevents a build up of energy in the PSI waves near the critical latitude. Atthe same time the energy transfer from the pump wave to the PSI waves is also lesse ective.The  -e ect 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 inCoriolis force to the south implies that the horizontal wave number must increase,77and therefore the wave energy is dissipated away more rapidly to the south. Thisprovides an alternative mechanism to the localization of PSI seen in MacKinnon& Winters (2005), where the localization was attributed to a vertical divergence inenergy  ux preventing a build up of energy to the south of the critical latitude.Depending on the horizontal scale, the horizontal eddy di usivity measuredfrom dye-release experiments ranges from 10 2 m2s 1 for a horizontal scale of 0.1kmto 103m2s 1 at a scale of 1000km (Ledwell et al., 1998). The intermediate values of100 101m2s 1 correspond to   10 2 10 1 since the eddy di usivity is scaledby  3 = 1:36 102m2s 1. For this range of  , our model predicts that the PSIactivity will be localized to between 15-30 units, or 300-600km, to the south of thecritical latitude (c.f.  gure 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 to22 N, or approximately between 750 to 1000 km south of the critical latitude. Thusthe meridional extent of PSI imposed by horizontal eddy di usion can be signi cantcompared to other length scales, such as the source location for the pump-wave.For an near-intertial wave, we can demonstrate that the horizontal eddydi usivity is signi cant compared to vertical eddy di usion. For a plane wave, thehorizontal and vertical eddy di usion can respectively be characterized by H(k2+l2)and  Vm2. We wish to consider the ratio between the two quantities. Using thedispersion relation for an internal wave with the near-inertial approximation ! = f0,f = f0 + y, and N f H(k2 +l2) Vm2 = H V!2 f2N2 f2   H V 2f0 yN2 : (5.1)Taking  V = 10 5m2s 1,  H = 10m2s 1 and the parameters from table 2.1, the78last term in (5.1) is: H(k2 +l2) Vm2 = y27km (5.2)What this tells us is that near the origin, the vertical dissipation is much largerthan the horizontal dissipation, which is expected since near the critical latitude, thevertical structure of near-inertial waves is much  ner than the horizontal structure;however as the wave packet moves south, the horizontal wave number increases, andat 27km south of the critical latitude the horizontal dissipation is the same orderof magnitude as the vertical dissipation. This justi es the include of the horizontaldissipation in our model.In summary, while horizontal dissipation was included originally as a curefor numerical issues (boundary modes), it is clear that it is in fact dynamicallysigni cant 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 eigenvalueproblem with  = l = 0, where a continuous spectrum of unstable modes are possi-ble when  > 0. On the other hand, while a continuous spectrum of eigenfunctionsalso exists for the inviscid problem with  = 0 but l6= 0, the eigenfunctions are allstable. This suggests that  = l = 0 is a singular limit as we cannot recover thesolutions by taking l!0. The inviscid solutions appear to be singular as well sincethe growth rate of the discrete normal modes in the dissipative problem tends to 1 as  !0.795.2 Transient growthAnother interesting mathematical property of our model is the non-normality ofthe spatial operator. As we approach the inviscid limit, the eigenfunctions for theoperator become parallel to each other, which caused di culties in our numericalcomputations for eigenvalues. A closer look at our spatial di erential operator re-vealed connections to the Orr-Sommerfeld operator for the stability problem in planePoiseuille and Couette  ows, 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 havedemonstrated that transient growth is indeed possible, despite the system beingasymptotically 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 unstablenormal 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 theregime where the growth rate increases with dissipation (c.f. panel (a) of  gure 3.3).For  = 0:05, the normal modes are all stable, yet the pseudo-spectra indicates thepossibility of large amplitude transient growth (see  gure 5.1). The implication isthat the transient growth we observed may be oceanographically relevant; transientgrowth may be the mechanism via which a transition into non-linear regime takesplace, rather than via normal mode instability.80γIγR−8 −6 −4 −2 0 2 4 6 8−1−0.500.51−8−7−6−5−4−3−2−1Figure 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 eigenvaluesare plotted as black dots. The operator is discretized using a spectral collocationmethod using Chebyshev polynomials as interpolant over the domain [ 75; 75]5.3 Outstanding issuesIn this thesis we have established that dissipation plays an important role in thestability of PSI. However a number of key questions remain:1. What are the implications of a realistically strati ed ocean?Through-out this thesis we have assumed a linearly strati ed ocean, and thevertical structure of the PSI wave and that of the pump wave are assumed tobe independent. Further more, the pump-wave is assumed to be a plane wave.81A more realistic pro le for the ocean is the Gill’s pro le, where the oceanis modeled by a mixed layer with uniform density near the surface, and thedensity increases smoothly below the mixed layer (Gill, 1984). Mathematically,the buoyancy frequency is given by:N(z) =8><>:0 if 0 z hmixN0=(z z0) if z hmix(5.3)The internal-tides associated with this pro le can be solved analytically, andthese can be used in place of the planar pump-wave in our model. YTBinvestigated PSI with the Gill’s pro le on thef-plane, and found intensi cationof PSI at the base of the mixed layer. On the other hand, the measurements ofAlford et al. (2007) found PSI activity extending much deeper into the ocean.By generalizing our model to a realistically strati ed ocean, we can perhapsbetter reproduce the observations.2. How is the vertical scale of the PSI-waves determined?Another aspect of PSI that demands clari cation is the vertical scale selection.The numerical results of MacKinnon & Winters (2005) and the  eld data ofAlford et al. (2007) both suggest that PSI favours a particular vertical scale.In the current model the vertical wavenumber m is a scaling parameter, andtherefore 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 signi cant transient growth. We have also demon-strated that the transient growth depends strongly on the initial perturbation.82Therefore 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 isclear in our model that perturbations are either transiently or asymptoticallyunstable, and small perturbations can be ampli ed to the point where a con-stant amplitude for the pump-wave is no longer valid. In addition the mainmotivation for the study of PSI is to understand the energy transfer out ofinternal tides, which requires a model with genuine three-mode coupling thatis 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-wavecan be derived if we allow the background variables to vary on the slow timescale. 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 thenumerical model of MacKinnon & Winters (2005), and the more moderate lossmeasured 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 oceanographically83relevant, one must study non-linear e ects such as saturation of the instabilityand the feedback on the pump-wave to fully understand the PSI in the ocean.84BibliographyAbramowitz, Milton & Stegun, Irene A. 1964 Handbook of mathematicalfunctions with formulas, graphs, and mathematical tables, National Bureau ofStandards Applied Mathematics Series, vol. 55. For sale by the Superintendent ofDocuments, U.S. Government Printing O ce, Washington, D.C.Alford, M. H, Mackinnon, J. A, Zhao, Zhongxiang, Pinkel, Rob, Kly-mak, Jody & Peacock, Thomas 2007 Internal waves across the paci c. Geo-physical Research Letters 34 (24), 6.Bender, Carl M. & Orszag, Steven A. 1999 Advanced mathematical methodsfor scientists and engineers. I . New York: Springer-Verlag, asymptotic methodsand perturbation theory, Reprint of the 1978 original.Berenger, J.P 1994 A perfectly matched layer for the absorption ofelectromagnetic-waves. Journal of Computational Physics 114 (2), 185{200.Bretherton, F. P 1964 Resonant interactions between waves. the case of discreteoscillations. Journal of Fluid Mech 20, 457{479.Butler, KM & Farrell, B 1992 Three-dimensional optimal perturbations inviscous shear  ow. Physics of Fluids A 4 (8), 1637{1650.Carlson, D. R., Widnall, S. E. & Peeters, M. F. 1982 A  ow-visualizationstudy of transition in plane Poiseuille  ow. Journal of Fluid Mechanics 121, 487{505.Carter, G & Gregg, M 2006 Persistent near-diurnal internal waves observedabove a site of m 2 barotropic-to-baroclinic conversion. Journal of PhysicalOceanography 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 laserranging- a continuing legacy of the apollo program. Science 265 (5171), 482{490.85Egbert, GD 1997 Tidal data inversion: Interpolation and inference. Progress inOceanography 40 (1-4), 53{80.Egbert, GD & Ray, RD 2000 Signi cant dissipation of tidal energy in the deepocean inferred from satellite altimeter data. Nature 405 (6788), 775{778.Egbert, GD & Ray, R 2001 Estimates of m 2 tidal energy dissipation fromtopex/poseidon altimeter data. Journal of Geophysical Research 106 (C10),22475{22502.Farrell, B & Ioannou, P 1996 Generalized stability theory. part i: Autonomousoperators. Journal of the Atmospheric Sciences 53 (14), 2025{2040.Garcia, Alejandro L. 1999 Numerical Methods for Physics (2nd Edition). UpperSaddle 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 progressreport. Journal of Geophysical Researchophys. Res. 80 (3), 291{297.Gill, A. E. 1984 On the behavior of internal waves in the wakes of storms. Journalof 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). Proceedingsof the Royal Society of London. Series A, Mathematical and Physical Sciences299 (1456), 77{103.Hibiya, T, Nagasawa, M & Niwa, Y 2002 Nonlinear energy transfer within theoceanic internal wave spectrum at mid and high latitudes. Journal of GeophysicalResearch-Oceans 107 (C11), 3207.Hibiya, T, Niwa, Y & Fujiwara, K 1998 Numerical experiments of nonlinearenergy transfer within the oceanic internal wave spectrum. Journal of GeophysicalResearch 103 (C9), 18715{18722.86Jeffreys, H 1921 Tidal friction in shallow seas. Philosophical Transactions of theRoyal Society of London. Series A 221, 239{264.Laurent, L St. & Garrett, C 2002 The role of internal tides in mixing the deepocean. Journal of Physical Oceanography 32, 2882{2898.Ledwell, JR, Watson, AJ & Law, CS 1993 Evidence for slow mixing across thepycnocline 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: Signi cant lossof low-mode tidal energy at 28.9 degrees. Geophysical Research Letters 32 (15),L15605.McComas, C & Bretherton, F 1977 Resonant interaction of oceanic internalwaves. Journal of Geophysical Research-Oceans 82 (9), 1397{1412.Miller, G 1966 The  ux of tidal energy out of deep oceans. Journal of GeophysicalResearch 71 (10), 2485{2489.M uller, 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 windmixing. Deep-Sea Res Pt I 45 (12), 1977{2010.Olbers, D & Pomphrey, N 1981 Disqualifying two candidates for the energybalance of the oceanic internal wave  eld. 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  nite amplitude.i. the elementary interactions. J Fluid Mech 9, 193{217.Phillips, O. M 1961 On the dynamics of unsteady gravity waves of  nite amplitude.ii. local properties of a random wave  eld. J Fluid Mech 11, 143{155.Rainville, L & Pinkel, R 2006 Baroclinic energy  ux at the hawaiian ridge:Observations from the r/p  ip. Journal of Physical Oceanography 36 (6), 1104{1122.87Ray, R & Mitchum, G 1997 Surface manifestation of internal tides in the deepocean: Observations from altimetry and island gauges. Progress in Oceanography40 (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 modi cation and geographic redistribution ofthe semi-diurnal internal tide. Ocean Model 21 (3-4), 126{138.Taylor, G 1920 Tidal friction in the irish sea. Philosophical Transactions of theRoyal 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 ofStandards and Technology.Thorpe, S.A. 2005 The Turbulent Ocean. Cambridge University Press.Tillmark, N. & Alfredsson, P. H. 1992 Experiments on transition in planeCouette  ow. Journal of Fluid Mechanics 235, 89{102.Trefethen, L, Trefethen, A, Reddy, S & Driscoll, T 1993 Hydrodynamicstability 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. Thebehavior of nonnormal matrices and operators.Weideman, J. A. & Reddy, S. C. 2000 A matlab di erentiation matrix suite.ACM Trans. Math. Softw. 26 (4), 465{519.Young, W. R, Tsang, Y. K & Balmforth, N. J 2008 Near-inertial parametricsubharmonic instability. J Fluid Mech 607, 25{49.Zheng, C 2007 A perfectly matched layer approach to the nonlinear schr odingerwave equations. Journal of Computational Physics 227 (1), 537{556.88


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics

Country Views Downloads
China 24 19
United States 12 2
Japan 8 1
Canada 8 0
Russia 4 0
India 3 0
Sweden 2 0
Poland 1 0
City Views Downloads
Unknown 12 3
Beijing 11 0
Shenzhen 10 16
Ashburn 7 0
Tokyo 5 0
Saint Petersburg 4 0
Ann Arbor 2 0
Tianjin 2 1
Wilmington 2 0
Stockholm 2 0
Fort Worth 1 0
Halifax 1 0
Chandigarh 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items