NON-SYMMETRIC HOLMBOE WAVESBySusan Patricia HaighB. Math. (Pure & Applied Mathematics) University of Waterloo, 1988M. A. Sc. (Aerospace) University of Toronto, 1990A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF MATHEMATICSINSTITUTE OF APPLIED MATHEMATICSWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAOctober 1995© Susan Patricia Haigh, 1995In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)_____________________Department ofThe University of British ColumbiaVancouver, CanadaDate 3 Oc:-DE-6 (2188)AbstractWhen two flows of different velocity and density meet, a shear layer with a densitygradient is formed. Under certain conditions this flow can be unstable. A staticallystable stratified shear flow in which the density interface is much thinner than the shearlayer thickness can be linearly unstable to two modes with equal growth rates and equaland opposite phase speeds. The superposition of these two modes is called a Rolmboeinstability. This instability is oniy possible when the flow is symmetric about the center ofthe shear layer. We examine the effect of breaking this symmetry by allowing the centerof the density interface to be displaced with respect to the center of the shear layer. Thereare three major components to this study: linear stability analysis, nonlinear numericalsimulations, and comparison with laboratory experiments.Linear stability analysis is used to examine the effect of the density interface offset onthe overall stability of the flow. Both inviscid theory with piecewise linear backgroundvelocity and density profiles, and viscous theory with smooth background profiles areused. As in previous studies, it was found that the growth rate of one mode increases andthat of the other mode decreases as the density interface displacement is increased. Theprecise behaviour depends on the relative thickness of the density interface with respect tothe shear layer thickness. For inviscid theory with piecewise linear background profiles,we show that the initial perturbations must be two-dimensional. When the effects ofviscosity and diffusion are included, however, it may be possible that the weaker mode isinitially three-dimensional. Detailed analysis of the energy transfer in the linear regimeindicates that when the background flow loses its symmetry, it is the mode with the largergrowth rate that is primarily responsible for the extraction of energy from the mean flow.11Two-dimensional numerical simulations are used to examine the nonlinear development of “non-symmetric” Holmboe instabilities. We start by pertiirbing the flow withthe stronger mode predicted by linear theory. We then examine the response of the flowto weaker mode. Finally, we impose both modes. By comparing the development of thesethree flows we are able to study the interaction between the two modes and the effect ofinitial conditions on the development of instabilities. Although the initial developmentof instabilities depends on the initial conditions, this dependence weakens as the densityinterface offset is increased. Also, preliminary results indicate that long term behaviourof the perturbations are independent of initial conditions.Results of the numerical simulations are compared to both symmetric and non-symmetric Holmboe instabilities that have been observed in laboratory experiments.Since we are unable to compute the flow at the high Prandtl number of the salt stratifiedexperimental flows, we run the simulation for a range of increasing Prandtl numbers todetermine how instabilities of thermally stratified flows (with low Prandtl number) andsalt stratified flows differ. Results of these simulations indicate that differences betweenthe experimental and numerical flows can be attributed to the thicker density interfaceand lower Prandtl number used in the numerical simulations. When the density interfaceis thicker and the Prandtl number lower, the waves or billows formed by the instabilitiesare not as sharply defined as in the experiments.111Table of ContentsAbstract iiTable of Contents ivList of Tables viiList of Figures viiiAcknowledgements xi1 Introduction 12 Background 82.1 Governing Equations .. 82.1.1 Small Density Variations 92.1.2 Stream Function Representation 112.1.3 Nondimensional Equations . . . 112.2 Linear Stability Analysis 132.3 Literature Review 163 Linear Stability Analysis 233.1 Derivation of Linear Equations 233.2 Piecewise Linear Profiles 243.2.1 Results: Piecewise Linear Profiles 263.2.2 Three-Dimensional Perturbations: Piecewise Linear Profiles . 28iv4 Nonlinear Numerical Simulations4.1 Numerical Method4.1.1 Numerical Stability4.2 Results4.2.1 Methodology4.2.2 Parameter Values4.2.3 Right Mode4.2.4 Left Mode4.2.5 Both Modes7375777979828489913.3 Smooth Profiles3.3.1 Numerical Method3.3.2 Results: Smooth Profiles .3.3.3 Three-Dimensional Instabilities:3.3.4 EigenfunctionsSmooth Profiles31323442435 Comparison with Laboratory Experiments5.1 Symmetric Case5.2 Non-Symmetric Case6 Summary and Conclusions6.1 Linear Stability Analysis6.1.1 Inviscid Results6.1.2 Viscous Results6.2 Nonlinear Simulations6.3 Comparison with Laboratory Experiments6.4 Future Work128129131145145146146149151151VList of Symbols 153Bibliography 158A Dispersion Relation for Piecewise Linear Profiles 164B Maximum Growth Rates for Holmboe’s Equation 168B.1 Ilolmboe Instabilities 169B.2 Kelvin-Helmholtz Instabilities 170C Energy Equations 172C.1 Kinetic Energy 172C.2 Potential Energy 176D Period of Perturbation Kinetic Energy 178viList of Tables3.1 Values of Prandtl number for water 364.1 Parameters used for the numerical simulations 845.1 Parameters of flows from various laboratory experiments 130viiList of Figures1.1 Schematic of a two layer flow with a thin density interface 61.2 Two-layer flow studied by Holmboe [23] 72.1 Two-layer flow considered by Helmholtz [22] and Kelvin [25] 223.1 Background flow for piecewise linear velocity and density profiles 513.2 Stability diagrams for inviscid case 523.3 Growth rates of most unstable modes for inviscid case 533.4 Phase speeds of most unstable modes for inviscid case 543.5 Linear stability diagrams for small J when = 0 553.6 oc/’J for Kelvin-Helmholtz and Holmboe waves 563.7 c//T for e = 0, 0.25, and 0.5 573.8 Background flow for hyperbolic tangent velocity and density profiles . 583.9 Maximum growth rate in Holmboe regime with varying R 593.10 Linear stability diagrams for smooth profiles: varying R 603.11 Maximum growth rates for smooth profiles: varying R 613.12 Linear stability diagrams for smooth profiles: varying H 623.13 Maximum growth rates for smooth profiles: varying H 633.14 Smooth profiles: varying Re 643.15 Eigenfunctions and correlations when = 0 653.16 Eigenfunctions and correlations when e = 0.5 663.17 Evolution of perturbation kinetic energy and correlations with e = 0. . . 673.18 Evolution of perturbation kinetic energy and correlations with 0.25 . 69viii3.19 Evolution of perturbation kinetic energy and correlations with = 0.5 714.1 Growth rates from linear stability analysis4.2 Absolute stability region of fourth order Runge-Kutta schemePerturbation kinetic energy for right modePerturbation kinetic energy for left mode . .Perturbation kinetic energy for both modes .Effect of grid size when e = 0.5Density and vorticity for right mode only with = 0Density and vorticity for right mode only with e = 0.25Density and vorticity for right mode only with e 0.5Linear growth rates from nonlinear simulationsPositions of waves for right mode onlyAverage phase speeds from nonlinear simulationsEnergy budgets for right mode oniy with e = 0 . . . .Energy budgets for right mode only with e = 0.5Evolution of mean profiles for right mode onlyDensity and vorticity for left mode only with = 0.25 . .Density and vorticity for left mode only with e = 0.5 . .Positions of waves for left mode onlyDensity and vorticity for both modes with 0Density and vorticity for both modes with e = 0.1Density and vorticity for both modes with e = 0.25Density and vorticity for both modes with = 0.54.23 Positions of waves for both modes4.24 Energy budgets for both modes with = 0 1234.34.44.54.64.74.84.94.104.114.124.134.144.154.164.174.184.194.204.214.22100101102103104105106107108109110111112113114115116117118119120121122ix4.25 Wave positions and energy budgets in linear regime for e = 0 . . . 1244.26 Wave positions and energy budgets in nonlinear regime for e 0 . 1254.27 Energy budgets for both modes with E = 0.5 1264.28 Evolution of mean density profiles for both modes . . . 1275.1 Experiments of Zhu [72] 1375.2 Results from numerical simulations for J = 0.3 1385.3 Experiments of Guez & Lawrence [20] 1395.4 Experiments of Lawrence et al. [38] 1405.5 Results from numerical simulations for J = 0.03 with R = 3 . 1415.6 Results from numerical simulations for J = 0.03 with R = 6 . 1425.7 Results from numerical simulations for J = 0.03 with R = 10 . 1435.8 Results from numerical simulations for J = 0.06 . . . 144B.1 Holmboe’s [23] stability diagram 171xAcknowledgementsThis thesis would not have been completed without the constant guidance and encouragement from my supervisor, Dr. Gregory Lawrence. I am indebted to Dr. Brian Wettonfor keeping me numerically honest (at least to order Ax2). In addition, I would like tothank the other members of my supervisory committee for their helpful suggestions inthe preparation of this document.The work in this thesis is based on previous work by Dr. Bill Smyth who was generousin his advice and who promptly responded to my pleas for help. I also greatly benefitedfrom my acquaintances with Dr. Coim Caulfield, Dr. Craig Stevens, and Dr. NoboruYonemitsu.Although I often found dealing with computers and their software frustrating, thefollowing people make the experience less overwhelming: Bob Bajwa, Djun Kim, DaveMoulton, Peter Newbury, Tom Nicol, George Phillips, and Craig Stevens.There are many friends from both the Department of Mathematics and the Environmental Fluid Mechanics group in the Department of Civil Engineering who made mystay at UBC enjoyable and the rain tolerable. I will miss you all.Finally, a special thanks to my husband James for putting up with my strange scheduleas I worked towards completing this thesis and for making sure that I came home to ahot meal and a cold beer at the end of a long day.I received support from the Natural Sciences and Research Council of Canada and aIzaak Walton Killam Memorial Pre-Doctoral Fellowship from the University of BritishColumbia.xiChapter 1IntroductionMixing of stably stratified shear flows is an important aspect of many problems inoceanography, meteorology, and several branches of engineering. Properly quantifyingthis mixing, in order to accurately model turbulence, remains a problem of great importance (see, for example, Gregg [19] and Fernando [14]). Studying mixing in a stablystratified flow has an advantage over other types of transition; stable stratification appears to sufficiently suppress the onset of turbulence to enable distinct events in thetransition to turbulence to be classified (Sherman et al. [50]). Wave-like instabilitiesare the first to occur during transition. Thus understanding the types of wave-like instabilities that can occur, and how they develop, is the first step in gaining completeunderstanding of mixing. The study of such instabilities remains one of the fundamentalproblems in the field of hydrodynamic stability.The traditional method used to determine the stability of a given flow is the normalmode approach to linear stability analysis. This method adds a perturbation (usually,but not necessarily, two-dimensional) to a background flow. The governing equations arethen linearized about this state. The normal mode approach allows one to determinethe wave number, growth rate, and spatial structure of the most unstable modes. Acomplete description of this method is given in section 2.2. The underlying assumptionin linear stability analysis is that, if the perturbation from the mean flow initially consistsof random noise, then it is the most unstable modes that grow fastest and dominate.All investigations into the stability of a flow start with determining its linear stability.1Chapter 1. Introduction 2On its own, however, linear stability analysis does not provide an adequate portrait ofthe evolution of the flow, since it is only valid when the perturbation amplitudes arevery small. It is well known that nonlinear effects rapidly become important and shouldnot be ignored. With today’s computers it is possible to use numerical simulations toexamine these nonlinear effects. One drawback, however, is that only specific cases canbe examined. The results of linear stability analysis can be used as a guide, indicatingwhich flow parameters may produce instabilities of interest. Also, since linear stabilityanalysis tell us which modes are most likely to grow, the nonlinear simulation can begiven a head start by using the results from linear stability analysis to initialize theperturbations.In this thesis, we restrict our attention to the stability of a statically stable stratifiedshear flow whose density interface is much thinner than the shear layer thickness, asdepicted in figure 1.1. Such background flows occur in many physical situations, for example, salt wedge flows (Yonemitsu [70]) and exchange flows (Zhu [72]). Yonemitsu [70]and Zhu [72] have conducted laboratory experiments to model such flows and have observed wave-like instabilities at the density interface of the two fluids. Since the densitygradient acts as a stabilizing force, the instabilities caused by the presence of the velocityshear do not always develop into fully turbulent flow. In this case, studying the onsetand development of instabilities deals directly with many practical issues: for example,determining the interface location and the interfacial friction are of interest in manyengineering problems. In addition, by determining the effect on the mean flow of smallscale structures, a study of the development of interfacial instabilities can lead to a betterunderstanding of large scale dynamics of two layer systems.We assume that the only non-zero component of the background velocity vector is inthe x-direction and that the backgrollnd flow does not vary in this direction. In addition,we allow the background flow to develop in time only, since temporally developing flowsChapter 1. Introduction 3are easier to compute than spatially developing flows. These assumptions result in thebackground flow having parallel streamlines, and hence the term parallel flow is used.A good first approximation to the background flow shown in figure 1.1 is the piecewiselinear inviscid flow studied by flolmboe [23] (figure 1.2). This flow is characterized by asingle parameter, the bulk Richardson number, J 2(p—pi)gh/(p-+-pi)(U_Uwhichis a measure of the relative strength of the stratification to the strength of the velocityshear. The flow depicted in figure 1.2 is linearly unstable to two types of instabilities.When J < 0.07, linear stability analysis predicts a Kelvin-Helmholtz instability. Thisinstability is characterized by zero phase speed and is easily identified by the rollingup of the density interface. Since the original work of Helmholtz [22] and Kelvin [25],much has been accomplished towards understanding this instability, including the firstnonlinear numerical simulations by Tanaka [61] and Patnaik et al. [44]. Recent researchin this area includes the study of secondary instabilities (see, for example, Klaassen &Peltier [31] and Staquet [58]), and understanding how the Kelvin-llelmholtz instabilitydevelops into fully three-dimensional turbulence (Caulfield & Peltier [8]).The second type of instability for the flow depicted in figure 1.2 is a Holmboe instability. This is defined as the superposition of two unstable modes of equal growth rateand equal, but opposite, phase speeds, and is named after Holmboe [23]. A Holmboeinstability is characterized by a right-moving wave whose energy is concentrated abovethe center of the shear layer and a left-moving wave whose energy is concentrated belowthe center of the shear layer. Although Holmboe instabilities occur for J > 0, they arethe only possible instabilities when J > 0.07. These instabilities are of interest as theyare the only mechanism through which the flow may become unstable when stratificationis strong compared to the strength of the velocity shear (J > 0.07).flolmboe [23] studied the case where the flow is symmetric about the center of theshear layer. In nature, the background flow is not necessarily symmetric, e.g. the saltChapter 1. Introduction 4wedge flows described by Yonemitsu [70]. There are two ways that the flow may loseits symmetry: either by having horizontal boundaries placed at different distances fromthe center of the shear layer (Yonemitsu [70]), or by displacing the density interface withrespect to the center of the shear layer (Lawrence et al. [38]). A more general case wouldallow both of these to occur. When the symmetry of the background flow is broken,one of the modes of the flolmboe instability dominates. Non-symmetric backgroundflows have been used to explain why “one-sided” instabilities often occur in laboratoryexperiments, instead of the Holmboe instabilities predicted by the symmetric model (see,for example, Maxworthy & Browand [40], and Koop & Browand [33]). Since relativelylittle research has been done to understand these “non-symmetric” Holmboe instabilities,the main purpose of the present study is to understand how the displacement of thedensity interface affects the development of instabilities in a flow which, for the symmetriccase, supports pure Holmboe instabilities.In chapter 2, we present the background needed for the present study, including anintroduction to the governing equations, a description of linear stability analysis, and areview of the relevant literature. Chapter 3 focuses on how displacing the position of thedensity interface with respect to the center of the shear layer affects the linear stabilityof the flow. This is separated into two distinct parts: inviscid theory for piecewise linearbackground profiles, and viscous theory for smooth background profiles. The inviscidtheory is an extension of the work of Lawrence et al. [38] who modified flolmboe’s originalmodel by allowing a displacement of the density interface. In section 3.2.1 we examine theeffect of imposing horizontal boundaries on this model. We also examine the possibilitythat the perturbations are initially three-dimensional (section 3.2.2).Since inviscid flows with piecewise linear profiles do not take into account such physicalparameters as the Reynolds number and the Prandtl number, section 3.3 is devoted tothe linear stability of viscous flows with smooth background profiles. In particular, weChapter 1. Introduction 5examine how varying the thickness of the density interface affects the stability of the flowwhen the displacement of the density interface is also varied. As with the inviscid case,we examine the effect of horizontal boundaries and the possibility of three-dimensionalinstabilities. We also examine, in detail, the shapes of the normal modes and how theygovern the growth mechanism of the instability.In chapter 4, numerical simulations are used to examine the nonlinear evolution ofnon-symmetric Holmboe instabilities. Since a Holmboe instability consists of two wavespropagating in opposite directions, all previous studies of the nonlinear evolution ofHolmboe instabilities have initialized the perturbations with both modes predicted fromlinear stability analysis. As a starting point in our nonlinear investigations, we initializethe flow with one mode only. We start by perturbing the flow with only the right movinginstability, because linear theory predicts that this mode dominates the flow as the offsetincreases. We then examine the response of the flow to left moving perturbations. Finally,we impose both a left moving and a right moving perturbation. By comparing theresults of these different cases, we are able to determine how initial conditions affect thedevelopment of an instability and the interaction between the two modes.As a final test of the validity of the numerical simulations, chapter 5 concentrates oncomparing the results of the nonlinear simulations with those of laboratory experiments.We start by comparing results of chapter 4 with the recent laboratory observations ofHolmboe instabilities by Zhu [72]. Next we run simulations which match the flow conditions of Lawrence et al. [38] and Guez & Lawrence [20]. These latter flows are clearexamples of “one-sided” instabilities. Comparing results of the numerical simulationswith both symmetric and non-symmetric instabilities observed in the laboratory helpsus identify strengths and weaknesses of the present model.Chapter 1. Introduction 6zxFigure 1.1: Schematic of a two layer flow with a thin density interfaceChapter 1. IntroductionzFigure 1.2: Two-layer flow studied by Holmboe [23]x7U1p1p2Chapter 2Background2.1 Governing EquationsA two-dimensional, unsteady, stratified flow can be described completely in terms of thefollowing variables: u and w (velocity components in the horizontal, x, and vertical, z,directions, respectively), p (density), 0 (concentration of stratifying agent, e.g. heat orsalinity), and p (pressure). Assuming that the fluid is Newtonian, incompressible, andthe equation of state linear, the governing equations are given by1. conservation of momentum:PUt + PUUx + PWUz——p + (u + (u + wz)x), (2.1)PWt + UW PWWz = Pz — pg + (w + (u + w)) (2.2)2. conservation of mass:Pt + PUv + P’Wz + UPx + W,Oz 0, (2.3)3. diffusion equation for stratifying agent:O + uO + wO = ,cA, (2.4)4. the equation of state:p = po[l — 7(0—0w)]. (2.5)8Chapter 2. Background 9In the above set of equations, g, t, i, and y are the gravitational acceleration, themolecular viscosity, the diffusivity of the stratifying agent, and the coefficient of expansionfor the stratifying agent, respectively, and are all assumed to be constant. In addition, Arepresents the two-dimensional Laplacian operator, i.e. A = - +-.The equation ofstate (2.5) and the heat equation (2.4) are often good approximations when the fluid isa pure liquid, where the effects of compressibility are minor (see Pedlosky [45]). In (2.5),p0 is the density at a standard state &o. The validity of the two-dimensional assumptionwill be examined in chapter 3.If p and p are expanded about a state of hydrostatic equilibrium Pa and Pa, as discussed, for example, in Turner [65], the equations of conservation of momentum can bewritten asPUt + PUUx + PWUz =—p + (Au + (u + wz)), (2.6)PWt + PUWx + PWWz = —p — p’g + i (Aw+(u + w)), (2.7)where p = pa(Z) + p(x,z,t), P = Pa(Z) +p(x,z,t), and -a-— —gp.2.1.1 Small Density VariationsFor the flows we consider, the maximum density variation Ap, is small compared tothe standard density p, i.e. Ap/po = 0(10—2), or less. Using this property, the aboveequations can be simplified considerably. For a fluid in which the effects of compressibilityare minor, changes in density have little effect on the conservation of mass, and (2.3) canbe approximated by the continuity equation for an incompressible and nondiffusive flow.u, + w 0. (2.8)In order to derive the above approximation, we use a scaling argument. For the meanbackground flow depicted in figure 1.1, motions in the x-direction are much larger thanChapter 2. Background 10those in the z-direction. We can write u = 0(U) and w 0(W), where W << U. Also,we let x 0(L) and z = 0(D). We expect the time scale in the vertical direction tobe of the same order of magnitude as the time scale in the horizontal direction. Hence0(L/U) = 0(D/W), which implies that D <<L. Finally, since variations in densityare small, the equation of state (2.5) suggests that we may write p = po(i + r(x, z,where p is constant and r = 0(e) < 0(10—2) << 1.Using the equation of state, we can write (2.4) asPt + ‘UPx + WPz = + ‘Pzz0 (PoU) 0 (PoU) 0 (PoW) 0 () 0 () (2.9)We notice that all terms on the left hand side of (2.9) have the same order of magnitude.The dominant term on the right hand side is lcpzz, since D << L. As diffusive terms onthe right hand side cannot dominate, we have 0 (P).Subtracting (2.9) from the conservation of mass (2.3), we obtainpu + P’Wz = Pxx — Pzzo() oQ°°) o() o() (2.10)Now,(2.11)The relations in (2.11) indicate that we can neglect the terms on the right hand side of(2.10) and obtain (2.8).We can now simplify the momentum equations (2.6) and (2.7) using the Boussinesqapproximation. This approximation assumes that variations in density are only importantin the gravitational term of (2.7). This assumption is valid providing the maximumdensity variation Lp is small compared to po, as is the case here (see Turner [65] for amore detailed discussion of the Boussinesq approximation and the types of flows for whichChapter 2. Background 11it is valid). Using the Bonssinesq approximation and the incompressibility condition (2.8)derived above, we can rewrite (2.6) and (2.7) ast + uu + wu2 = — + vAu, (2.12)P0Wt + + = —— yAw, (2.13)where v=,u/po is the kinematic viscosity.Equations (2.12) and (2.13) along with (2.8) and (2.9) form or new set of governingequations. These equations are valid when the fluid studied is a pure liquid for whichdensity variations are small and the effects of compressibility are negligible.2.1.2 Stream Function RepresentationWe can use the continuity equation for incompressible flow (2.8) to define a streamfunction such thatthb th1L’u=— ; w=———. (2.14)Taking 8(2.12)/8z — 8(2.13)/8x and using (2.14) one obtains the following equation forthe stream function.(A + u(A) + w(A) = + (2.15)We have reduced the original system of five equations, (2.1) to (2.5) with five unknowns,‘u, to, p, p, and 6, to a system of two equations (2.15) and (2.9), with two unknowns,‘(x,z,t) and p(x,z,t) = pa(z) + p’(x,z,t).2.1.3 Nondimerisional EquationsIt is convenient to examine the nondimensional equations. This has the advantage ofclassifying various flows by the means of certain nondimensional parameters. A schematicChapter 2. Background 12of the background flow that we will be examining is depicted in figure 1.1. Precisedefinitions of the background profiles are discussed in chapter 3. Defining the meanvelocity as U (U1 + U2)/2 and using the length scale L, the velocity scale 6U, andthe density scale 6p, where L, 6U, and Sp will be defined as needed, the nondimensionalvariables, denoted by ‘, are defined as= 6Uu* + U b = 6ULb* x = Lx*= 6Uw* p = 6pp + po Z = Lz* (2.16)L*— suSubstituting the above definitions into the governing equations (2.15) and (2.9), we obtain(A) + u(A) + w() = Jp + (2.17)andPt + UPx + WPz = RPAP. (2.18)In the above equation the * notation for nondimensional quantities has been dropped. Thethree nondimensional parameters, namely, the bulk Richardson number J, the Reynoldsnumber Re, and the Prandtl number Pr, are defined as6p gL, Re =SU L, Pr =. (2.19)po(6U) vThe Prandtl number, Pr=v/ii, is the ratio of kinematic viscosity to thermal diffusivity, whereas the Schmidt number, Sc=v/i, is the ratio of kinematic viscosity to thecoefficient of mass diffusivity (we are interested in the diffusion of salt). Throughout theremainder of this study, however, we shall refer to the ratio v/it as the Prandtl numberregardless of the stratifying agent (heat or salt).Chapter 2. Background 132.2 Linear Stability AnalysisThe normal mode approach to predicting linear stability is the most common analytictool used to investigate the stability of a flow. For a parallel shear flow, we start with asteady flow in the x-direction which satisfies the governing equations. Perturbations arethen superimposed onto the steady flow and linearized equations for the perturbationsare derived. The normal mode approach assumes that the perturbations are proportionalto exp[ic(x — ct)]. For example, for the governing equations (2.15) and (2.9) we writez, t) = (z) + ‘(x, z, t), (2.20)p(x, Z, = Pa(Z) + p’(x, z, t), (2.21)wherez, t) = {(z)ei_ct)}, (2.22)p’(x,z,t) = {(z)ei_ct)}. (2.23)Here the nondimensional wave number, c, is real and positive, and c = Cr + ic is thecomplex nondimensional wave speed. It is the amplification factor cc which determinesthe linear stability of the problem. If c > 0, the flow is linearly unstable since theperturbations will grow exponentially. Conversely for c < 0 it is stable. If c = 0 then theflow is said to be neutrally stable. Choosing c real and c complex means that the solutionwill grow in time (known as the temporal modes). This is in contrast to the physical casewhere solutions tend to grow in space (known as the spatial modes). Spatial modes canbe found by choosing c real and o complex (see Drazin Reid [12] for a discussion onthe use of spatial modes). We will restrict our study to that of temporal modes, however,since the nonlinear evolution of these is more easily computed. Gaster [17] provides thefollowing approximate transformation between the parameters of the temporal problemChapter 2. Background 14and those of the corresponding spatial problem. If the temporal solution has parametersa = ar and ac = a,cr + iarcj, then the corresponding spatial solution has parametersa = ar — iarcj/cg and = arcr, where Cg = ã(arcr)/ãar is the group velocity. Thisapproximation is valid providing amplification rates are small.When linear stability analysis was first used, many investigators hoped that transitionto turbulence would follow directly from linear instability. It was quickly shown, however,that this is not necessarily the case. Linear analysis is only valid for a very short timebefore nonlinear effects become important. Nevertheless, linear instability does correctlydescribe the onset and early evolution of infinitesimal perturbations. It also seems to givea qualitatively correct indication of the overall stability of the flow (Maslowe [39]). Forthis reason, much of the literature has been devoted to linear stability theory includingbooks by Chandrasekhar [9] and more recently by Drazin & Reid [12], as well as reviewarticles by Maslowe [39] for viscous shear flows, Pellacani [46] for stratified shear flows,and Gage [16] for viscous stratified shear flows.The normal-mode approach of linear stability analysis reduces the system of nonlinear partial differential equations into a system of linear ordinary differential equationswhere the boundary conditions depend on the flow configuration. These new differentialequations describe an eigen-problem where c is the eigenvalue. The coefficients of theresulting differential equations are, in general, not constant. For this reason these equations can only be solved exactly for a handful of special cases. These are mostly limitedto steady flows approximated by either piecewise constant or piecewise linear velocityand density profiles. Although these flows are not physically realistic, they often admitsimple solutions. For more complicated flow profiles it is usually necessary to use eitherasymptotic approximations, or to solve the equations numerically.Having described the linear stability problem, we will now discuss how it is applied.For a given flow configuration, the wave speed c depends on the value of the wave numberChapter 2. Background 15a, which can take on any real value. Depending on the problem, there is, for each value ofa, either a continuous spectrum of wave speeds c, a collntable infinite number of discretevalues of c, a finite number of discrete values of c, or some combination of these cases(see Drazin & Reid [12]). The complete solution for the stream function perturbation b’is of the formz, t) J A(z)ei_t)da (2.24)where c is the flth wave speed associated with the wavenumber a and A are coefficientswhich are determined by the boundary conditions. If there is a continuous spectrum ofwave speeds, then the summation in (2.24) is replaced by an integration. Obviously forthe solution to be stable, we need E{c} < 0 for all a and n whereas for instability weonly need {c} > 0 for one value of a and n.In view of the above form of the general solution, it is easy to see that solving thecomplete linear problem is not a trivial task. Fortullately, one does not usually need todetermine “ completely to gain some understanding of the stability mechanism. In fact,for the most part, the use of linear analysis is restricted to determining the dependenceof the wave speed c on the wavenumber a. This dependence is usually described in theform of a dispersion relation. Once this relation has been obtained, either exactly, withasymptotics, or numerically, one of two problems is usually addressed. The first dealswith transition from laminar to turbulent flow and is concerned with determining at whatvalue of a given parameter, say the Reynolds number, does the flow first become linearlyunstable. This value is often called a critical value. The second problem deals with aflow which is known to be unstable. In this case, given a fixed value of a parameter, forexample, the Reynolds number, we wish to know at what wave number a does the mostunstable mode occur (i.e., the largest value of act). Both approaches are used in thepresent study.Chapter 2. Background 162.3 Literature ReviewIn this section we review the relevant literature. We concentrate on important resultsthat hold for stratified shear flows in general and those that are specific to Holmboeinstabilities.Many of the major results obtained in the study of the stability of stratified shearflows have been for inviscid flows. An advantage of performing linear stability analysison inviscid flows is that exact solutions to the dispersion relation can often be found.The simplest model of a stably stratified shear flow is depicted in figure 2.1. Instabilityin such a flow was first remarked upon by Helmholtz [22] and calculated by Kelvin [25].Details of the calculations are given in section 232 of Lamb [36] and section 4 of Drazin& Reid [12]. Results of linear stability analysis indicate that a mode is unstable ifg(p— p) < kpip2(Ui — U2), where k is the wave number of the normal mode. WhenU1 U2 we can always find a value of k which satisfies this condition, and hence theflow is unstable. If, in addition, P1 = P2, then the flow is unstable at all wave numbers.This type of instability is called a Kelvin-Helmholtz instability and is characterized by a‘rolling up’ of the shear layer into periodic train of symmetric vortices which travel withthe mean flow (Smyth [51]).In 1962, Holmboe [23] examined the stability of a discontinuous density interfacecontained within a piecewise linear velocity shear of finite thickness (figure 1.2). Thestability of such a flow is governed by the bulk Richardson number J. When J < 0.07, thevelocity shear dominates and Kelvin-Helmholtz instabilities are predicted. He discovered,however, that when J > 0.07, the flow is unstable to two modes having equal growthrates and equal, but opposite, phase speeds. Using the method of symmetric waves hewas able to determine the time dependence of the phase speeds for the two modes. Thesuperposition of these two modes is usually referred to as a Holmboe instability, a termChapter 2. Background 17first used by Browand & Wang [4].Smyth [52] offers the following physical explanation for the difference between a Holmboe instability and a Kelvin-Helmholtz instability. If one thinks of the density interfaceas a flexible boundary between the two fluids, then the boundary is almost rigid whenstratification is strong (i.e., large J). In this case, the two modes of the Holmboe instability, i.e., the right moving mode located in the upper layer and the left moving mode in thelower layer, propagate through the flow with very little interaction. When stratificationis weak (small J), however, the interface is very flexible and disturbances in the upperlayer interact with those in the lower layer as they pass through each other. Below somecritical value of J (J < 0.07 for the case studied by Holmboe [23]), the phase speeds of theright and left moving modes become zero and there is a phase locking of the two modes.The two modes then wrap around each other, forming the Kelvin-Helmholtz instability.Hazel [21] developed a numerical method for solving the linear stability problem forcontinuous profiles. He examined the stability of flows where the density interface thickness has a finite value and is characterized by R, the ratio of the shear layer thickness tothe density interface thickness. When R = 1, the only possible instabilities are KelvinHelmholtz instabilities. When R = 5, unstable modes with non-zero phase speed, i.e.,Holmboe instabilities, were found. This indicates that there is some critical value ofR below which Holmboe instabilities do not exist. By numerically computing stabilitydiagrams for various values of R, Smyth [51] found the critical value of R to be approximately 2.4. Smyth [51] and Smyth et al. [53] used numerical simulations to examine thenonlinear evolution of Holmboe instabilities and were able to confirm Holmboe’s prediction of varying phase speed as the right and left modes interact. In addition, Smyth andPeltier have conducted two studies on the transition from Kelvin-Helmholtz to Holmboeinstabilities. The first concentrated on linear theory [54], while the second examined theresults of nonlinear numerical simulations [56].Chapter 2. Background 18Observations of both Kelvin-Helmholtz and Holmboe instabilities have been documented. Perhaps the best known examples of Kelvin-Helmholtz instabilities are thebillow clouds near Denver, Colorado (Colson [11]) and the tilting tube experimentsby Thorpe [63]. Holmboe instabilities are, in general, less well-known than KelvinHelmholtz instabilities. There are numerous examples, however, of laboratory observations of Holmboe instabilities. These illclude the experimental results of Browand &Winant [5], Koop [32] (reproduced in part in Tritton & Davies [64] figure 8.8), and,more recently, Poiiliquen et al. [48], and Zhu [72]. Browand & Wang [4] measured theposition of the stability boundary experimentally and found good agreement with thepredictions of Holmboe [23]. The calculated growth rates for their low Reynolds number flows (Re = (U1 — U2)h/v = 40—300), however, were an order of magnitude lower.There are two possible explanations for this discrepancy. The first lies in the piecewiselinear profiles used by Holmboe [23]. Smyth [51] has computed the stability of an inviscid flow with smooth background profiles. He observed that as the thickness of thedensity interface increases, the maximum growth rate in the Holmboe regime decreases.Secondly, Nishida & Yoshida [43] have computed the stability boundaries for a two-layerstratified shear flow at various Reynolds numbers. They show that the region of instability decreases significantly with decreasing Reynolds number, indicating a decrease inthe maximum growth rates. These results suggest that comparison with a model thatincludes the effects of viscosity and has smooth background profiles might yield betteragreement between experiments and theory.In addition to Kelvin-Helmholtz and Holmboe instabilities, a third type of instabilityhas been observed in stratified shear flows with thin density interfaces. This instabilityoccurs for larger Richardson numbers when Holmboe waves are expected. Instead, theinstability exhibits a “one-sidedness” with the wave protruding into one layer only andis believed to be either, the result of a shift of the density interface with respect to theChapter 2. Background 19center of the shear layer, or, the result of horizontal boundaries being non-symmetricallyplaced with respect to the centers of the shear layer and density interface. This has beenobserved by Keiilegan [26], Yoshida [71], Koop & Browand [33] (see also Maxworthy &Browand [40]), Lawrence et al. [38], and Yonemitsu [70]. In order to better understandthis “one-sideness”, Lawrence et al. [38] modified Holmboe’s model by allowing the density interface to be displaced with respect to the center of the shear layer. They foundthat, as the density interface is displaced, one of the two modes of a Holmboe instability dominates (we will refer to this as a non-symmetric Holmboe instability). For aviscous flow, Yonemitsu [70] examined the effect of non-symmetrically placed horizontalboundaries. Similar results were obtained as discussed below. Smyth [51] made somepreliminary computations on the nonlinear evolution of non-symmetric Holmboe waves.He found that these waves evolve, after many oscillation periods, into weakly nonlinearinterfacial waves which are qualitatively similar to waves that have been observed in thelower atmosphere by Gossard & Richter [18] and Emmanuel et al. [13].Inviscid theory has been used to determine several general criteria for instability. In1880, Rayleigh [49] proved his famous inflection-point theorem: a necessary condition forinstability is that the background velocity profile have an inflection point. Howard [24]proved the semi-circle theorem which states that, for any unstable mode, the complexwave speed, c, must lie inside a semi-circle in the upper half-plane whose diameter isgiven by the maximum amplitude of the background velocity. Miles [41] and Howard [24]showed that a sufficient condition for stability in a stratified parallel shear flow is thatthe gradient Richardson number Ri —(gpa(z)/(paUz(z)2),where U(z) is the horizontalbackground velocity profile, be greater than 1/4 everywhere in the fluid. It is possible,however, to have unstable flows for Ri > 1/4 when the effects of viscosity and thermaldiffusion are included (see Gage [16]). Also, for viscous, thermally dissipative flows,unstable waves whose phase speed lie outside the semi-circle predicted by Howard [24],Chapter 2. Background 20for inviscid flows, have been found (Gage [16]).Since boundaries are always present in the physical world, it is of interest to know howthey affect the stability of a flow. Hazel [21] studied the effects of moving the boundariesin from infinity on an inviscid flow with hyperbolic tangent background velocity and density profiles. Two effects were observed. When the total domain height is greater thanor equal to ten times the thickness of the shear layer, longer wavelengths are destabilizedand are the only wavelengths affected by the presence of the horizontal boundaries. Thisis reasonable as we only expect waves of lengths comparable to the domain height tobe affected. In addition, as the domain height is decreases further, shorter wavelengthsbecome more stable. If the boundaries are moved in closer than some critical distanceapart, the latter effect dominates and the flow becomes stable for all wavenumbers andRichardson numbers. Smyth [51] studied the effect of imposing boundaries on the profilestudied by Holmboe (figure 1.2). Transition from Kelvin-Helmholtz to Holmboe instabilities occurs at a lower value of J when boundaries are present. When the boundaries arenot symmetrically placed, pure Holmboe instabilities are not found. Instead one modeof the Holmboe instability has a larger growth rate. For a viscous flow, Yonemitsu [70]examined how the presence of a lower boundary affects Holmboe instabilities. The wavethat protrudes into the upper layer is little affected by the presence of the boundaryunless the boundary is very close to the shear layer. The wave that protrudes into thelower layer, however, is significantly altered by the presence of the boundary indicatingthat perhaps the two modes are independent.In section 2.2 linear theory was developed under the assumption that perturbationsare two-dimensional. At some point in the transition from laminar to turbulent flowthe perturbations do indeed become three-dimensional, but one must consider whetherthey are initially two-dimensional during the time when linear theory is valid. In 1933,Squire [57] showed that, for an unstratified flow, a three-dimensional disturbance at aChapter 2. Background 21given Reynolds number is equivalent to the problem of a two-dimensional disturbance ata lower Reynolds number. Since a flow becomes more stable with decreasing values ofthe Reynolds number, this implies that all instabilities are initially two-dimensional for ahomogeneous flow. Koppel [35] derived a generalized version of Squire’s result, showingthat a three-dimensional disturbance at given Prandtl, Reynolds and Richardson numbers is equivalent to a two-dimensional disturbance at the same Prandtl number, a lowerReynolds number and a higher Richardson number. Thus any three-dimensional problem can be transformed into a two-dimensional problem. Since, for Holmboe’s instability,there is a region where the flow is destabilized as the Richardson number increases, itmay be possible that the most unstable modes are three-dimensional. This was firstsuggested by Browand & Wang [4]. Two approaches have been used to examine this possibility. Smyth & Peltier [55] solved the three-dimensional problem and showed that, fora class of dissipative, stratified, parallel shear flows, the perturbations evolve directly intothree-dimensional flows without going through the two-dimensional stage. Caulfield [7]calculated maximum growth rates to show that primary three-dimensional instabilitiesare possible in a generalized version of the inviscid three-layer flow studied by Taylor [62].In chapter 3 we use the approach used by Caulfield [7] to examine the possibility of theinitial perturbations being three-dimensional for the inviscid piecewise linear backgroundflow studied by Lawrence et al. [38].Chapter 2. Background 22zU1p1nr2U2Figure 2.1: Two-layer flow considered by Helmholtz [22] and Kelvin [25]Chapter 3Linear Stability Analysis3.1 Derivation of Linear EquationsWe wish to examine the linear stability of the flow governed by (2.17) and (2.18). Asdescribed in Section 2.2, we split the flow field into a parallel mean component and aperturbation field, i.e.,z, t) (z, t) + b’(x, z, t), (3.1)p(x, z, t) = (z, t) + p’(x, z, t), (3.2)where ‘ and p’ are O() <<1. Substituting these expressions into the governing equations(2.17) and (2.18), and collecting terms of equal magnitude, we obtain equations for boththe background flow and the perturbations.1Ut =(3.3)RPPt Pzz,erand(‘) + u(’) + w’u = Jp +(3.4)p + Üp + WZ RePr’where u(x, z, t) = ü(z, t) + u’(x, z, t).Equations (3.3) indicate that the mean velocity and density profiles diffuse over timedue to viscosity and the diffusivity of the stratifying agent. For the purpose of solving23Chapter 3. Linear Stability Analysis 24the linear problem (3.4), however, we will use the initial mean profiles U(z) = ü(z, 0)and Pa(Z) = ö(z, 0). Hence st” and p’ are the perturbations from the initial mean state.Using the method of normal modes, we write‘(x,z,t) = {(z)ci_ct)}, (3.5)z, t) {(z)ei_ct)}. (3.6)Substituting these into the linear equations (3.4) we obtain the eigen-problem:(U — — a2)— Uzz = J——2a +aRe (37—z(U_c)7pa (a2)aRePrNotice that the above set of equations reduce to the Orr-Sommerfeld equation when nostratification is present. Furthermore, the Taylor-Goldstein equation is obtained if thereis no diffusion. Similar equations have been derived by Koppel [35], and Baldwin &Roberts [1].3.2 Piecewise Linear ProfilesThe starting point of our investigation is to examine the stability of an inviscid flowwith two layers of constant density and a layer of constant horizontal shear (figure 3.1).This background flow has the advantage of admitting an exact dispersion relation to thenormal mode analysis. Also, other than the density interface displacement d/h, the onlyadditional parameter is the bulk Richardson number J. Finally, results of inviscid theorywith piecewise linear profiles often compare favourably with experimental results. Forthe reasons given above, an inviscid flow with piecewise linear background profiles is anobvious first approximation to our problem.Chapter 3. Linear Stability Analysis 25We define the background velocity, (u, ) = (U, 0), and density profiles asU1 h/2<z<H/2,U= U1hU2z+UU2 —h/2<z<h/2,(12 —H/2 < z < —h/2, (3.8)1 P1 z>—d,Pa =1P2 z<—d,where h is the shear layer thickness, H is the total depth of the flow, and d is thedisplacement of the density interface from the center of the shear layer (figure 3.1). If wenondimensionalize using L = h/2, SU (U1 — U2)/2, U = (U1 + U2)/2, 5p (p2 — p1)!2,and po = (p2 + P1)!2 then the nondimensional background profiles are1 1<z<H/2,U z —1<z<1,—1 —H/2<z<—1,(3.9)1—iPa =z<—e,where H is now the nondimensional depth of the flow, and e = 2d/h. Since the flow isinviscid, (3.3) indicates that there is no diffusion of the mean flow. Furthermore, (3.7)reduces to the nondimensional Taylor-Goldstein equation:(U— c)( — a2)— + 2J6(z + E) U = 0, (3.10)where 6(z) is the delta function.Assuming the rigid wall boundary condition 0 at z = +H/2 and using thetechniques outlined in Drazin Reid [12] (see appendix A for details of the calculations),Chapter 3. Linear Stability Analysis 26we obtain the following dispersion relation.c4+a132+a= 0, (3.11)where a = a(e, c, J, H) and are given in appendix A. As H —* oc we retrieve the resultsof Lawrence et al. [38]. Also, this is simply a special case of the dispersion relation derivedby Smyth [51] for piecewise linear profiles with non-symmetric horizontal boundaries.3.2.1 Results: Piecewise Linear ProfilesSymmetric Case (s = 0)For e 0 and H oo, we reproduce the results of Holmboe [23] (figure 3.2). WhenJ > 0.07, Rolmboe instabilities are the only instabilities that occur. In the unstableregion, there are two modes of equal growth rates and equal but opposite phase speeds(figures 3.3 and 3.4). The growth rate of the Holmboe instability increases with increasingJ until it reaches a maximum value (figure 3.3) after which it decreases as J increases.When J < 0.07, the flow is unstable to both Kelvin-Helmholtz and Holmboe instabilities.An enlargement of the stability diagram in this region is shown in figure 3.5. When s 0,(3.11) reduces toc4+a2=0 (3.12)Solutions to (3.12) are oscillatory, i.e., Cr 0, when a — 4a < 0. In this case there arefour solutions, two stable and two unstable, of the form c = +Cr + ic. The superpositionof the two unstable modes is Holmboe’s instability. When a— 4a > 0 and a2 > 0,there are two stationary (Cr = 0) unstable solutions to (3.12) with different growth rates.These are Kelvin-Helmholtz instabilities. The contour a — 4a = 0 gives the stabilityboundary for the Holmboe instabilities.Results in figures 3.2 and 3.3 show the growth rates for the dominant Kelvin-HelmholtzChapter 3. Linear Stability Analysis 27instability only. It is clear from figure 3.5 that the Kelvin-Helmholtz instability dominates for J close to 0 whereas the Holmboe instability dominates when J is near 0.07.It is not obvious at which value of the bulk Richardson number, JT, the flow changesfrom having a Kelvin-Helmholtz instability as the most unstable mode to a Holmboeinstability. To determine the value of JT, the maximum growth rate for a fixed value of Jmust be found for both Kelvin-Helmholtz and Holmboe regimes. In appendix B the ci—Jcurves along which the maximum growth rates occur are calculated analytically for eachregion. These results are used to determine the maximum growth rates in each regime(figure 3.6). It is evident that transition from Kelvin-Helmholtz to Holmboe instabilitiesoccurs at JT 0.046.By varying the value of H, we examine the effect of horizontal boundaries on the stability of the flow (figure 3.2). Decreasing H reduces the growth rate of the most amplifiedmode (figure 3.3). Moving the horizontal boundaries closer together, however, also causesthe flow to become unstable at small values of c. This mechanism was initially observedby Hazel [21]. The presence of horizontal boundaries does not significantly affect the stability of the flow, however, until H 10. Since Kelvin-Helmholtz instabilities occur atsmaller wave numbers than Holmboe instabilities, they are stabilized more rapidly as Hdecreases. This results in the transition from Kelvin-Helmholtz to Holmboe instabilitiesbeing shifted to a smaller value of J as the distance between boundaries decreases. WhenH 4, the flow is unstable to Holmboe instabilities only. As the boundaries continue tomove together, the flow becomes stable at small J (see figure 3.2, H = 3). These resultsare consistent with those of Smyth [51].Non-Symmetric Case ( 0)When e > 0 (figures 3.3, and 3.2), the flow has no pure Kelvin-Helmholtz or Holmboeinstabilities. When H = oc, there are two unstable modes, one moving to the right andChapter 3. Linear Stability Analysis 28one moving to the left. As e increases, the growth rate of the right mode increases whileits phase speed decreases. The left mode exhibits the opposite behaviour. Lawrence etal. [38] showed that, for fixed a and J, 4i 1) where (r) and (1) correspond toright and left moving waves, respectively. For a fixed value of J, the range of a for whichthe flow is unstable decreases for the right moving wave and increases for the left movingwave as increases. On the stability diagram, the limb corresponding to the right modestilts towards the J-axis and thus the most unstable right mode for a fixed value of Joccurs at smaller values of a with increasing s. The limb corresponding to the left modetilts towards the a-axis and its most unstable mode occurs at larger a with increasing(figure 3.2). This is the behaviour described by Lawrence et al. [38].We now examine the effect of horizontal boundaries (i.e., a finite value of H) on thestability of non-symmetric flows. To date, this effect has not yet been studied in anydetail (Smyth [51] studied the case corresponding to e = 0.05 and a = 0.2 with H = 00and H = 8). As with the = 0 case, the presence of horizontal boundaries causes theflow to become unstable at small wave numbers. The boundaries have little effect on thefastest growing mode until H drops below 10, at which point growth rates of the fastestgrowing right and left modes decrease with decreasing H. Since right moving instabilitiesoccur at smaller wave numbers than left moving instabilities, they are stabilized morequickly by the presence of boundaries. The result is that growth rates of right movingwaves decrease more rapidly than those of left moving waves. This effect becomes evenmore pronounced as E increases and causes the left mode to dominate when H 5(figure 3.3).3.2.2 Three-Dimensional Perturbations: Piecewise Linear ProfilesUntil now, we have assumed that the perturbations are two-dimensional. Here, we discussthe relationship between two-dimensional and three-dimensional linear instabilities.Chapter 3. Linear Stability Analysis 29The most general form of a three-dimensional disturbance (u’,p’,p’), where u’ =(it’, v’, w’), is that of a traveling wave whose amplitude varies with z and which propagatesat an angle t9 with respect to the x-axis (see White [68]), i.e.,(u’, p’, p’) {(ü(z), (z), (z))ei c0s+ sin-ct) } (3.13)Starting with a background flow of the form u (U(z), 0, 0) with density distributionPa(Z) in hydrostatic eqililibrium with the pressure pa(Z), we superimpose perturbationsof the form (3.13). The equations of motions for an inviscid flow with the Boussinesqapproximation under the linear assumption reduce to(u — c ) — a2) — — Paz1 = 0, (3.14)cost cos2(U_ )coswhich, for the background profile described by (3.9), becomes(u — C ) — 2) — + 2J6(z + e) = (3.15)cos cos29(U_ C )cosThe above is simply the Taylor-Goldstein equation for two-dimensional perturbations,(3.10), with bulk Richardson number J/ cos2 and complex phase speed c/ cos 9. Findingthe growth rate for the three-dimensional problem is equivalent to finding the growth rateof a related two-dimensional problem and multiplying it by cos , i.e.,ac(U, 9, a, J) = cos acj(U, 0, a, J/ cos2 9). (3.16)This is simply a special case of Koppel’s [35] result for the more general case of a stratifiedshear flow with viscosity and dissipation. Koppel’s result is given in a form similar tothe above by Smyth & Peltier [55].The assumption that perturbations are initially two-dimensional is valid providingthe fastest growing two-dimensional mode is more unstable than the fastest growingChapter 3. Linear Stability Analysis 30three-dimensional mode, i.e.,cos9 a*c(U0*J/2t9)< a*ci(U0c*J) (3.17)where is the wave number at which the maximum growth rate occurs for a given valueof the bulk Richardson number. For J > 0, this is equivalent to havingcosi9 a*c.(U 0, , J/ cos2) ac(U, 0, a, J). (3.18)In order to validate the two-dimensional assumption, we must show that a*cj/ decreases with increasing J along the curve of maximum growth rate for a given J. Inappendix B we derive analytic expressions for c—J curves for the case originally studiedby Holmboe (H oc and E 0). These are used to find the maximum growth ratesin each region (figure 3.6). In the Kelvin-Helmholtz regime, the maximum growth ratedecreases with increasing J and thus (3.18) is automatically satisfied. In the Holmboeregime, however, there is a region where the maximum growth rate increases for increasing J. Thus it may be possible that (3.18) is violated in this region. Figure 3.6 shows,however, that a’c//T does indeed decrease with increasing J. This implies that theinitial instabilities are two-dimensional.For the case studied by Lawrence et al. [38] (H = oc, e > 0), we were not able tofind exact expressions for the a*_J curves along which the growth rate is maximized. Wecan, nevertheless, use the results used to produce the stability diagrams (figure 3.2) todetermine these curves for right and left moving modes. These are shown in figure 3.7as are the curves a*cj/../J for e 0, 0.25, and 0.5. For larger values of , the maximumgrowth rate of the right mode decreases with increasing J and thus (3.18) is automaticallysatisfied. For smaller values of e, however, and for the left mode, there are regionswhere the maximum growth rate increases with J. Nevertheless, c*c//J decreases withincreasing J and thus we expect the initial instabilities to be two-dimensional.Chapter 3. Linear Stability Analysis 31Although we have shown that, for the inviscid case, the two-dimensional assumptionis valid, these results should be used with caution. Smyth & Peltier [55] have shownthat it is possible for the most unstable initial instabilities to be three-dimensional whenviscosity and thermal dissipation are included. We will discuss this possibility further insection 3.3.2.3.3 Smooth ProfilesAlthough results for the piecewise linear profile case often agree with experimental results,two parameters are neglected: the Reynolds number and the Prandtl number. Furthermore, initially smooth background profiles are required in the nonlinear numerical simulations used in chapter 4. We now consider background flows with smooth profiles andinclude the effects of viscosity and diffusion. The initial mean velocity, (a, a) = (U, 0),and density profiles are depicted in figure 3.8 and are2z U1+U2(i(z) = tanh—+2 h 2(3.19)pa(z) = _P2_P1tanh(z+d)+P1+P2,where h is the shear layer thickness and is the density interface thickness. Wang [66] andYonemitsu [70] have found that the hyperbolic tangent function is a good approximationto the mean velocity fields measured in their laboratory experiments of salt stratifiedflows with thin density interfaces. Also, for large values of R = h/i’, the above densityprofile is a good representation of thin density interfaces.If we nondimensionalize (3.19) using L = h/2, SU = (U1— U2)/2, U = (U1 + U2)/2,Chapter 3. Linear Stability Analysis 32= (p2 — p’)/2, and po = (p2 + pi)/2, as described in section 2.1.3, then the nondimensional background profiles areU(z) = tanhz,(3.20)Pa(Z) = —tanhR(z+e), Jwhere, as before, is the nondimensional distance of the density interface displacement.With these background profiles, (3.7) become(U — c)( — a2) — Uzz = J—— +a4),aRe (3 21)(U —— Paz = p(zz — a2)where, as defined for the piecewise linear case, J= 2(p2 — pi)gh/(p2+ pi)(Ui — U2) isthe bulk Richardson number. We solve (3.21) numerically using the method describedin the next section.3.3.1 Numerical MethodIn this section we describe the method used to solve (3.21) for a viscous flow with diffusionand initial background profiles given by (3.20). For the inviscid problem with piecewiselinear profiles, it was possible to have an unbounded domain (i.e., —oc < z < oc). Sincewe solve the present problem numerically, we must restrict ourselves to a finite domain.The boundary conditions used are = = 0 at z = +H/2, where H is the totalheight of the domain (figure 3.8). These are equivalent to setting u = w’= p’ = 0 atz = +H/2.Following the method of Klaassen & Peltier [30] we use a truncated Fourier sine seriesto approximate q5 and and seek solutions of the formI \ — N . n7r(z+H/2)—Zn=l fl H‘ ( )— N n7r(z+ /2)pZj ThZZ1PTh S111 HChapter 3. Linear Stability Analysis 33Substituting (3.22) into (3.21) we obtain the following.N N—z{—(U—c)DnqnSn—UzzqnSn—JpnSn} = In=1N N (3.23)z{(U—c)pnSn—pabnSn} = DnpnSn, IaPrRen=1 n=1n7r(z + H/2)where Dn (nK/H)2+aand Sn = sin In order to solve for the coefficientsHg and Pn we use Galerkin’s method. We multiply the above equations by sin mK(z +H/2)/H for rn = 1 .. . N and integrate from z = —H/2 to z = H/2. We obtain 2Nequations that can be written as________JN D ti(m,n)+t2( ,n)Dmcbm+pm = CEbm In1 Dm+aRe iN (3.24){—t3(m, fl)n + ti(m, fl)Pn} + aReprDmPm = CPm Jn=1for m = 1,... , N, and where2 H/2 mTr(z+H/2) nlr(z+H/2)dti(m,n) = — I Usin sin z, (3.25)Hi-H/2 H H2 pH/2 m7r(z + H/2) n7r(z + H/2)t2(m, n)= j sin H sin H dz, (3.26)-H/22 H/2 mr(z + H/2) n7r(z + H/2)t3(m, n) = j I Paz Sill sin —dz (3.27)J-H/2 H HWe have reduced the system of ordinary differential equations (3.21) to a simple matrixeigen-problem (3.24) which is of the formAV = cV, (3.28)where A is a 2Nx2N matrix defined byDmAn,m = (Dnti(m,n)+t2m,n))e(3.29)An,m+N = 8nm (3.30)DmAn+N,m = —t3(m,n), (3.31)jDmAn+N,m+N = ti(m,n) + 6nm, (3.32)aRePrChapter 3. Linear Stability Analysis 346mn is the Kronecker delta function, andV= (3.33),0NIt should be noted that although the above method finds 2N values of the complexphase speed c and the corresponding eigenvalues, it is only the two most unstable modesthat are of interest. These correspond to the values of c with the two largest values ofac.It was determined that N = 100 is usually sufficiently large to resolve eigenvalues andeigenvectors and is the value used, unless otherwise stated.3.3.2 Results: Smooth ProfilesWe are interested in how displacing the density interface affects the stability of the flow.This corresponds to varying the parameter e. Lawrence et al. [38] have studied this effect for the inviscid case originally examined by Holmboe [23] where the only additionalparameter is the bulk Richardson number J. Their results have been reproduced in section 3.2.1 where the effects of boundaries were also studied. For smooth profiles we havethe following additional parameters: the ratio R between the shear layer thickness andthe density interface thickness, the Reynolds number Re, and the Prandtl number Pr.Yonemitsu [70] studied the viscous case without diffusion, but only looked for the mostunstable mode (as opposed to two unstable modes, one moving to the right and anothermoving to the left, which corresponds with the results of Lawrence et al. [38]).Chapter 3. Linear Stability Analysis 35Smyth et al. [53] have shown that an initial mean flow with arbitrary initial valuesof R and Pr will eventually reach a state where R = For this reason, we will setPr = R2, linking the effects of Pr and R on the stability of the flow. Below, we discussthe effect of specific parameters on the linear stability of the flow.Varying RThe most significant difference between the piecewise linear profiles (3.9) examined insection 3.2 and the smooth profiles (3.20) examined here is that we now allow the densityinterface to have a finite thickness, which we change by adjusting the value of the parameter R. Since the expected value of R can vary greatly depending on the stratifying agentand the temperature of the fluid (table 3.1), it is of interest to see how this parameteraffects the stability of the flow. In our study of the effect of R on the stability of theflow with varying , we use the values R 3, 4, 6, and 8 (with Pr = 9, 16, 36, and64, respectively). Although salt stratified flows have much larger values of the Prandtlnumber than examined here, we are unable to numerically resolve such a flow as thereis insufficient diffusion to compute with a feasible value of N (see section 3.3.1). Nevertheless, with the range of Prandtl numbers chosen, we are able to observe trends in howthe stability of the flow changes from thermally stratified to salt stratified flows.Results of Hazel [21] and those of section 3.2.1 indicate that, when H 10, horizontalboundaries do not significantly affect the stability of an inviscid flow. A similar result forthe viscous case is discussed below. Based on the above, we set H 10 to insure thatthe horizontal boundaries do not unduly influence the stability of the flow. We also useRe=300 which Smyth [51] found to be sufficiently large to give good agreement with hisinviscid solutions, facilitating comparison with results of section 3.2.1.For e = 0 we have computed the maximum growth rate in the Holmboe regime forthe above values of R (figure 3.9). As R decreases, the maximum growth rate decreases.Chapter 3. Linear Stability Analysis 36Temperature p Pr°C kg/rn3 kg/rn. s m2/sMolecular diffusion of heat in pure water0 999.843 1.79 x i0 13.4 x 10_a 13.45 999.967 1.52 x 10 13.6 x 10_a 11.210 999.702 1.31 x 10 13.8 x 10_a 9515 999.102 1.14 x 10 14.0 x 10a 8.220 998.206 1.01 x 10 14.2 x 10_a 7.125 997.048 0.89 x 10 14.4 x 10_8 6.230 995.651 0.80 x i0 14.6 x 10_a 55Molecular diffusion of NaC1 (0.01 molar) in water0 1000.322 1.79 x i0 0.1415 x 10_8 12645 1000.436 1.52 x 10 0.1441 x 10_8 105410 1000.162 1.31 x 10 0.1467 x i0 89315 999.554 1.14 x i0 0.1493 x 10_8 76420 998.653 1.01 x i0 0.1519 x iD_a 66625 997.490 0.89 x 10 0.1545 x 10a 57630 996.089 0.80 x i0 0.1571 x 10—8 511Molecular diffusion of salts in sea water0 1028.106 1.89 x i0 0.1360 x 10_8 13525 1027.676 1.61 x i0 0.1386 x 10_a 113010 1026.952 1.39 x 10 0.1409 x 10_a 96115. 1025.973 1.22 x 10 0.1434 x 10_a 82920 1024.763 1.09 x i0 0.1474 x 10 72225 1023.343 0.96 x 10 0.1484 x 10_a 63230 1021.729 0.87 x i0 0.1509 x 10_8 564Table 3.1: Values of Prandtl number for 1) diffusion of heat in pure water, 2) diffusionof NaC1 (0.01 molar) in water, and 3) diffusion of salts in sea water. In all cases p wascalculated using the international equation of state of sea water (see appendix 3 of Pond& Pickard [47]) with salinities of 0, 0.5844, and 35 pss, respectively. 1u was determinedfrom table 16 in Sverdrup et al. [60] using salinities of 0, 0.5844, and 35 ppt, respectively.Values for the coefficient of thermal diffusivity (i for case 1) are from Fischer et al. [15].For cases 2 and 3, c corresponds to the diffusivity of NaC1 in water. Values were takenfrom CRC Handbook of Chemistry and Physics [67] and are for molarity of 0.01 and0.6 for cases 2 and 3, respectively. Where necessary, linear interpolation was used todetermine coefficients that are not listed in referenced tables.Chapter 3. Linear Stability Analysis 37It is well known that Holmboe instabilities do not exist when R = 1 and so there is somecritical valne of R below which Holmboe instabilities do not occur. Using the methodof least squares, we fit a curve of the form (acj)max = a/(b — R) + c, where a, b, and care constants, to the data in figure 3.9. This gives us a critical value of R 2.39. Thisagrees with the results of Smyth [51] who found a critical value of R 2.4 by computingthe stability diagrams for several values of R.We start by examining the results for e = 0 (figure 3.10). We observe the samebehaviour found by Smyth [51] for the inviscid case with smooth profiles. IncreasingR increases the region of instability for both right and left travelling waves: at a givenvalue of J, the range of a for which the flow is unstable increases as R increases; ata given value of a, the range of J for which the flow is unstable also increases withincreasing R. In addition, the growth rates of the fastest growing Holmboe instabilities(figure 3.11) increase with R. Growth rates of Holmboe instabilities when R is finite aresmaller than the piecewise linear case which corresponds to R —* cc. When J is small,stratification is weak and we expect the flow to behave like the homogeneous case forwhich the value of R has no effect. Consequently the stability of the flow for small J doesnot change significantly as R varies. For larger values of J, the wave number at whichthe most unstable mode occurs increases with R. Also, transition from Kelvin-Helmholtzinstabilities to Holmboe instabilities occurs at smaller values of J as R increases.When 6 > 0, the region of instability for both left and right modes increase as Rincreases (figures 3.10). Also, the growth rate of the fastest growing mode increases withincreasing R (figure 3.11). These results are the same as the e = 0 case. The wavenumber at which the most unstable mode occurs, however, changes with both 6 and R.If we fix 6, then the wave number at which the fastest growing mode occurs increaseswith increasing R; conversely, if we fix R and vary e the behaviour depends on the valueof R, as described below.Chapter 3. Linear Stability Analysis 38For a piecewise linear inviscid flow, Lawrence et al. [38] found that as e increases, themost unstable right travelling wave occurs at a smaller wave number, whereas the mostunstable left travelling wave occurs at a larger wave number. On the stability diagrams,this corresponds to the tilting of the limb corresponding to the right moving wave (theright limb) towards the J-axis, and the limb corresponding to the left moving wave (theleft limb) towards the a-axis (figure 3.2). This behaviour remains unchanged for largefinite values of R (figure 3.10, R = 6 and 8). As R decreases below some critical valueof approximately 4, however, this behaviour starts to change (figure 3.10, R = 4). ForR = 4, there is very little tilting of the right limb and the wave number at which the mostunstable right moving wave occurs remains almost constant as increases. When R = 3,however, the wave number of the most unstable right mode increases with increasing e(figure 3.10, R = 3). This change in behaviour with changing density interface thicknesswas documented by Caulfield [7] for three layer inviscid flows. It is seen that for all valuesof R, the left limb tilts towards the a-axis. Also, the left limb lifts off the a-axis asincreases (figure 3.10). This effect becomes more noticeable as R decreases. This did notoccur in the findings of Lawrence et al. [38] since their case corresponds to R —* c. Thewave number of the most unstable left mode increases as e increases, except when R = 3and 0 < < 0.25. In this case, a* decreases as E increases. This is a result of the largeseparation of the left mode from the a-axis. The separation of the left limb from thea-axis affects the range of J for which both left and right moving unstable modes exist.As E decreases, the range of J for which both modes exist decreases.Varying HAbove, it was assumed that H = 10 is sufficiently large for the boundaries to have littleaffect on the stability of the flow. In section 3.2.1 we found this to be true for an inviscidflow with piecewise linear profiles. We now examine the validity of this assumption forChapter 3. Linear Stability Analysis 39smooth profiles. Having studied the effect of R in the previous section, we now set R = 3since it is sufficiently large to allow Holmboe instabilities, yet small enough to convergewithout requiring an excessive number of Fourier modes. Also Pr = R2 = 9 is a physicallyrealistic value of the Prandtl number, corresponding approximately to the diffusion ofheat in water at 10°C (see table 3.1). As above we set Re = 300.When e = 0 (figure 3.12), the effect of H on the stability of the flow differs markedlyfrom the piecewise linear case. For H 10 the boundaries do not significantly affect thestability of the flow except at small values of both a and J. This is not surprising in lightof the results from the piecewise linear case. The fastest growing mode is not affectedby the presence of horizontal boundaries until H < 10, at which point its growth ratedecreases as H decreases (figure 3.13). An interesting difference between the effect of Hon the piecewise linear profiles and smooth profiles is how the transition between KelvinHelmholtz instabilities and Holmboe instabilities changes as H varies. For the piecewiselinear profiles, it was observed that as H decreases the transition occurs at smallervalues of J until only Holmboe modes are present at H 4. This was explained by theKelvin-Helmholtz instabilities being stabilized first as they occur at smaller wave numbersthan the Holmboe instabilities. In contrast to the piecewise linear case where Holmboeinstabilities exist for all J > 0, for smooth profiles there is a range of J between 0 andsome critical value J (which depends on the value of R) for which Holmboe instabilitiesdo not exist. In addition, for smooth profiles, instabilities near the transition betweenKelvin-Helmholtz and Holmboe instabilities occur at smaller wave numbers and have thesmallest growth rates. These are the first modes to be stabilized by the boundaries. Thismanifests itself by the separation of the two regions of instability with a stable regionin between. As H decreases, the region separating the two modes increases while thegrowth rates of the two modes decrease.When e = 0.25 with H = 20 (figure 3.13), the growth rate of fastest growing left modeChapter 3. Linear Stability Analysis 40is much smaller than that of the fastest growing right mode. As H decreases, the growthrates of the two modes decrease. Overall, the growth rate of the right mode decreasesmore rapidly than that of the left mode as H decreases. Nevertheless, the right modealways remains the faster growing of the two modes. This is in contrast to the piecewiselinear case where for small H the left mode dominates. For the piecewise linear casewith e = 0.25 and large H (figure 3.3), the difference between growth rates of left andright modes is not great, and thus the rapid drop in growth rate of the right mode as Hdecreases allows it to fall below the growth rate of the left mode. This is not possible forthe case examined here since the initial difference between the growth rates of right andleft modes is so large.Varying ReIn this section we examine the effect of varying the Reynolds number, Re, on the linearstability of a mean flow described by (3.20). We study the effect of Re by computing thestability diagram for Re = 300 and 100 with R 3. These are shown in figure 3.14.As Re decreases, viscosity becomes more important and we expect an overall stabilization of the modes with decreasing Re. Since viscosity has the greatest effect on smallwave lengths, it is primarily the instabilities at large wave numbers which are stabilizedas Re decreases. This has the effect of decreasing the wave number at which the mostunstable mode occurs. When e = 0, growth rates of the most unstable modes decreasewith decreasing Re and are affected primarily in the transition between Kelvin-Helmholtzand Holmboe instabilities. When E > 0, the fastest growing right mode decreases as Redecreases. There is a slight increase in the growth rate of the most amplified left mode asRe decreases. Although this is surprising, the difference in growth rates is not significant.It is difficult to draw any conclusions about the left mode without further investigation.Chapter 3. Linear Stability Analysis 41Varying e: summaryWe now discuss the overall effect of varying. It was seen that precise behaviour of theflow stability with increasing r depends on the thickness of the density interface (whichis specified by the parameter R). Nevertheless, there are trends that were observed in allthe cases examined. As with the piecewise linear case, there are two types of instabilitieswhen 0. For small J, the flow is unstable to Kelvin-Helmholtz instabilities, whereasfor larger J it is unstable to Holmboe instabilities. When e > 0, there are no pure KelvinHelmholtz or Holmboe instabilities. The Kelvin-Helmholtz instability is replaced by aninstability with positive phase speed. The growth rate of this instability increases with. For larger J, there are still two unstable modes, one moving to the right and anotherto the left. As e increases, however, the growth rate of the right moving wave increaseswhereas that of the left moving wave decreases. This is consistent with the results ofLawrence et al. [38]. On the stability diagram this corresponds to a splitting into twolimbs. The limb corresponding to right travelling waves gets larger. Thus as e increases,the maximum value of J for which the flow is unstable increases. The limb correspondingto left travelling waves shrinks. Hence, the interval of J over which there are two unstablemodes with opposite phase speeds decreases. It was also seen that as E increases, thephase speed of the right mode decreases whereas that of the left mode increases. Themost significant difference between the smooth profile case and the piecewise linear caseexamined in section 3.2 is that when > 0, the difference in growth rates between rightand left modes is much larger for the smooth profile case studied here. We expect thisdifference to be largely due to the different background profiles used in the two studies.This is supported by the similarity between our results for a smooth background flowwith the effects of viscosity and diffusion included and with = 0 and those of Smyth [51]who studied the stability of an inviscid flow with the same smooth background profiles.Chapter 3. Linear Stability Analysis 423.3.3 Three-Dimensional Instabilities: Smooth ProfilesWe wish to examine the possibility of the instabilities initially being three-dimensionalwhen viscosity and diffusion are included. In this case, the relationship between two-dimensional and three-dimensional instabilities for an inviscid flow (3.16) becomescc(U, ,a, Re, J, Pr) = cos ac(U, 0, c, Recos 9, J/ cos2 9, Pr). (3.34)This is Koppel’s result [35] which reduces to Squire’s theorem [57] when J = 0. In orderfor the initial instabilities to be two-dimensional, we requirecos 9. g*cj(U 0, c, Recos 9, J/ cos2 , Pr) < c*cj(U, 0, a, Re, J, Pr), (3.35)where c is the wave number at which the growth rate is maximized for given values of J,Re, and Pr. We can conclude from the above that the initial instabilities may be three-dimensional if the growth rate of the fastest growing mode increases with decreasing Reand/or increasing J.Above, we examined the effect on the stability of the flow of decreasing the Reynoldsnumber. It was observed that this tends to decrease the growth rates of the right modewhen e > 0, and both modes when 0. Therefore the only way that the three-dimensional instabilities can dominate is if the growth rate of the fastest growing modeincreases with increasing J. When e = 0, there is a region where the most amplified growthrate increases with increasing J. This indicates that three-dimensional instabilities maygrow faster than two-dimensional instabilities. In order to verify this, one would haveto solve the three-dimensional linear stability problem to see if the initial perturbationsare indeed three-dimensional. Smyth & Peltier [55] have done this for a range of Rewith R = 3 and 0. They were able to show that, for smaller values of Re, theflow is initially unstable to three-dimensional perturbations. This includes the value ofRe = 300 examined here. When > 0, the growth rate of the fastest growing right modeChapter 3. Linear Stability Analysis 43decreases with increasing J (figure 3.11) implying that the most unstable instabilities aretwo-dimensional.The behaviour of the left mode when e > 0 is different than that of the right mode.First, we observed a slight increase in the growth rate of the most amplified left modewhen the Reynolds number decreased from 300 to 100. Also, there is a region where thegrowth rate of the left mode increases with increasing J (figure 3.11). Thus it may bepossible that the left mode is most unstable to three-dimensional perturbations. Furtherinvestigation is required to determine whether or not this is true. Although the leftmode may be unstable to three-dimensional instabilities, its small growth rate comparedto that of the right mode when e > 0 make it unlikely that we will initially see three-dimensional instabilities. Therefore, for the remainder of this study, we will concentrateon the development of two-dimensional instabilities. It may be possible, however, thatthe eventual growth of the left mode is one mechanism through which the perturbationsbecome three-dimensional at later times.3.3.4 EigenfunctionsSo far, we have only examined the linear growth rates to determine the overall stabilityof the flow. In order to better understand the mechanism of the instability in the linearregime, we now consider the eigenfunctions found from linear stability analysis. Since,for this study, we are interested in flows that have two unstable modes, a right movingwave and a left moving wave, we must examine not only the eigenfunctions of each modeseparately, but also the superposition of the two. We restrict our attention to the studyof the modes found by linear analysis of viscous flows with hyperbolic tangent profilesas these are the starting point of the simulations described in chapter 4. Since linearstability analysis only determines the eigenfunctions up to a constant, we normalize theeigeufunctions so that they have the same amount of perturbation kinetic energy in eachChapter 3. Linear Stability Analysis 44mode.In appendix C we derive the equation governing the evolution of the horizontallyaveraged perturbation kinetic energy, given by (C.14). In the linear regime, we canneglect higher order terms in the perturbations and obtain the following equation for thedevelopment of the perturbation kinetic energy.1 (u12 + w12) =— (w’p’) — + _ (u’An’ + w’’). (3.36)The left hand side of (3.36) is the rate of change of perturbation kinetic energy which isequal to the sum of the following terms (see Bradshaw [3] for homogeneous viscous flowsand Smyth & Peltier [54] for stratified inviscid flows).1. —ai7 represents the extraction of energy from the mean flow by the perturbation,2. (W’p’)z represents the change in perturbation kinetic energy associated with thetransfer of perturbation kinetic energy by the pressure fluctuations,3. —Jw’p’ represents the conversion of perturbation kinetic energy into potential energy,4. u’Au’ + w’Aw’ can be divided into a viscous diffusion term and a viscous dissipationterm.In order to gain a better understanding of the development of the instabilities, we examineterms on the right hand side of (3.36) by calculating the correlations andu’Au’ + w’z.w’. In addition, we examine the eigenfunctions , ii’, I, ‘, and determinedfrom linear stability analysis. Although we do not graph the complex amplitude of thestream function perturbation , it can be inferred from the graph of i’ since the two arerelated by = i’th/o.Chapter 3. Linear Stability Analysis 45Eigenfunctions of Individual ModesWe start by examining the eigenfunctions for the Holmboe instability corresponding toR = 3, Pr = 9, Re = 300, J 0.3, and E = 0. This corresponds to the case studied bySmyth & Peltier [54] for an inviscid fluid and gives us a basis for comparison for the 0case. The eigenfunctions and correlation terms in (3.36) for the right and left modes ofthe Holmboe instabilities are shown in figure 3.15. Since e 0 the growth rates of the twomodes are equal and the phase speeds are equal and opposite. The horizontal velocityperturbation amplitude, ü(z), is concentrated in a narrow region centered near the criticallevel (the height z, such that U(z) = Cr) of the mode. Also, the density perturbationamplitude, ,(z), is concentrated near z 0, the center of the density interface, and isdue to the strong gradient of Pa at this point. The vorticity perturbation amplitude &has two peaks, one just above (below) the density interface and another just above thecritical level (below) for the right (left) mode.For both left and right modes the maximum value of the perturbation kinetic energy,u’2 + w’2, occurs near their critical levels. Since 7 < 0, energy is being extracted fromthe mean flow. This indicates that the modes are unstable. It can be shown (see Klaassen& Peltier [29]) that iii can be expressed in terms of the slope of the phase of z’, i.e.,= (3.37)where ?1’ = r(z) exp(i(z)). A positive slope of indicates that the mode is growing, ifit is negative then it is decaying. For the eigenfunctions shown in figure 3.15 the slope ofthe phase of ti’ is positive and thus the perturbations are unstable. This is in agreementwith the positive growth rate &c = 0.030561 and the correlation The minimum ofoccurs near the critical level and so the perturbations extract energy from the meanflow most efficiently at this level.The correlation represents a vertical flux in the perturbation kinetic energy. IfChapter 3. Linear Stability Analysis 46iF > 0 then the kinetic energy is directed upward. It is directed downward when< 0. For the right mode, 7j < 0 and hence perturbation kinetic energy is beingmoved downward from the critical level towards the centre of the density interface. Forthe left mode, the perturbation kinetic energy is being moved upward from the criticallevel to the density interface.The quantity corresponds to a vertical displacement of mass. When i77> 0 massis moved against the gravitational restoring force, causing an increase ill the potentialenergy of the flow. Potential energy is converted into perturbation kinetic energy whenmass is displaced in the direction of the gravitational restoring force (i < 0). Forboth the right and left modes > 0 indicating that perturbation kinetic energy isbeing converted to potential energy. From the above we can describe the growth of theperturbations as follows. The perturbations extract kinetic energy from the mean flow atthe critical level. The perturbation kinetic energy is then transported towards the centerof the density interface where a portion of it is converted to potential energy.When e = 0.25 (not shown), the general mechanism for growth and the shapes of theeigenvalues are similar to the e = 0 case. The horizontal velocity perturbation has ajet concentrated near the critical level. Also, the density perturbation amplitude, ,3, isconcentrated near the centre of the density interface which now occurs at z —e. Thereare, however, some differences in the relative amplitudes associated with the two modeswhich did not occur for the = 0 case. The magnitude for the density perturbation ofthe left mode is greater than that of the right mode. The vorticity & has two peaks, onenear the density interface and one near the critical level, but instead of being equal inamplitude, as in the = 0 case, the peak near the density interface is greater. Also, theleft mode is starting to develop vorticity near the critical level of the right mode.As with the 0 case, the perturbation kinetic energy, u’2 + w’2, is largest near thecritical levels. The magnitude of iF is larger for the right mode, indicating that it isChapter 3. Linear Stability Analysis 47extracting energy from the mean flow more efficiently and hence will grow more quickly.Finally, and are larger for the right mode than the left. This indicates that theright mode is responsible for a larger portion of the vertical displacement of perturbationkinetic energy. Also, most of the conversion from perturbation kinetic energy to potentialenergy occurs through the right mode.When = 0.5 (see figure 3.16), the differences between the right and left modes areeven more pronounced. The overall shapes of the eigenfunctions for the right mode arenot significantly different from the e 0 case. The eigenfunctions of the left mode,however, have been altered dramatically by the displacement of the density interface.This is similar to results of Yonemitsu [70] who found that for non-symmetrically placedhorizontal boundaries only one mode was affected. it is clear, from 7j and Z7that the right mode is the main mechanism through which energy is exchanged betweenthe various reservoirs (mean kinetic, perturbation kinetic, and potential energies). Forthe left mode, the most efficient extraction of energy from the mean flow no longer occursat the critical level. Instead it has been shifted upward towards the density interface.Superposition of Right and Left ModesTo better understand how right and left modes interact in the linear regime, we superimpose the two modes and examine the evolution of the perturbation kinetic energythroughout one cycle of its period. In appendix D it is shown that the period of thehorizontally averaged kinetic energy is 2K/(c) — 41)). We examine the evolution ofthe terms in (3.36) for one period of the kinetic energy of the superimposed modes inincrements of Ax = r/4 where x = a(4r) — cO)t. For convenience, we suppress thelinear growth components of the perturbations.The evolution of the terms in (3.36) throughout a cycle is shown in figure 3.17 for the= 0 case. In the linear regime, growth and decay of the perturbation kinetic energy isChapter 3. Linear Stability Analysis 48controlled by (which determines the transfer of kinetic energy from the mean flowto the perturbation). The perturbation extracts energy most efficiently at x = 0. Thiscorresponds to a local maximum in the perturbation kinetic energy. ir correspondsto the time when the most amount of perturbation kinetic energy is being returned to themean flow, indicating a local minimum in the perturbation kinetic energy. In betweenthese local extrema there is a period of decay for 0 < < r and a period of growthfor r < x < 2ir after which the cycle repeats itself. Although the perturbations spendequal amounts of time in the growth and decay portions of the cycle, the amplitudes ofare largest when the perturbations are extracting kinetic energy from the mean flow.Thus if growth of the modes were included, we would observe a net gain in perturbationkinetic energy over one cycle. This oscillation in the growth pattern of Holmboe waveswas first predicted by Holmboe [23] and later verified by Smyth et al. [53], and Smyth &Peltier [54]. In chapter 4, we also observe this result numerically, and examine how theexchange of energy changes in the nonlinear regime.Examination of the correlations w’p’, and w’p’, enable us to determine how energy is transferred between the three energy reservoirs: mean kinetic energy, perturbationkinetic energy, and potential energy. When energy is transferred from the mean flow tothe perturbations (i.e., i7 < 0), the critical levels act as sources for the perturbationkinetic energy. Perturbation kinetic energy is then transported vertically away from thecritical level towards the density interface. The density interface then acts as a sink,converting perturbation kinetic energy into potential energy. This is as described abovefrom the results of the eigenfunctions for the individual modes. The opposite behaviouroccurs when the perturbations return kinetic energy to the mean flow. The critical levelsact as sinks and the density interface as a source. Potential energy at the density interface is converted back into perturbation kinetic energy which is transported towards thecritical levels where the perturbations return energy to the mean flow. There is a phaseChapter 3. Linear Stability Analysis 49difference, however, between i7 and At = 37r/4 perturbation kinetic energy isstill being transferred to potential energy even though it is also returning energy to themean flow. Similarly at x = 37r/2, potential energy is being converted to perturbationkinetic energy as is mean kinetic energy.Figure 3.18 shows the evolution of the perturbation kinetic energy and the correlationsfor = 0.25. There is essentially the same oscillation between growth and decay asdescribed for 0. The perturbation kinetic energy is maximum when x = 0 andminimum when x ir with corresponding periods of decay and growth in between theselocal extrema. As with the e = 0 case, the growth and decay of the perturbation kineticenergy is in phase with which determines the exchange of energy between the meanflow and the perturbations. Eigenfunctions shown in figure 3.16 indicate that we can,for the most part, use the line z = —e to separate contributions to left and right modes,particularly if we focus on the behaviour near their respective critical levels. In figure 3.18we observe that the right mode extracts energy from the mean flow more efficiently thanthe left mode. During the time when the perturbations are returning energy to the meanflow, the right mode returns less than the left mode. These two features indicate that theright mode grows more rapidly than the left mode. This is indeed the case, as determinedby the values of oc predicted from linear analysis. Upon examination of it is evidentthat more perturbation kinetic energy is being moved around in the upper layer thanthe lower layer. Thus the right mode is responsible for much of the vertical transport ofperturbation kinetic energy. We believe that the right mode is probably responsible formost of the exchange between perturbation kinetic energy and potential energy, althoughthis is difficult to determine from the evolution ofThe evolution of u’2 + w’2, and the corretations in equation 3.36 for e = 0.5 are shownin figure 3.19. As with the two previous cases, e 0 and 0.25, the perturbation kineticenergy oscillates between a maximum value at x —ir/4 and a minimum at x = 3ir/4.Chapter 3. Linear Stability Analysis 50There is a phase difference, however, between n’w’ and it’2 + w’2, whereas for e = 0and 0.25 they was not. The most efficient transfer of energy from the mean flow to thepertnrbations occur at x = 0. The perturbations return the most energy to the mean flowwhen x Also, the perturbations spend more time throughout the cycle extractingenergy from the mean flow than in the previous cases. If we continue to use z = —eto separate the behaviour associated with the right mode from that of the left mode,we observe that right and left modes no longer exhibit the same behaviour throughoutthe cycle. There is a slight phase shift between the two modes. The right mode startto release energy to the mean flow and resumes its extraction of energy from the meanflow at different times than the left mode. As with the e 0.25 case, the right modeextracts energy more efficiently from the mean flow and returns less than the left mode.Also, the right mode is the primary means through which perturbation kinetic energy istransported vertically towards the center of the shear layer.Chapter 3. Linear Stability Analysis 51Figure 3.1: Background flow for piecewise linear velocity and density profilesChapter 3. Linear Stability AnalysisS00 0.5 1.0 1.5H =4.0__________________________Out. —‘- U0.0 0.5 1.0 1.5_________________________H=30.400.300.200.10__ __ __ _ __ __ __ __ __0.00__________________________________0.0 0.5 1.0 1.5 2.0Figure 3.2: Stability diagrams for inviscid caseLinear stability diagrams for piecewise linear profiles with = 0, 0.25, and 0.5 for H = oo,10, 5, 4, and 3. Contours are of constant growth rate For > 0, solid contourscorrespond to right travelling waves and dashed contours to left travelling waves. When0, growth rates of left and right modes coincide in the Holmboe regime whereasin the Kelvin-Helmholtz regime only dominant modes are shown. Stability boundaries(ac = 0) are the outer most contour in each limb.e=o £O.250.500.400.300.200.10H=ou H=ou H=ou52u.Du0.400.200.10U0.0 0.5 1.0 1.5 2.00.500.400.300.200.10YJ4’Juhf Illwiiii iii — —dill/f !A —-0.0 0.5 1.0 1.5 2.0 C 1 0.5 1.0 1.5 2aH=1O0.0 0.5 1.0 1.5 2.0a0.500.400.300.200.10H=50.500.400.300.200.10H=50.500.40vi,00.500.400.300.200.10)0.00! ;E0.400.300.200.1000 0.5 1.0 1.50.010.0 0.5 1.0 1.5 21a:rQ.OO0.00Il/Il —I .Q,I/A / —/ //‘‘2,;I // ‘2/IIIaj0.5 1.0 1.5 -H=30.400.300.200.101 0— — — —— / // / 1. —‘%- ____.5.c0.400.300.200.10n on.0-1’l////,7•‘7‘ ——I——U=0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0— — — —— __u.(/77/i’l’ I —:—--1.5 2.0Chapter 3. Linear Stability Analysis 530.200.150.100.050.00 0.10 0.20 0.30 0.40 0.508=0.25 8=0.5H=30.250.200.150.10 7 — —/0.050.000.00 0.10 0.20 0.30Figure 3.3: Growth rates of most unstable modes for inviscid caseGrowth rate of most unstable mode at a given J for stability diagrams shown in figure 3.2.Maximum growth rate in the flolmboe regime (e = 0) is indicated by the horizontaldotted line.8=0H=co0.200.15S0.100.05H=oo H=aa0.200.150.10 /0.050.0c0.00 0.10 0.20 0.30 0.40 0JH=1O0.00 0.10 0.20 0.30 0.40 0JH=1O0.200.15a0.100.050/0.200.15a0.100.050.20- 0.15a0.100.050 0.00 0.10 0.20 0.30 0.40 0.50H=1O/u.zD0.200.15a0.100.050.00 0.10 0.20 0.30 0.40 0,H=50.20/0 0.00 0.10 0.20 0.30 0.40 0.50R=50.200.200.150.100.050.00 0.10 0.20 0.30 0.40 0.50H=50.250.200.150.100.050.00 0.10 0.20 0.30 0.40 0.50H=4 H=4 H=40.00 0.10 0.20 0.30 0.40 0.0.200.200.15a0.100.05a0.200.15a0.100.05 /*—0.00 0.10 0.20 0.30 0.40 0.50H=30.200.150.100.050.200.150.10.100.050 0.10 0.20 0.30 0.40 0.0.00 0.10 0.20 0.30 0.40 0.50H=30.250.200.150.10 — — — — — — —0.00 0.10 0.20 0.30 0.40 0.500.00010 0.20 0.30 0.40 0.50a0.40 0.50Chapter 3. Linear Stability Analysis 548=0 8=0.25 8=0.5H=co H=co H=co__ __::__0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J JH=iO H=1O H=1O__0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J JH=5 H=5 H=51.0 . .,—0.5—0.5—ç — ——Is —is. .. ..0.00 0.10 0.20 0.30 0.40 0.50 0. 0 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J SH=4 H=4 H=41.0 1.0—0.5—0.5—0.5—is’—1.00.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J S SH=3 H=3 H=30.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J S JFigure 3.4: Phase speeds of most unstable modes for inviscid casePhase speed of most unstable mode at a given J for stability diagrams shown in figure 3.2.Chapter 3. Linear Stability Analysis 550.0 0.2 0.4 0.6Figure 3.5: Linear stability diagrams for small J when E = 0Growth rates (solid contours) and phase speeds (dotted lines) for the case studied byHolmboe [23]. When J < 0.07 there are two regions of instability: Kelviri-Helmholtzinstabilities (cr = 0) and Holmboe instabilities (the superposition of right and left movingwaves).0.0 0.2 0.4 0.6 0.8a1.0a0.10 -0.08 Dominant Kelvin—Helmholtz Modea0.8 1.0Chapter 3. Linear Stability Analysis 560.10 0.200.200.15-g 0.100.05-0.00 L.0.004-3-7000.00Li0.30 0.40 0.500.10 0.20 0.30J0.40 0.50Figure 3.6: ac/Vi for Kelviu-Helmholtz and Holmboe wavesGrowth rate, ac,j and ac/i/J along a-J curve for which growth rate is maximized forH = oc, e = 0. Results for both Kelvin-Helmholtz (dotted line) and Holmboe (solid line)instabilities are shown.Chapter 3. Linear Stability Analysis 578=0 8=0.25 6=0.5_0.10 0.10 0.100.05 0.05 0.050.00 o.oc0.00 0.10 0.20 0.30 0.40 0.E 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J J2.0 2.G 2X1.5 1.5 1.51.0 l.o 1.0_____0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J JFigure 3.7: ac//J for = 0, 0.25, and 0.5Growth rate, cc and c/\/J along a-J curve for which growth rate is maximized forH = oo, with = 0, 0.25, and 0.5. Both right (solid line) and left (dashed line) modesare shown.Chapter 3. Linear Stability Analysis 58zFigure 3.8: Background flow for hyperbolic tangent velocity and density profilesChapter 3. Linear Stability AnalysisFigure 3.9: Maximum growth rate in Holmboe regime with varying RMaximum computed growth rate in flolmboe regime (+) for R 3, 4, 6, and 8 (Pr=9,16, 36, and 64, respectively) with e = 0, Re=300, and H = 10. The solid line is theapproximation (cj)max = 0.2234/(0.5544 — R) + 0.1218 determined using the methodof least squares.590E-ThC)0.100.080.060.040.020.000 2 4 6 8R10Chapter 3. Linear Stability Analysis0.0 0.5 1.0 1.5 2.0a-38O.25Figure 3.10: Linear stability diagrams for smooth profiles: varying R60Linear stability diagrams for e = 0, 0.25, and 0.5 with Re=300, H = 10, and R 3, 4,6, and 8 (Pr = 9, 16, 36, and 64, respectively). Contours are of constant growth rateact. For e > 0, the solid contours correspond to right travelling waves and the dashedcontours to left travelling waves. When 0, growth rates of left and right modescoincide in the flolmboe regime whereas in the Kelvin-Helmholtz regime only dominantmodes are shown. Stability boundaries (oc, 0) are the outer most contour in eachlimb.a0.0 0.5 1.0 1.5 2.0a0.100.000.0 0.5 1.0 1.5 2.0a1.0a2.00.100.000.0 0.5 1.0 1.5 2.0aR=60.500.400.300.200.10fi nfl0.0 0.5 1.0 1.5 2.0a0.0 0.5 1.0 1.5 2.0aR=8a-30.500.400.300.200.10n finR80.0 0.5 1.0 1.5a0.500.400.300.200.10fi fin2.0 0.0 0.5 1.0 1.5 2.0a0.0 0.5 1.0 1.5 2.0aChapter 3. Linear Stability Analysis 610.150.100.050.000.00 0.10 0.20 0.30 0.40 0.500.15J0.100.050.000.00 0.10 0.20 0.30 0.40 0.50Jn fin0.10 0.20 0.30 0.40 0.500.00 0.10 0.20 0.30 0.40 0.50JFigure 3.11: Maximum growth rates for smooth profiles: varying RMaximum growth rate at a given J for stability diagrams shown in figure 3.10. Maximumgrowth rate in the Holmboe regime (E = 0) is indicated by the horizontal dotted line.=O.25R=3UR30.050.001 “.— :i.. —0.00 0.10 0.20 0.30 0.40 0.50JR=40.001 .-.,— ——0.00 0.10 0.20 0.30 0.40 0.50J:0.050.00 L0.000.150.100.05()00R=6 R=60.000.00 0.10 0.20 0.30 0.40 0.50J0.100.050.000.00 0.10 0.20 0.30 0.40 0.50J0.150.100.05R= 60.15/0.150.100.0500.15:i 0.100.05R=80R=80.150.100.050.00R=80.150.100.05fi fin/0.00 . .. .N0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50J J JChapter 3. Linear Stability Analysisé=O.2562H =200.100.Oc0.0 0.5 1.0 1.5 2.0::0.100.000.0 0.5 1.0 1.5 21E0.200.100.00 ‘.,—‘——-.——_Figure 3.12: Linear stability diagrams for smooth profiles: varying HLinear stability diagrams for R = 3, Pr=9, Re=300, = 0, and 0.25 with H = 20,10, 5, 4, and 3, (N 150, 100, 100, 100, and 100, respectively). For > 0, the solidcontours correspond to right travelling waves and the dashed contours to left travellingwaves. When = 0, growth rates of left and right modes coincide in the Holmboe regimewhereas in the Kelvin-Helmholtz regime only dominant modes are shown. Stabilityboundaries (ac = 0) are the outer most contour in each limb.aH=100.0 0.5 1.0 1.5 2.0aH=42,,0.5 1.0 1.5 20.000.0 0.5 1.0 1.5 21H=30.51,\0.00 .2 —0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0Chapter 3. Linear Stability AnalysisS=O.2563Figure 3.13: Maximum growth rates for smooth profiles: varying HMaximum growth rate at a given J for stability diagrams shown in figure 3.12. Maximumgrowth rate in the flolmboe regime (e 0) is indicated by the horizontal dotted line.H=20 H=20o.010.00 0.10 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 -- 0.40 0.H=1O0.150.100.050.00 0.10 0.20 0.30 0.40H50.150.100.0510 0.10 0.20 0.30 0.40 0.H=40.150.10a0.05n n0.50 0.00 0.10 0.20 0.30 0.40 0.H=50.15a0.050.00 ..-.0.00 0.10 0.20 0.30 0.40 0.H =40.150.100.050.000.00 0.10 0.20 0.30 0.40 0.0H=30.150.100.05non — —. —.10 0.10 0.20 0.30 0.40 0.H=30.0.150.10a0.05n on0.00 0.l0 0.20 0.30 0.40 0.50 0.00 0.10 0.20 0.30 0.40 0.50Chapter 3. Linear Stability Analysis8=0Figure 3.14: Smooth profiles: varying ReLinear stability diagrams for R = 3, Pr=9, H 10, Re=100, and 300 with E = 0, 0.25,and 0.5. Here, N = 50 and 100, for Re = 100 and 300, respectively. Also shown isthe maximum growth rate at a given J. When > 0, solid lines correspond to righttravelling waves and dashed lines to left travelling waves. When E = 0, growth rates ofleft and right modes coincide in the Holmboe regime whereas in the Kelvin-Helmholtzregime only dominant modes are shown. Stability boundaries (ac = 0) are the outermost contour in each limb. Maximum growth rate in the Holmboe regime (E = 0) isindicated by the horizontal dotted line.6480.25 8=0.5-30.500.400.300.200.0 0.5 1.0 1.5 2.0a0.100.000.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5a a0.500.400.30Ro=300 Re=3002.0“I0.200.10(1 Ora0a0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0aR= 1000.150.10a0.050.00 0.10 0.20 0.30 0.40 0.50 )O 0.10 0.20 0.30 0.40Re3000. 500.0.150.100.050.00 0.10 0.20 0.30 0.40Chapter 3. Linear Stability Analysis 65Na) pphase/rrb) wphase/Trc) uphase/7T0.0 0.3amplituded) c— phase/ITe) p— phase/Tr0.00 0.26amplitude5f) (u2+w)/2 g) uw-5__ __0.00 0.02 —0.011 0.011Right mode, ac = 0.030561, cr = 0.598220, and Zc = 0.69037.5f) (u2+w)/2 g) uw h) wp i) wp j) uu+wtw-5 ..__ __ ______0.00 0.02 —0.011 0.011 —0.003 0.003 —0.015 0.015—3x10 3x105Left mode, ac = 0.030561, Cr —0.598220, and z —0.69037.Figure 3.15: Eigenfunctions and correlations when E = 0Eigenfunctions and correlations when e = 0, R = 3, Re = 300, Pr = 9, H = 10, J = 0.3,and a 0.54, a), bth, c) ü, d)’, e), f) (u’2 + w’2)/2, g) iF7, h) i) i)u’L\u’ + w’/w’. z = —e is indicated by the dashed horizontal line and the critical levelsare the horizontal dotted line. For a) through e), both amplitude (thick line) and phase(thin line) are shown. In a) the phase is only given when the amplitude is nonzero.0.0 1.5 0.00 0.15amplitude amplitude0.00 0.65amplitudeNQh) wp i) wp j) uu+ww::.*—0.003 0.003 —0.015 0.015—3x10 3x105a)p b)w c)u d)i e)pphase/IT phase/IT phase/IT phase/IT phase/IT5N—0—50.0 1.5amplitude0.00 0.15amp1 it u d e0.0 0.3 0.00 0.65 0.00 0.26amplitude amplitude amplitudeN—0Chapter 3. Linear Stability Analysis 66i) wp j) uu+ww0.015—3x10 3x105Right mode, cc = 0.110458, cr = 0.296927, and z 0.30615.0.0 1.5 0.00 0.15 0.00 0.26amplitude amplitude amplitudef) (u2+w)/2—0.003 0.003 —0.015 O.015—3x10 3x105Left mode, ac = 0.006855, Cr —0.943582, and z = —1.76975.Figure 3.16: Eigenfunctions and correlations when = 0.5Eigenfunctions and correlations when E = 0.5, R = 3, Re = 300, Pr = 9, H = 10,Na)p b)w c)u d)c e)pphase/it phase/it phase/it phase/it phase/it5—0—50.0 1.5amplitude0.00 0.15amplitude0.0 0.3 0.00 0.65amplitude amplitude0.00 0.26amplitudef) (u2+w)/2 g) uwNOz:0.00 0.02 —0.011h) wp0.011 —0.003 0.003 —0.015Na)p b)w c)u d)c e)pphase/it phase/it phase/it phase/it phase/it0.0 0.3amplitude0.00 0.65amplitude5N—0—50.00 0.02 —0.011 0.011g) uw h) wp i) wp j) uzu+wzwJ = 0.3, and o = 0.67.Chapter 3. Linear Stability AnalysisFigure 3.17: Evolution of perturbation kinetic energy and correlations with = 0.67Evolution of perturbation kinetic energy and correlations when both modes are presentfor e = 0, R = 3, Re = 300, Pr 9, H 10, J 0.3, and c = 0.54.5uw wp wp’0-5____0.000 0.035 —0.024 0.024 —0.012uu+Ww’0.18 4i0bS0.012 —0.1805—0.024 0.024—0.024 0.024—0.024 0.024+—0.012 0.012—0.012 0.012—0.012 0.012—5p0.000 0.0350.000 0.0350.000 0.035—0.18 0.18—0.18 0.18—0.18 0.18+_4i0b4x105—4x1054x105—4x1054x105Chapter 3. Linear Stability Analysis 68—0.024 0.024 —0.012 0.012 —0.18—0.012 0.012 —0.180.18 —4x10 4x1050.180.18 —4x10 4x105uw wp w’p—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x1054x10(u2+w)/20.000 0.0350.000 0.0350.000 0.0350.000 0.035LO(Nr)L—0.024 0.024—0.024 0.024 —0.012 0.012 —0.18Figure 3.17 continued.Chapter 3. Linear Stability Analysis—0.012 0.012 —0.18 0.1869Figure 3.18: Evolution of perturbation kinetic energy and correlations with 0.25Evolution of perturbation kinetic energy and correlations when both modes are presentfor e = 0.25, R = 3, Re 300, Pr 9, H = 10, J = 0.3, and a = 0.52.0.000 0.035 —0.024 0.024uu+ww50ac’1II—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x1054x105—0.024 0.024 —0.012 0.012 —0.18—5110.000 0.035-0.000 0.0350.000 0.0350.18 —4x—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x 4x105Chapter 3. Linear Stability Analysis 70uw wp wp uu+wW—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x10 4x105—0.024 0.024 —0.012 0.012 —0.18 0.18—0.012 0.012 —0.18 0.18—0.012 0.012 —0.18 0.18 —4xFigure 3.18 continued.—4x1054x105—4x1054x105LO-0.000 0.035-0.000 0.0350.000 0.0350.000 0.035—0.024 0.024—0.024 0.024 4x1Chapter 3. Linear Stability AnalysisFigure 3.19: Evolution of perturbation kinetic energy and correlations with e = 0.571Evolution of perturbation kinetic energy and correlations when both modes are presentfor e = 0.5, R = 3, Re = 300, Pr = 9, if 10, J = 0.3, and a = 0.67.uw wp uu+WW500.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18-50.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.18:0.000 0.035 —0.024 0.024 —0.012 0.012 —0.18 0.180IIII><N)II—4x1054x10550—50.000 0.035—4x1054x105—4x1054x105—0.024 0.024 —0.012 0.012I—0.18 0.18Chapter 3. Linear Stability Analysisuru+ww—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x1054x1—0.024 0.024 —0.012 0.012 —0.18 0.18—0.024 0.024 —0.012 0.012—4x1054x105—0.024 0.024 —0.012 0.012 —0.18 0.18 —4x10 4x105Figure 3.19 continued.72u,w wp w ‘p(u2+w)/250—50.000 0.0350.000 0.0350HE0.000 0.0350.000 0.035LI)r)N.—0.18 0.18 —4xChapter 4Nonlinear Numerical SimulationsIn this chapter we wish to study the nonlinear evolution of perturbations in a flow withinitial background profiles given by (3.20). When variations in density are small, (2.17)and (2.18), derived in chapter 2, govern the two dimensional evolution of the flow:(A) + u() + w() = Jp + (4.1)andPt + Px + WPz RePr(4.2)Since we wish to calculate the perturbations from the initial mean state, we writeb(x,z,t) = (z,0) + ‘(x,z,t),p(x,z,t) Pa(Z)+P’(X,Z,t),u(x, z, t) = U(z) + u’(x, z,w(x, z, t) w’(x, z, t),where /“, p’, u’ aild w’ are the perturbations from the initial mean state. Substituting(4.3) into (4.1) and (4.2) we obtain the following equations to be solved.(‘) + (U + u’)(A’) + w’(U + (‘)) = Jp + (AA’ + (4.4)andp + (U + u’)p + w’(pa + P’)z= RePr + pazz).(4.5)It is assumed that the flow is periodic in the direction of the flow (x-direction) with==p’ = 0 at the horizontal boundaries z = +H/2.73Chapter 4. Nonlinear Numerical Simulations 74The initial perturbations are of the form (sb’, p’) = J{(q, ) exp(i&x)}, where d is thereal wave number which determines the period of the flow, and (q, ) are the complexamplitudes determined using the normal modes approach to linear stability analysis. Asdiscussed in chapter 3, the most unstable modes for the right and left moving wavesoccur at different wave numbers when e 0 (see figure 4.1 for the case J = 0.3, R = 3,Re 300, Pr = 9, and H = 10). If these modes were used to initialize the flow, thecomputational domain would have to be very large in order to maintain periodicity inthe x-direction. Since it is believed that it is the linear growth rate that determinesthe relative stability of the modes, we use the following criteria for selecting the wavenumber, &, of the flow. We choose & so that the ratio of the left to right growth ratesat this wave number equals the ratio of the maximum left growth rate to the maximumright growth rate. The values of these are indicated in figure 4.1 for the above mentionedcase. The growth rates at the selected wave numbers are close to the maximum values.Although both modes have the same period in the x-direction, the complex phase speed,c, and complex amplitudes, (, ), are different for left and right moving waves. Theprecise form of the initial perturbations were discussed in section 3.3.4.The initial amplitudes of the perturbations are determined by setting the initialamount of wave kinetic energy in a given mode, K’, sufficiently small (0.0005 for thesimulations examined here) to guarantee that the nonlinear simulation spends some timein the linear regime. The wave kinetic energy K’ is defined asK’(t) — (u/2 + w2). (4.6)Since we are calculating the perturbations from the initial mean flow we write u(x, z, t) =ü(z, t) + u’(x, z, 1) where ü(z, 1) = U(z) + iZ(z, t), and u’(x, z, t) = u’(x, z, t) — F(z, t).Here represents the change in the mean flow due to diffusion. (j= (j dx)/L and() = LH/2 •dz indicate horizontal averaging and vertical integration respectively, whereChapter 4. Nonlinear Numerical Simulations 75L = 2K/& is the length of the domain in the x-direction.Equations (3.3) indicate that the mean velocity and density profiles diffuse over timedue to the presence of viscosity and thermal diffusion. This means that the values of thebulk Richardson number J and the Reynolds number, Re, will increase over time (dueto the increase in the shear layer thickness, Ii). Hence, when the simulation enters thenonlinear regime, the actual values of J and Re will have increased from their initial value.We will continue to describe, however, the flow by the initial values of these parameters.4.1 Numerical MethodIn order to numerically solve (4.4) and (4.5) we use second order centered differencingin space and a fourth order Rnnge-Kutta scheme in time. We let Ax and Az denote thestep sizes in the x- and z-directions respectively. Using the following centered differencingapproximations—()i+,i“-i —Dx 2AxD — ()+i —— 2AzD . (.)j+lj — 2(.) tJJrhL 4 7(Ax)2 ( . )_____D— 2(.) + (‘)i,j—l8z2 (Az)2= +we can write the semi-discretization of (4.4) as= —(Ui + + D’J4(U, + DAhW,) +Chapter 4. Nonlinear Numerical Simulations 76+ += F(’I”, ‘, U) (4.8)and (4.5) as= (Uj + + Dx,j(pazj + Dze,) + RePr’ += G(1’, e’, U, p),j (4.9)where ‘(iAx,jAz,t) and p’(iAx,jLz,t) are the approximate values ofthe stream function and density perturbations at a given grid point, and () indicatesdifferentiation with respect to time. Before discretizing in time, (4.8) must be written inthe form= f(’, e’, (4.10)To this end, we use a method described by Strang [59] which makes use of periodicity ofthe function in the x-direction. The first step is to take the discrete Fourier transformof the semi-discretization (4.8) in the x-direction. Writing .F{.} = , we obtain thefollowing.+ { 2 cos(2kAx) -2 ()2 } + =, e’, U)k, (4.11)Now for each Fourier component, k, (4.11) describes a tridiagonal system. This is easilysolved, giving= )(‘, ‘, U) (4.12)Taking the inverse discrete Fourier transform of the above yields an equation of thedesired form (4.10) which along with (4.9) can be discretized in time using the fourthorder Runge-Kutta scheme.Chapter 4. Nonlinear Numerical Simulations 774.1.1 Numerical StabilityIn this section we discuss the stability of the numerical scheme described above. Canutoet al. [6] give the stability region of the fourth order Runge-Kutta scheme used for thetime discretization in the numerical solution of (4.4) and (4.5). A schematic of theregion is shown in figure 4.2. In order to have stability, we require that the amplificationfactor of the semi-discretized scheme (4.8) and (4.9) lie inside the kidney-shaped region.For simplicity, we will require that the amplification factor, A, lie within the indicatedrectangle, i.e.,o ?{A}(4.13)o < {A} <where At is the size of the time step.To determine an approximate range for the amplification factor A, we first considerthe original equations (4.4) and (4.5). Since these are nonlinear equations in ‘ and p’,we use the method of frozen coefficients to linearize the problem. We know that U 1.Also, since the velocity perturbations u’ and w’ are much smaller than the mean flow wecan use u’ < 1 and w’ < 1 (note that these bounds are much larger than the physicalbounds). Using the above bounds, we can reduce (4.4) and (4.5) to the following linearsystem.(A’) + 2(A’) + U + (A’)2 Jp + (AA’ + U2), (4.14)andp + 2p + Paz + P RePr + pazz). (4.15)Since the quantities U, U, Paz, and Pazz are independent of the grid size, we canrestrict our attention to the stability of the following equations.(A’) = —2(A’) — (A’) + Jp + AA’, (4.16)Chapter 4. Nonlinear Numerical Simulations 78and= —2p + RePA (4.17)Using the notation described by (4.7), the semi-discretization of (4.16) and (4.17) is= — + + (4.18)and= —2D,— + RePr’ (4.19)We use von Neumann analysis to determine the amplification factors of (4.18) and (4.19).Taking the two-dimensional discrete Fourier transform of (4.18) and (4.19) we obtainhk,l= [(Ah)2+ iAhD],l + iJ-k,l(4.20)ek,l= RePr +D gk,lwhere = 27rkx, = 2irlz.z, 15 = —2sin/Ax — sin/Az, and /q, = 2(cos —1)/(Lx)2+ 2(cos — 1)/(L\z)2,the discrete Fourier transform of the Laplacian operator.We can write (4.20) in matrix form,-‘. sine /+zD1hAX I k,l (4.21)0 ek,lRePrThe amplification factors of the semi-discretization are given by the eigenvalues of theabove matrix, ) = Ah/Re + iD, and A2 = Lh/RePr + iD. The restrictions (4.13) givethe following conditions for stability., 2 (cos—1 cosC—1” —4 ( 1 1__\\ —2.79=Re (x)2 + (z)2 ) Re (Ax)2 + )2} 2t 4.22and{A}= —2sin sinCLx Az Ax AzAt’(4.23)Chapter 4. Nonlinear Numerical Simulations 79providing Pr> 1 (and hence {A1} >From (4.22) and (4.23) we can deduce the following restriction on the time step, At,I 2.79Re 1At < mm < >. (4.24)— I_L+J_’o 1 1x tzAlthough (4.24) only guarantees a stable numerical scheme for the linear equations (4.14)and (4.15) with constant coefficients, in general it is still used to determine an appropriatevalue of At for the numerical solution to the nonlinear equations (4.4) and (4.5).4.2 Results4.2.1 MethodologyFor the background flow defined by (3.20), linear theory predicts that the right modedominates when > 0. We start our investigation of the nonlinear evolution by examiningthe development of a. flow when the initia’ perturbation contains only the right mode.Next, we perturb the flow with the mode corresponding to the more slowly growing leftmoving wave. Finally, in order to examine the interaction between the two waves, weimpose both modes.To examine the evolution of the perturbations, we examine the density and vorticitydistributions. At different times throughout the simulation, we plot contours of constantdensity ranging from -0.72 in the upper layer to 0.72 in the lower layer. The contours areof equal increments of 0.18. We superimpose vorticity (w = — w) shading upon thedensity contours, where white corresponds to -1.8 (counter-clockwise vorticity) and blackto 1.8 (clockwise vorticity). These plots enable us to visualize how the perturbationsevolve. They also allow us to qualitatively compare results with those of laboratoryexperiments.It is desirable to compare results of the nonlinear simulation with those of linearChapter 4. Nonlinear Numerical Simulations 80stability analysis. Specifically, we wish to compute the growth rates and phase speedsof the perturbations. In appendix D, we showed that K’(t) oscillates in time when twomodes are present. A lack of oscillations would indicate that the perturbation is composedof one mode only. To compute the growth rate of the simulations, any oscillations in K’(t)are removed using a low-pass filter. The linear growth rate ac can then be determinedfrom the resulting perturbation kinetic energy (see Smyth [51]) K’.ldac = j-1nK’. (4.25)To determine the phase speed cr, we track the position of maximum/minimum values ofdensity at the critical Jevels determined from linear analysis of the right/left modes. Thechange in position with respect to time gives the phase speed. This method of trackingthe right and left moving waves has been used by Smyth et al. [53] in their study ofsymmetric Holmboe waves.In appendix C, the equations governing the evolution of mean kinetic energy, perturbation kinetic energy and potential energy were derived and are given by (C.22), (C.23),and (C.29), respectively. Adopting the notation of Klaassen & Peltier [29], we have= —C(K, K’) — D(K) (4.26)C(k, K’) + C(P, K’) — D(K’) (4.27)= —C(P, K’) — D(P) (4.28)whereC(k,K’) —(u7), D(k) =C(F,K’) -J(), D(K’)=((u) + ()2 + (w)2 + (w)2), (4.29)D(P)= RePrChapter 4. Nonlinear Numerical Simulations 81The notation C(R1,R2) indicates a conversion of energy from reservoir R1 to reservoirR2. D(R1) > 0 represents a loss of energy from reservoir R1 due to diffusion. We shallexamine the evolution of the three energy reservoirs, K, K’, and P over time as well asthe total energy T = K + K’ + P. In addition, we will examine the terms that governthe evolution of energy, i.e., the exchange of energy between reservoirs and the loss ofenergy due to dissipation.Equations (3.3) indicate that the mean velocity and density profiles diffuse over timedue to the presence of viscosity and thermal diffusion. It is of interest to see how thepresence of perturbations in the flow affect the development of the mean flow. To thisend we compare the evolution of the computed mean flow ü(z, t), and ,ö(z, t) (where(z,t) pa(Z) +‘(z,t)) with the exact solutions to (3.3). Solutions to (3.3) with initialconditions given by (3.20) and boundary conditionsTn H — h2( H’—— sec 2) —‘H ‘ — i2l’H —U4— 1) — secn j,j —2’ 2 (4.30)pa(,t) —tau1hR(—*+) Pa_pa(,t)= —tanhR(*+6) = Pa÷can be found using separation of variables (see, for example, Boyce DiPrima [2]). Theyare given byA0 °° I—n27rt1 frz 7r\ü(z,t) = Uz+--+Aexp1ReH2 1c0snff+), (4.31)00 2(z, t)= a + PaZ + B exp { :P:H2 } ( + ), (4.32)where2 Id/2 irz 7rA =H I H/2 {U(z) — uz} cos n + dz, (4.33)f1{Pa(Z)_(a +Apaz)}sinn (+ ) dz, (4.34)Chapter 4. Nonlinear Numerical Simulations 82with = (c-- + pa_)/2 and LPa = (f3a — pa_)/H. Truncating the above series (4.31)and (4.32) at n = 15 is an excellent approximation for t > 100.4.2.2 Parameter ValuesAs already discussed in chapter 3, Smyth [51] found that R> 2.4 is required in order forHolmboe instabilities to exist. Since higher values of R mean a thinner density interface,more grid points are required in the vertical direction in order to numerically resolve theflow. We choose to study the case R = 3 which is sufficiently large to support Holmboeinstabilities yet small enough to allow computations without requiring too fine a grid.This corresponds to Pr = 9 which is the approximate value for the diffusion of heat inwater at 10°C (table 3.1).In the numerical simulations described we set Re = 300 which is sufficiently large toallow the instabilities to fully develop before they are damped by viscosity. Furthermore,it is of the same order of magnitude as that of many laboratory experiments. Smyth etal. [53] have examined in detail the nonlinear evolution of flolmboe waves when e = 0with J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. We choose the same parametervalues in our study of non-symmetric flolmboe instabilities. This allows us to verify ourresults for e = 0 against those of Smyth et al. It also gives us a good basis for comparisonfor the non-symmetric case. J = 0.3 corresponds approximately to the most unstableHolmboe instability (figure 3.11). Also, figure 3.12 indicates that by setting H = 10 weavoid the influence of the horizontal boundaries on the development of the instabilities.We examine the cases e = 0, 0.1, 0.25, and 0.5.To determine the appropriate number of grid points we compute the evolution ofthe perturbed flow at several grid levels and compare the results. As discussed in section 4.1.1, the time steps selected satisfy (4.24). Perturbation kinetic energies are shownChapter 4. Nonlinear Numerical Simulations 83in figures 4.3, 4.4, and 4.5 for the right mode, left mode, and both mode cases, respectively. Since the numerical scheme used is dispersive, it is not surprising that, for e = 0and 0.1, there is a phase shift between the results computed on the nx=nz=64 grid,where nx and nz are the number of grid points in the x- and z-directions, respectively,and those from the nx=nz= 128 grid. Comparing density contours (not shown), we foundthat, although the phase speeds differed between the two grid levels, the flows lookedqualitatively the same. Since there is little phase difference between the 128 by 128 and256 by 256 grids used for e = 0.25 and 0.5, we believe that we can compute on a 128 by128 grid with confidence up to t 400 for e = 0 and 0.1.When = 0.25 and 0.5 we encountered difficulties in numerically resolving the flow.This appeared to be associated with the extreme gradients in the density field and isdemonstrated in figure 4.6. The perturbation in the density gets very steep. Denser fluidfrom the lower layer separates from the wave and is released into the less dense upperlayer. For e 0.5, this results in a decrease in the perturbation kinetic energy (see, forexample, figure 4.5 for E = 0.5) and occurs at different times for the nx=nz=128 andnx=nz=256 cases. Although this type of behaviour has been observed in experimentalstudies, we were not able to resolve the exact time or manner in which it occurs. It wasobserved that this ejection of denser fluid into the upper layer is most likely to occurwhen e is large. This is a result of the rapid growth rate of the right mode and explainswhy we have resolution problems when = 0.5.When e = 0.25, ejection of denser fluid into the upper layer only occurs on thecoarser grid when the initial perturbations only contain one mode. We also notice thatthere is no sudden drop in the perturbation kinetic energy. Since the density contoursremain qualitatively the same, except during the brief period of fluid breaking off fromthe perturbation, and the perturbation kinetic energy compares favourably between thetwo grid levels, we believe that we can compute with confidence on the finer grid upChapter 4. Nonlinear Numerical Simulations 84to t = 400, when oniy one mode is initially present. When the flow is perturbed withboth modes, however, fluid breaks off from the perturbation at both grid levels, but atdifferent times. Although the perturbation kinetic energies do not differ significantly, thedensity contours quickly loose their similarities. Therefore we only have confidence inour results up until t = 170.A summary of the parameters used is given in table 4.1.e c lix liZ At tmaxright mode left mode both modes0 0.54 128 128 0.0625 400 400 4000.1 0.53 128 128 0.0625 400 400 4000.25 0.52 256 256 0.03125 400 400 1700.5 0.67 256 256 0.03125 140 240 140Table 4.1: Parameters used for the numerical simulations when J 0.3, R = 3, Pr 9,Re = 300, and H = 10. The time step size is At. trnax is the largest time for a givensimulation for which the results are reliable.4.2.3 Right ModeWhen the flow is initialized with the right mode only, we expect to see the evolution of aperturbation protruding into the upper layer and moving to the right. This is indeed thecase as seen by the density contours shown in figures 4.7, 4.8, and 4.9 for 230 t 245(e = 0), 190 < t 205 ( 0.25), and 125 t 140 ( 0.5), respectively. This isfurther supported by the lack of oscillations in the initial evolution of the wave kineticenergy (figure 4.3). When e = 0, however, oscillations begin to develop at t 200. Thisindicates the presence of another mode. Figure 4.7, 385 t < 400, shows that it is theleft moving wave that is developing. This is not surprising as left and right modes havethe same linear growth rate for this case and one expects the left wave to eventuallyobtain the same amplitude as the right wave. Another interesting observation is that,Chapter 4. Nonlinear Numerical Simulations 85even though right and left waves have the same linear growth rates, once the left modestarts to grow it is in the nonlinear regime and grows much more slowly than if it hadstarted its growth in the linear regime. As e increases, the right mode dominates for alonger amount of time before the development of the left moving wave. This is due tothe reduced growth rate of the left mode as well as the increased growth rate of the rightmode.An interesting phenomenon is the tilting of the wave in the opposite direction towhich it is travelling (figures 4.7 and 4.8). This tilting is caused by a concentrationof counter-clockwise vorticity that develops at the wave crest and has been observed inlaboratory experiments by Lawrence [37]. The movement of the wave to the right isclosely linked to the concentration of clockwise vorticity that proceeds the perturbationin the density interface and also moves to the right. Ripples in the vorticity to the leftof the wave are a resolution problem caused by the dispersive nature of the numericalscheme and become less noticeable as the grid is refined.We calculate the effective linear growth rate of the perturbations. Since only theright mode is present in the initial development of the perturbations, no oscillationsare observed in the perturbation kinetic energy. Thus we can compute the growth ratedirectly without filtering. In figure 4.10 these are compared to the growth rates predictedby linear stability analysis. When = 0 we have excellent agreement between lineartheory and the simulations. As increases, the growth rate of the simulation is lessthan that of linear theory. The maximum amplitude achieved appears independent ofthe value of (figure 4.3). Consequently, the perturbations spend less time in the linearregime as e increases becailse of the increasing growth rate of the right mode.Holmboe [23] predicted from linear theory, and Smyth et al. [53] showed numerically,that the phase speeds of the two modes in a Holmboe wave are not constant but movefastest when they pass each other and slowest when they are widely separated. WhenChapter 4. Nonlinear Numerical Simulations 86the flow is perturbed with the right mode only, this oscillation in the phase speed isnot initially observed since the left wave is not present to interact with the right wave(figure 4.11). In all cases, the perturbations at both critical levels are initially movingto the right. At later times, however, the perturbation moving at the lower critical levelstarts to exhibit oscillations in its phase speed. This coincides with the development ofthe left moving wave, as indicated by the oscillations in the perturbation kinetic energy(figure 4.3), and occurs at later times as e increases. When = 0, the perturbation atthe lower critical level behaves erratically between 330 < t 370 after which it appearsto be moving to the left. At this point, the phase speed of the right moving wave alsohas oscillations in its phase speed as the two waves pass through each other.In order to compare the initial phase speed with that predicted by linear theory,we compute the average phase speed over a cycle. Differences between predicted andcomputed phase speeds are shown in figure 4.12. In general, the initial phase speedis slightly faster than that predicted by linear theory. Also, the phase speed increasessteadily with time. It appears that the phase speed may be asymptoting to a constantvalue, although further investigations are required to verify this conjecture.By studying the energy budgets given by (4.26), (4.27), and (4.28), we can gain abetter understanding of the mechanism that controls the growth of the instability. Wewill start by discussing the trends that are present for all cases. First, we notice that thetotal energy decreases with time (figures 4.13 and 4.14) and that this is caused by thediffusion of the mean flow. In addition, diffusion of the mean flow causes the potentialenergy P to increase over time. During the initial linear growth of the instabilities, thedecay in mean kinetic energy, K is primarily due to diffusion. As the perturbationsenter the nonlinear regime, however, we see a decrease in the mean diffusion term D(K)and an increase in C(K, K’). At this point, the growth of the perturbation has a moresignificant effect on the total energy budget.Chapter 4. Nonlinear Numerical Simulations 87Since the flow was perturbed with oniy the right mode, we initially see no oscillationsin any of the quantities in the energy budget. During this period, C(K, K’) > 0 andC(P,K) <0 with C(f(,K’) > C(P,K’), indicating a transfer of mean kinetic energyto perturbation kinetic energy and a smaller transfer of perturbation kinetic energy topotential energy. Thus the perturbations are undergoing a net growth as indicated bythe increase in K’. For = 0 oscillations occur at later times in the quantities K, K’, F,C(k,K’), C(P,K’), D(]?) and D(K’). As discussed above, this is due to the eventualgrowth of the left mode and indicates that the two waves are interacting. Notice thatC(K, K’) remains almost always positive and hence the perturbation is constantly extracting energy from the mean flow. Conversely, C(F, K’) oscillates between positive andnegative values. Therefore perturbation kinetic energy is converted to potential energyfor some time and then the process is reversed. The cumulative transfer f C(F, K’)dt isalways negative, however, and so there is a net transfer of perturbation kinetic energy topotential energy. When the left mode becomes sufficiently large, there is a drop in D(K)which is reflected by a flattening of the mean kinetic energy curve. A more detaileddescription of the energy budget when both waves are present will be discussed for the= 0 case initially perturbed with both modes as this permits easy tracking of the twowaves.As e increases, the oscillations occur at later times. For = 0.5 (figure 4.14), thereare no rapid oscillations in the short time for which the simulation was resolved. Thebehaviour is very different from the e = 0 case. We now see oscillations of much longerperiods in C(]?, K’) and C(P, K’). This is reminiscent of the energy budgets observed byKlaassen &i Peltier [29] for the Kelvin-Helmholtz instability. One significant difference isthat, in Klaassen & Peltier [29], C(]?, K’) and C(F, K’) oscillate between positive andnegative values, whereas, in our case, C(K, K’) remains positive and C(P, K’) is negativefor t < 130. Since the Kelvin-Helmholtz instability is characterized by the presence ofChapter 4. Nonlinear Numerical Simulations 88interwinding fingers of heavier and lighter fluid, Klaassell & Peltier [29] attribute thetransfer of potential energy to mean kinetic energy (via the perturbation kinetic energy)to the fingers of heavier fluid being pulled back down and the figures of lighter fluidbeing lifted upward. In the flow examined here, the perturbation is always displacingdenser fluid upward and so there is no transfer of potential energy to mean kinetic energy.When 130 < t < 140, C(K, K’) is negative. During this time, there is a decrease in thewave kinetic energy. The wave amplitude, as seen by the density contours (figure 4.9),also decreases. The transfer quantities C(K, K’) and C(P, K’) are in phase and thusthere is a continuous transfer of energy from mean kinetic energy to perturbation kineticenergy to potential energy. When 0.1 and 0.25 (not shown) the energy budgets showcharacteristics of both the E = 0 and e = 0.5 cases. Eventually rapid oscillations appearas the left wave grows, but there are more slowly varying oscillations superimposed uponthese.One of the main concerns in the study of the evolution of instabilities is to determinehow they influence the development of the mean flow. Mean velocity and density profilesat the final time in the numerical simulation and those given by (4.31) and (4.32) forthe evolution of the mean flow with no perturbations are shown in figure 4.15. Thus,we are able to analyse how the presence of the perturbations affects the development ofthe mean flow. When e = 0, the mixing of momentum (as seen from the mean velocityprofile) is primarily in the upper layer. As e increases, the amount of mixing in theupper layer does not change significantly, but it spreads into the lower layer. Mixing ofmomentum occurs throughout the upper layer, whereas changes in density due to mixingare concentrated near the crest of the wave with a smaller change occurring near thedensity interface. As e increases, the change in density due to mixing increases at boththese locations.Chapter 4. Nonlinear Numerical Simulations 894.2.4 Left ModeExcept for the 0 case, the flow behaviour is markedly different when it is initializedwith the left mode only than when it is initialized with the right mode only. This is adirect consequence of the right mode having the largest linear growth rate when e > 0.As the difference between the growth rates of the left and right modes increases, i.e.,as increases, the right mode influences the flow much more quickly and results in theonset of oscillations in the wave kinetic energy occurring at earlier times as increases(figure 4.4). In fact, for = 0.25 and 0.5, the right mode affects the flow while it is still inthe linear regime. Furthermore, density contours when e = 0 and 0.1 (not shown), showthat only the left mode is present in the early stage of the nonlinear regime, whereas for0.25 and 0.5 the right mode, having a larger growth rate, has had sufficient time togrow and is clearly present (figures 4.16 and 4.17) even at early times in the nonlinearregime. When = 0 the flow behaves exactly the same way as the right mode only case,except that the roles of the right and left moving waves are reversed. Results for > 0indicate that both left and right modes will eventually develop in the flow. We also seethat both left and right moving waves have a concentration of counter-clockwise vorticityin their peaks once their amplitudes become sufficiently large.When > 0 the left mode only case spends longer in the linear regime than theright mode only case due to the smaller growth rates. The linear growth rates fromthe nonlinear simulation are given in figure 4.10 and are only slightly larger than thosepredicted by linear stability analysis. As is seen from the wave kinetic energy (figure 4.4),the = 0.5 case behaves differently from the smaller values of e. It only spends a veryshort amount of time growing at the rate of the left mode. After a short transitionperiod, it continues to grow linearly but with a growth rate of 0.045. As this value liesmidway between the growth rates for right and left modes predicted from linear analysis,Chapter 4. Nonlinear Numerical Simulations 90we conclude that the change in growth rate is due to the presence of the rapidly growingright mode.As with the right mode only case, the phase speed of the left mode does not oscillateinitially (figure 4.18). Also, perturbations at both critical levels move to the left. It isonly when the right wave affects the flow that we observe oscillations in phase speed atthe critical level of the right wave. When e 0.1, the motion of the perturbation atthe critical level corresponding to the right moving wave exhibits interesting behaviour(figure 4.18). Initially the wave moves to the left, as expected. As the wave evolves,oscillations in the phase speed become more dramatic until a turn about occurs at t 200.The wave then travels to the right. Notice that the wave in the lower layer is unaffectedby the turn about and continues to move to the left. Similar behaviour is observed for= 0.25 and 0.5 except that the turn about occurs sooner as increases. Althoughonly the right mode is noticeable in the density contours in figure 4.17 ( = 0.5) the factthat the perturbation at the lower critical level is still moving to the left would suggestthat the left moving mode is growing and would eventually become noticeable in theflow. This is further supported by the results shown in figure 4.16 where both modes areclearly present.Figure 4.12 indicates that the initial average phase speed of the left mode is slightlysmaller than that predicted by linear theory. Also, for e = 0, 0.1, and 0.25, the phasespeed increases steadily over time. This effect becomes less pronounced as increasesand reverses when E = 0.5 for which the phase speed decreases with time. This supportsthe hypothesis that the phase speeds are asymptoting to a single value.When = 0, the energy budget for the nonlinear simulation initially perturbed withthe left mode only is identical to the right only case shown in figure 4.13. This result isnot surprising in view of the symmetry of the flow in this case. As e increases, oscillationsbegin sooner since, with increasing e, the right mode becomes important more quickly.Chapter 4. Nonlinear Numerical Simulations 91Because the difference in linear growth rates between right and left modes increases withincreasing , the interactions between the two waves is weaker for larger values of e. Thisis probably due to the difference in amplitudes of the two waves for the times consideredand is evident in the reduced amplitude of oscillation in C(P, K’). In general, however,the energy budgets look similar to the E = 0 case.When 0, momentum mixing, as indicated by the difference between the evolutionof the mean flow with and without perturbations, is the same as the right only case(figure 4.15) with the amounts of mixing in upper and lower layers reversed. As eincreases, there is a decrease in the amount of mixing in the lower layer and and increasein the amount of mixing in the upper layer. To explain this phenomenon, we examinehow the growth rates of right and left modes change with increasing e. Since the growthrate of the left mode decreases as increases, it grows less as increases which resultsin less mixing in the lower layer. A similar argument holds in the upper layer which isinfluenced by the right moving wave. The right wave causes more mixing with increasinge since it is the dominant mode when e> 0. The change in density due to mixing exhibitsthe same trends observed for the momentum mixing, although, as with the right modeonly case, the magnitude of the change in density is much smaller than that of the changein the mean velocity field.4.2.5 Both ModesFor the symmetric case (e = 0) when the flow is initialized with both modes, growth ofthe wave kinetic energy (figure 4.5) is similar to those of the left and right mode onlycases except that oscillations are present at all times. These oscillations are expectedsince both modes are present from the start and grow at the same rate. The flow lookssimilar to that computed by Smyth et al. [53] with a right moving wave protruding intothe upper layer and a left moving wave protruding into the lower layer. One interestingChapter 4. Nonlinear Numerical Simulations 92difference is that, as seen in figure 4.19, our simulations have not lost the symmetryat later times whereas those of Smyth et al. [53] did. We expect the flow to remainsymmetric even at large times since the background flow is initially symmetric, as arethe horizontal boundary conditions. Results when the flow is initially perturbed withonly one mode further supports the results observed here. It was seen that, when e = 0,the flow tends towards a state where both modes are present.As e increases, oscillations in the perturbation kinetic energy decrease (figure 4.5), especially in the linear regime. This observation indicates that, although the perturbationsare initialized with the same amount of kinetic energy in both left and right modes, theright mode quickly dominates the flow for large values of e. Furthermore, the effectivelinear growth rate lies in between the values predicted by linear theory for the right andleft modes, but is dominated by that of the right mode (figure 4.10). When e = 0 wehave excellent agreement between the predicted linear growth rate and that of the linearsimulation.Even for relatively small values of e, i.e. e = 0.1, linear theory predicts a difference inthe growth rates for left and right modes (figure 4.1). This has a significant effect on thedevelopment of the perturbations. The right mode grows more quickly and initially has alarger amplitude (figure 4.20, 195 t 210). The left mode, however, continues to growand is almost as large as the right mode at the end of the simulation. The numericalsimulation is starting to have resolution problems at later times. In order to compute forlonger times, a finer grid should be used.When 0.25, there is a rapid growth of the right mode which dominates the flow inthe early development of the perturbations (figure 4.21, 100 t 115). Since the flow isinitialized with both modes, the left moving wave is still present but, due to its smallergrowth rate, takes longer to develop. At t = 170 (figure 4.21), we observe both right andleft modes, although the right mode still has the larger amplitude. When e = 0.5, theChapter 4. Nonlinear Numerical Simulations 93flow initially evolves in a similar manner (figure 4.22) with only the right mode present.In fact, for the time in which we have confidence in our results, the flow looks similar tothe right mode only case. This is probably a direct result of the reduced growth rate ofthe left mode at larger values of e. Although the left mode is still present, it has not yethad a chance to grow to an observable level for the results presented here.We discuss the behaviour of the flow for t > 140 when E = 0.5. This is beyond therange of time for which the flow was adequately resolved, but is still instructive from aqualitative point of view. The right mode continues to grow until t 185. At this time,the protrusion into the upper layer separates from the density interface and diffusesinto the ambient fluid (figure 4.6), causing a drop in the perturbation kinetic energy(figure 4.5). After this sudden decrease in magnitude of the perturbation, the growth ofthe right mode appears to be stunted. The left mode, however, continues its slow growthand starts to dominate the flow (figure 4.22). Although we were not able to resolve thisbehaviour numerically, it occurred in several of the simulations. We believe that thenumerical simulation helps us understand the overall mechanism of the instability forlarge values of 6 even though the precise details were not resolved. Furthermore, Koopet al. [34] have observed a similar phenomenon in laboratory experiments. Their flow isinitially unstable to Kelvin-Helmholtz type instabilities. The billows formed mix quicklyinto the ambient fluid. Further downstream, weaker instabilities of a Holmboe natureappear at the interface. This leads us to believe that the phenomena observed in oursimulations are physically possible.The positions of the perturbations over time in both top and bottom layers are shownin figure 4.23. For 6 = 0 our results are in agreement with those observed by Smyth etal. [53] and predicted by Holmboe [23] for the linear regime. During this stage, the wavesspeed up as they pass through each other, and slow down when they are far apart. Asthe perturbations enter the nonlinear regime, however, there is a shift in the behaviour ofChapter 4. Nonlinear Numerical Simulations 94the phase speed. The phase speed still fluctuates during a cycle, but now the waves movefastest as they approach each other, slowing down as they pass through each other. Themethod for tracking the phase speed breaks down as the waves start to interact in thenonlinear regime. Just after the waves pass through each other, there is a jump in theposition of maximum/minimum density (figure 4.23). This jump occurs when densitycontours within the wave have more than one local extremum (figure 4.19, t 245),and indicates an inaccurate measurement of the instantaneous phase speed. The averagephase speed throughout the cycle appears to increase with time (figure 4.12). Thisobservation is consistent with the results of the right mode only and left mode onlycases, but the interaction between the two waves slows down their average phase speed.In light of the above discussion, further investigations are required to better resolve thisbehaviour.When = 0.1, the flow initially behaves the same as the e = 0 case. It has a rightmoving perturbation in the upper layer and a left moving perturbation in the lowerlayer. It also exhibits oscillations in the phase speed as the two waves interact. For40 <t <100, the right moving wave appears to he unaffected by the presence of the leftmoving wave and does not have any oscillations in its phase speed throughout the cycle.On the other hand, the left moving wave has large changes in its phase speed throughoutthe cycle. This is most likely a result of the difference in amplitude of the two modes atthis stage of their development. The right mode has had sufficient time to grow and isclearly present in the density contours, whereas the left mode has not yet grown enoughto be seen (not shown). After t = 100 the two waves behave as in the nonlinear regimefor e = 0, speeding up as they approach each other. For t> 330, the two waves appear tobe moving in a non-smooth way. This behaviour is due to the resolution problem alreadydiscussed above. As for the = 0 case, the average phase speed of the right moving waveincreases over time (figure 4.12) and is initially the same as for the right mode oniy case.Chapter 4. Nonlinear Numerical Simulations 95As the left wave grows, however, its presence starts to affect the phase speed of the rightwave which is slower than the right mode only case. The phase speed of the left waveinitially increases but then decreases as the influence of the right wave strengthens.For larger values of (i.e., E = 0.25 and 0.5), perturbations in the upper and lowerlayers exhibit the same characteristics in the initial linear regime as the two previouscases. Oscillations in the phase speeds as the waves pass through each other are onlypresent for a very short time, however, after which the right moving wave behaves verymuch like the case of the right wave only, with no oscillations in the phase speed. Theperturbation at the left critical layer looses it distinctive motion to the left. It bouncesbetween moving to the left and moving to the right in no particular fashion, indicatillga complicated interaction between the slowly growing left moving wave and the rapidlygrowing right moving wave. For t > 140 in the e = 0.25 case, the left moving wave hasgrown large enough for the two waves to interact in the nonlinear regime as describedabove. The reason it takes longer to reach this state than the e = 0.1 case is the leftmode has a much slower growth rate and hence takes longer to achieve a large enoughamplitude to affect the motion of the right mode. For e = 0.5, the simulation does notlast long enough for the left mode to grow sufficiently large to affect the evolution ofthe right mode, but the above results suggest that the flow would eventually exhibit thesame behaviour as the other cases with smaller values of e.For both r 0.25 and 0.5, the average phase speed of the right mode increases steadilyover time (figure 4.12) and appears unaffected by the left moving wave. The phase speedof the left mode, when = 0.25 appears to be slowed down by the presence of the rightwave when compared to the left only results. The simulation for E — 0.5 is too short todraw any conclusions about how the phase speed of the left mode varies with time.The energy budget for e = 0 is shown in figure 4.24 and was discussed in detail bySmyth et al. [53]. We will review the results here and discuss how they tie into theChapter 4. Nonlinear Numerical Simulations 96previous cases. As before, we see that the total energy T decreases with time due todiffusion. Although there are rapid oscillations in K, K’, and F, there are no oscillationsin the total energy indicating a balance in the exchange of energy between the reservoirs.The long term behaviour of K, K’, and P are as in the right mode and left mode onlycases. Since j C(k, K’)dt > 0 for any value of t, there is a net transfer of energyfrom the mean flow to the perturbations. The greatest transfer of energy, however, isbetween the potential energy and the perturbation kinetic energy. Although C(P, K’)oscillates between positive and negative quantities, we have that J C(P, K’)dt < 0 atany instant. Thus there is a net conversion of perturbation kinetic energy to potentialenergy. Nevertheless, f C(K, K’)dt >— J’ C(P, K’)dt indicating that a portion of theenergy transferred to the perturbation kinetic energy from the mean flow remains in theperturbation reservoir. Diffusion of the mean kinetic energy decreases as K decreases.Also, since both modes are present at all times, there is no sudden drop in D(K), aswhen the flow is initially perturbed with one mode only, but a steady decrease overtime. Initially D(K) > G(K, K’) but as the flow enters the nonlinear regime, the twoquantities have the same order of magnitude and both contribute to the overall decreasein K. Diffusion of the perturbation kinetic energy increases as K’ increases whereasD(P) remains constant throughout the simulation.From figure 4.23 we can determine when the right and left wave pass through eachother and when they are maximally separated (for example, figure 4.19, at t 235and t 240, respectively). Studying the energy budgets we can determine how theinteractions of the two waves affect the growth of the perturbations. We will examine thisinteraction during the linear stage and in the fully nonlinear regime (see figures 4.25 and4.26, respectively, for enlargements of positions and energy budgets during these times).While the waves are growing linearly, the growth and decay of the perturbation kineticenergy is controlled by the transfer of mean kinetic energy to the perturbation kineticChapter 4. Nonlinear Numerical Simulations 97energy, and vice versa. Extraction of energy from K is optimal just after the waveshave passed each other. This corresponds to a maximum in the perturbation kineticenergy. As the waves approach each other, there is a brief period when C(K, K’) < 0,corresponding to a minimum in K’. In general, the transfer of energy is from meankinetic energy to perturbation kinetic energy and then from perturbation kinetic energyto potential energy (or in the reverse direction). There is, however, a slight phase shiftbetween C(K, K’) and C(P, K’). This was observed when studying the evolution of theeigenfunctions in chapter 3.Once the waves enter the nonlinear regime (figure 4.26), transfer of energy from Kto K’ has two maxima. One occurs just before the waves pass through each other andthe other just after. During the earlier stages of the nonlinear regime the most efficientextraction of energy from the mean flow occurs just after the waves have passed througheach other. At later times, however, it occurs as the two waves are approaching each other.When the waves are maximally separated, there is a brief period when C(k, K’) < 0and thus the perturbations are giving up energy to the mean flow. As left and rightmoving waves approach each other, C(P, K’) > 0. This corresponds to a growth periodfor K’ and a decay period for P. During the part of the cycle when the waves are movingaway from each other, there is a reversal in the behaviour of C(P, K’) corresponding toa decrease in the perturbation kinetic energy and an increase in the potential energy. Incontrast to the linear regime where oscillations in K’ are governed by C(K, K’), in thenonlinear regime they are linked to C(P, K’). Also, in the linear regime, maxima in K’occur just after the waves have passed through each other, minima occur as the wavesare approaching each other. In the nonlinear regime, however, maxima in K’ occur asthe waves pass through each other whereas minima occur when the waves are maximallyseparated.When E = 0.1 and 0.25 (not shown), the energy budgets look similar to those of theChapter 4. Nonlinear Numerical Simulations 98e 0 case. Therefore, the two waves need not have equal amplitudes for the nonlinear interaction to exhibit the same behaviour described above. If e is sufficiently large,however, the energy budget can look very different even when the flow is initially perturbed with both right and left modes, as seen by examining the = 0.5 case. Theenergy budget when 0.5 (figure 4.27) looks similar to that of the right mode only(figure 4.14). This similarity is not surprising as the density contours (figures 4.9 and4.22) are virtually identical for the two cases. One difference between the two energybudgets is that when the flow is initialized with both modes, there are small oscillations inthe transfers between the reservoirs, C(]?, K’) and C(P, K’). These oscillations indicatethat the left mode, however small it may be, is interacting with the right mode. Initially,these interactions are very weak and we see a continuous transfer of energy from K to K’with a smaller transfer from K’ to P. With increasing time, C(P, K’) oscillates betweenpositive and negative values. This behaviour in the energy budgets looks similar to thatof the right mode only when = 0.1 and 0.25, suggesting that, given enough time, theleft mode will grow and the energy budget will behave as those described above whenthe amplitudes of the two waves are similar.Mean velocity and density profiles are shown in figure 4.28. When E = 0, the changesdue to mixing in the mean velocity and density profiles are the same in both upper andlower layers. This is as expected since the perturbations remain symmetric even near theend of the simulation (figure 4.19). As observed earlier, momentum mixing tends to occurthroughout the entire domain whereas changes in density due to mixing are concentratednear the waves’ crests. When the perturbations initially contained the right or left modeonly, there was a small change of density due to mixing near the center of the densityinterface. Changes were in opposite directions for right and left modes. When bothmodes are present, there is no change in density due to mixing at the density interface.Hence, any changes mean density at the interface are due to diffusion and changesChapter 4. Nonlinear Numerical Simulations 99in density due to mixing for left and right modes cancel at this level. As E increases,the amount of momentum mixing increases in the upper layer and decreases in the lowerlayer. This result is consistent with results of the right and left mode only cases. For theprofiles shown it is difficult to determine the effect of varying E on the density changesdue to mixing. Based on previous results, we expect the same type of behaviour as forthe momentum. It also appears that as e increases, there is a larger change in densitydue to mixing at the interface. When is larger, the left mode is not sufficiently strong toproduce an opposing change in density at the interface and thus there is no cancellationof the change in density due to the right mode at this level.Chapter 4. Nonlinear Numerical Simulations 1000.1Figure 4.1: Growth rates from linear stability analysisGrowth rates from linear stability analysis for J = 0.3 with Re = 300, Pr = 9, R 3,and H 10. The solid line corresponds to the right moving wave and the dotted line tothe left moving wave. When E = 0 the right and left growth rates coincide. The verticaldashed line indicates the value of ö 0.54, 0.53, 0.52, and 0.67 for E = 0, 0.1, 0.25, and0.5, respectively, the wave number used in the simulations.000.120,000.120.000.0 1.5a0.25000.120.000.120.000.0 1.5a0.50.0 1.5a0.0 1.5aChapter 4. Nonlinear Numerical Simulations 1012 -planeFigure 4.2: Absolute stability region of fourth order Runge-Kutta schemeSchematic of absolute stability region of fourth order Runge-Kutta scheme. The schemeis stable inside the kidney shaped region. See Canuto et al. [6] for the precise shape ofthe stability boundary.2f2-2.79Chapter 4. Nonlinear Numerical Simulations 102—4—68=0::nxnz=128, dtO.O625nx=nz=64, dt=O.125 —0 200 300 4C8=0.1U-2-nx=nz=128, dt=O.0625nx=nz=64, dt=O.1250 100 200 300 408=0.250-2-nx=nz=256, dt=O.03125nx=nz=128, dt=O.0625I—100 100 200 300 408=0.5—8—4—6—8U00 100 200 300 400Figure 4.3: Perturbation kinetic energy for right modePerturbation kinetic energy when flow is initialized with the right mode oniy with varyinggrid size for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10.—2—4—6—8Chapter 4. Nonlinear Numerical Simulations 103o ioo 200 300 400=O.2500 100 200 300 4000 100 200 300 400Figure 4.4: Perturbation kinetic energy for left modePerturbation kinetic energy when flow is initialized with the left mode only with varyinggridsizeforJ=0.3,R=3,Pr=9,Re=300, andH=10.8=0.1Chapter 4. Nonlinear Numerical Simulations 104C,i-6nxnz=128, dt=O.0625:: nx=nz=64, dt=o.1250 100 200 300 408=0.1C)::—6——: nx=nz=128, dt==O.0625-8 nx=nz=64, dt=O.125—Ic’0 100 200 300 4C=0250-2- ,,..nxnz=256, dt=O.03125—nx=nz=128, dt=O.0625ci—100 100 200 300 4C—4—6—8008=0.50 100 200 300Figure 4.5: Perturbation kinetic energy for both modesPerturbation kinetic energy when flow is initialized with both modes with varying gridsizeforJ=0.3, R=3, Pr=9, Re=300, andH=10.—2—4—6—8—In400Chapter 4. Nonlinear Numerical Simulations—2—4105Density contours for both modes for e = 0.5 with J = 0.3, R = 3, Pr = 9, Re = 300,and H = 10. x has been scaled with respect to its period 2ir/. Top and bottom rowsare results for nx=nz= 128 and nx’nz=256, respectively. Breaking off of fluid parcel wasnot resolved numerically.420.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0Figure 4.6: Effect of grid size when e = 0.5Chapter 4. Nonlinear Numerical Simulations 106Figure 4.7: Density and vorticity for right mode oniy with = 0Results of nonlinear simulations when the flow is initialized with the right mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, if = 10, and e = 0. Contours are of constant densityp = p(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respectto its period 2ir/&.0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0Chapter 4. Nonlinear Numerical Simulations 107Figure 48: Density and vorticity for right mode only with e = 0.25Results of nonlinear simulations when the flow is initialized with the right mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and a = 0.25. Contours are of constantdensity p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled withrespect to its period 2K/&.0.0 0.2 0.4 0.6 0.6 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4Chapter 4. Nonlinear Numerical Simulations 108Results of nonlinear simulations when the flow is initialized with the right mode oniy forJ = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and = 0.5. Contours are of constantdensity p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled withrespect to its period 2ir/.Figure 4.9: Density and vorticity for right mode only with e = 0.5Chapter 4. Nonlinear Numerical Simulations 1090.15 I I I I I I+0.10U+0.050.00 I I I I I I0.0 0.2 0.4 0.6Figure 4.10: Linear growth rates from nonlinear simulationsLinear growth rates from: 1. linear stability analysis (A), 2. nonlinear simulations withleft or right mode only (+), 3. nonlinear simulations with both modes () when J 0.3,R 3, Pr = 9, Re = 300, and H = 10. Upper values correspond to the right mode andlower values to the left mode.Chapter 4. Nonlinear Numerical SimulationsFigure 4.11: Positions of waves for right mode only110400Horizontal position of maximum density at critical level for right wave (solid line) andminimum density at critical level for left wave (dashed line) when the flow is initiallyperturbed with the right mode oniy for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10.The horizontal position has been scaled with respect to its period 27r/a*.0 100 200 300Chapter 4. Nonlinear Numerical Simulations 111:0.0—cycle1•0:0.5 =0 0c3 0.0—-o ocycle1.08=0.250.5 Eto’ 0.0- 00+30cycle1.00.5-ci- 0-0—0.5- 00 0-1.0 tt.t.t0 10 20 30cycleFigure 4.12: Average phase speeds from nonlinear simulationsAverage phase speed over a cycle for: 1. nonlinear simulations with right or left modeonly (+), 2. nonlinear simulations with both modes (0) when J = 0.3, R = 3, Pr = 9,Re = 300, and H = 10. Phase speed predicted by linear theory is indicated by a dashedline. Positive phase speeds correspond to right moving waves and negative wave speedsto left moving waves.Chapter 4. Nonlinear Numerical Simulations 1121.00.5 —0.00 100 200 3000.04000—0.02 I ......... I0 100 200 300 45Figure 4.13: Energy budgets for right mode only with e = 0Energy budget for nonlinear simulation initially perturbed with the right mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, and H 10 with e = 0. Reservoirs: T + 6.1 (solidline), K—2.1 (dashed line), K’ (dotted line), and F+9.1 (dash-dot-dash line). Transfers:C(K, K’) (solid line), and C(F, K’) (dotted line). Diffusions: D(K) (dashed line), D(K’)(dotted line), and D(F) (dash-dot-dash line).2.52.01.540000.004 —0.002—0.0000100 200 300 400Chapter 4. Nonlinear Numerical Simulations 1132.5 -2.0 — -,,1.5 -0 70 1400.02 I0.000 70 1400.004 —0.002 —0 -0.000— —0 70 140tineFigure 4.14: Energy budgets for right mode oniy with e = 0.5Energy budget for nonlinear simulation initially perturbed with the right mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, and H = 10 with e = 0.5. See figure 4.13 fordescription.Chapter 4. Nonlinear Numerical Simulations 114Figure 4.15: Evolution of mean profiles for right mode onlyMean velocity and density profiles at 1. t = 0 (dashed line), 2. from numerical simulationsinitialized with right mode only (solid line), and 3. diffusion of mean flow withoutperturbations (dotted line). Parameter values are J = 0.3, R = 3, Pr = 9, Re = 300, andH=l0.U4U Up—1.0 —0.5 0.0 0.5 1.0pChapter 4. Nonlinear Numerical Simulations 115Figure 4.16: Density and vorticity for left mode only with e = 0.25Results of nonlinear simulations when the flow is initialized with the left mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and e = 0.25. Contours are of constantdensity p Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled withrespect to its period 27r/&.0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0Chapter 4. Nonlinear Numerical Simulations 116Figure 4.17: Density and vorticity for left mode only with = 0.5Results of nonlinear simulations when the flow is initialized with the left mode only forJ = 0.3, R = 3, Pr = 9, Re = 300, H = 10, and e = 0.5. Contours are of constantdensity p = pa(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled withrespect to its period 2ir/&.Chapter 4. Nonlinear Numerical Simulations 117400Figure 4.18: Positions of waves for left mode oniyHorizontal position of maximum density at critical level for right wave (solid line) andminimum density at critical level for left wave (dashed line) when the flow is initiallyperturbed with the left mode only for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10.The horizontal position has been scaled with respect to its period 21r/a*.0 100 200 300Chapter 4. Nonlinear Numerical Simulations 118Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3,R = 3, Pr = 9, Re 300, H 10, and e = 0. Contours are of constant densityp p(z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respectto its period 27r/&.Figure 4.19: Density and vorticity for both modes with E = 0Chapter 4. Nonlinear Numerical Simulations 119Figure 4.20: Density and vorticity for both modes with e = 0.1Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3,R = 3, Pr = 9, Re = 300, H = 10, and = 0.1. Contours are of constant densityp = Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respectto its period 27r/.0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.Chapter 4. Nonlinear Numerical Simulations 120Figure 4.21: Density and vorticity for both modes with 0.25Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3,R = 3, Pr = 9, Re = 300, H = 10, and = 0.25. Contours are of constant densityp = Pa(Z) + p’ in the x-z plane and shading is of vorticity. x has been scaled with respectto its period 27r/&.0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0Chapter 4. Nonlinear Numerical Simulations 121Figure 4.22: Density and vorticity for both modes with e = 0.5Results of nonlinear simulations when the flow is initialized with both modes for J = 0.3,R = 3, Pr = 9, Re = 300, H = 10, and E = 0.5. Contours are of constant densityp Pa(Z) + p’ iii the x-z plane and shading is of vorticity. x has been scaled with respectto its period 27r/&.0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0Chapter 4. Nonlinear Numerical Simulations 122400Figure 4.23: Positions of waves for both modesHorizontal position of maximum density at critical level for right wave (solid line) andminimum density at critical level for left wave (dashed line) when the flow is initiallyperturbed with both modes for J = 0.3, R = 3, Pr = 9, Re = 300, and H = 10. Thehorizontal position has been scaled with respect to its period 27r/a*.100 200 300Chapter 4. Nonlinear Numerical Simulations 1230 100 200 300 4000.02000—0.02 ——0.04 I I3 100 200 300 4000.004J/J\/r J\/ 1/J\1/\0.002 —0.000 ——0 100 200 300 400Figure 4.24: Energy budgets for both modes with E = 0Energy budget for nonlinear simulation initially perturbed with both modes for J = 0.3,R = 3, Pr = 9, Re = 300, and H = 10 with e = 0. See figure 4.13 for description.Chapter 4. Nonlinear Numerical Simulations0.40.20.001240.0020>a)(I)a)0.00010 20 30 40 5000.00040.0002a-)—0.0000C)F——0.0002—0.000410 20 30 40 500Figure 4.25: Wave positions and energy budgets in linear regime for e = 0Wave positions and energy budgets in linear regime for 0 as described in figure 4.24.Times of maximum wave separation and wave intersection are indicated by vertical lines.1.00.80.610 20 30 40 50(I)0U)D0.0040.0020.0000 10 20 30 40 50tChapter 4. Nonlinear Numerical SimulationsFigure 4.26: Wave positions and energy budgets in nonlinear regime for e = 0Wave positions and energy budgets in nonlinear regime for e = 0 as described in figure 4.24. Times of maximum wave separation and wave intersection are indicated byvertical lines.1.00.80.60.40.20.0200125250210 2202.52.0230 2401.51.0 - - - -0.5Eo.0_ -.200 210xU)0>0(I)00.040.020.00—0.02—0.040.0040.0020.000220 230 240 250-.U)G)U)C0F-U)C0U)D200 210 220 230 240 250200 210 220 230 240 250tChapter 4. Nonlinear Numerical Simulations 1262.5 -2.0 —1.5 -0.0• - -0 70 1400.04 I0.02 —0.00 - - - - - -. . - - - - - - - - - -.‘ - - - - - - - - - - - - - - - - - - ‘.- - - - - - - - - - - - -- - - -- - S.. -- - - - -—0.02 ——0.04 I0 700.0040.002 —0.000 — -- —0 70 140Figure 4.27: Energy budgets for both modes with e = 0.5Energy budget for nonlinear simulation initially perturbed with both modes for J = 0.3,R = 3, Pr = 9, Re = 300, and H = 10 with = 0.5. See figure 4.13 for description.Chapter 4. Nonlinear Numerical Simulations 1278O.2521.0 -1.0-0.5 0.51.0Figure 4.28: Evolution of mean density profiles for both modesMean velocity and density profiles at 1. t = 0 (dashed line), 2. from numerical simulations initialized with both modes (solid line), and 3. diffusion of mean flow withoutperturbations (dotted line). Parameter values are J = 0.3, R = 3, Pr = 9, Re = 300, andH=10.U4—1.0 —0.5 0.0 0.5 1.0 —1.0 —0.5 0.0 0.5 1.0p p0—2—1.0—0.5 0.0 0.5pChapter 5Comparison with Lab oratory ExperimentsAlthough numerical simulations provide a vast amount of quantitative information abouta flow, we are using it to study an isolated phenomenon. In nature, such a phenomenonis unlikely to occur on its own. Also, due to numerical restrictions, it is often difficult, ifnot impossible, to reproduce the high Reynolds number characteristic of natural systems.These issues lead to the question of whether or not the results of the numerical simulationsare indeed physically realistic.There have been field observations of Kelvin-flelmholtz billows, for example, in cloudformations, Colson [11], and in the Mediterranean Sea, Woods [69]. There is, however,a paucity of quantitative data on Holmboe instabilities in the field and in the laboratory. Also, it is difficult to make field observations of Holmboe instabilities due to thenon-stationary nature of the instabilities. Thus, although we would like to make detailed quantitative comparisons between the numerical results and both laboratory andfield observations, we restrict our attention to qualitative comparisons with laboratoryobservations.In order to compare our results with those of laboratory experiments, we draw on threesources: the results of Lawrence et al. [38], Zhu [72], and Guez &? Lawrence [20]. First,we examine the results of Zhu who provides one of the best sequences of photographs of anearly symmetric Holmboe instability. This instability occurs at a relatively large valueof the bulk Richardson number. Next, we compute the evolution of two flows with weakstratification and large density interface displacement and compare these with the results128Chapter 5. Comparison with Laboratory Experiments 129of Lawrence et al. [38] and Guez & Lawrence [20]. Iii all three cases, the stratificationis a result of having a layer of fresh water over a layer of salt water. Full descriptions ofthe experimental setups are given in the individual references.5.1 Symmetric CaseThe results from Zhu [72] are shown in figure 5.1. It is evident that there is a waveprotruding into the upper layer which is moving to the right and one into the lower layermoving to the left. Since the mean flow is moving to the right, it appears that the rightwave is moving faster than the left wave. Flow conditions for this experiment are given intable 5.1. From the data, we calculate the relative speeds of right and left moving waveswith respect to the mean flow. The relative phase speeds are 0.6 cm/s and -0.6 cm/sfor right and left moving waves, respectively. The fact that the two waves are movingin opposite directions at the same speed relative to the mean flow is indicative that thedensity interface is at the same level as the center of the shear layer (i.e., = 0). Sincethe measured waves speeds are approximate, these wave speeds do not contradict themeasured range of as it is still possible that the density interface displacement is smalland negative. For the purpose of comparison, however, we set e = 0 in the numericalsimulation.Due to the large values of both the Reynolds and Prandtl numbers (table 5.1), wecannot match the flow conditions in the numerical simulations since doing so would resultin insufficient diffusion to adequately resolve the flow. Also, selecting 10 R 15 withRe = 460 requires a large number of grid points in the z-direction. For the purpose ofcomparison, we use the same parameters used in chapter 4 to examine the evolutioll ofthe symmetrie case. The values of these parameters are J 0.3, Re 300, Pr = 9,R = 3, and = 0.54. Parameters relating to the numerical scheme are given in table 4.1.Chapter 5. Comparison with Laboratory Experiments 130J Re Pr g’ Ui U2 h A R c$!)cm/s2 cm/s cm/s cm cm cm/s cm/s0.45 460 700 1.6 2.5 -1.7 5 14 [-0.1,0] [10,15] 1.0 -0.20.06 25 700 0.4 3.7 1.7 0.58 3.9 0.5 10 0 *0.03 25 700 0.13 3.1 1.5 0.6 4.1— 10 0 *Table 5.1: Approximate flow conditions of laboratory experiments from Zhu [72](J = 0.45), Lawrence et al. [38] (J 0.06), and Guez & Lawrence [20] (J = 0.03).A is the observed wavelength of the instability. g’ is the modified gravitational acceleration. When J = 0.03 and 0.06 there are no observable left moving waves (indicated by *in the table).The lower value of the initial bulk Richardson number used in the numerical simulationsis not an unreasonable choice since J increases as the flow evolves. Also, there is someuncertainty in determining the value of J in the laboratory experiments.In order to qualitatively compare our results with those of Zhii [72], we must outputour results at the same frequency as the photographs were taken. During the 0.5 secondbetween photographs, a wave has travelled (cr — U) 0.5/A 0.02 of its total wavelength.In non-dimensional time, we must output our results at intervals of t = 0.02 A*/c =0.02 . 2n/c 0.389. Here * indicates the non-dimensional quantities used in thenumerical simulation.Results of the numerical simulation for J = 0.3 with e = 0 are shown in figure 5.2. Asin the experimental observations of Zhu [72], there is a wave moving to the right whichprotrudes into the upper layer and one moving to the left protruding into the lower layer.There are, however, several qualitative differences between the two results. The wavepeaks are sharper in the laboratory results of Zhu [72] (reproduced in figure 5.1) thanill the results of the numerical simulations. A possible explanation for this difference isthat the density interface is much thinner for the experimental flow than in the numericalChapter 5. Comparison with Laboratory Experiments 131simulation. Also, it is difficult to determine at what stage of the perturbation evolutionthe flow shown in the photographs occurs. As seen in figure 4.19, the shape of the waveschanges throughout the evolution of the perturbation.Another observable difference is that in the experiments, the amplitude of the leftmode appears larger than the amplitude of the right mode, whereas in the numericalsimulations they have the same amplitude. There are two possible explanations for thisdifference. First, as shown in figure 4.20, the two waves can have significantly differentamplitudes even when the density interface displacement is small. Thus it may be possiblethat the density interface shift is a small negative number. This possibility is supportedby the range of e given in table 5.1 for the resilits of Zhu [72]. Second, in chapter 4 it wasdiscovered that the relative amplitudes of right and left modes during early stages of theirnonlinear evolution depends on the relative initial amplitude of the two modes. Sincethe Holmboe instability observed by Zhu [72] ocdllrred in a facility designed to studyan exchange flow through a channel with an underwater sill, the Holmboe instability isjust one of several mechanisms that are taking place (see figure 7.la of Zhii [72]). It ispossible that one of the other mechanisms has an influence on the amplitude of the twomodes. In spite of these differences, we feel that the comparison between the numericaland experimental results are favourable.5.2 Non-Symmetric CaseExperimental results from Guez & Lawrence [20] and Lawrence et al. [38] are shownin figures 5.3 and 5.4 respectively. Approximate flow conditions are given in table 5.1.These laboratory observations are clear examples of flows exhibiting “one-sidedness”. Inboth cases, perturbations take the form of fluid being dragged from the lower layer intothe upper layer. This behaviour is characteristic of the non-symmetric case with > 0.Chapter 5. Comparison with Laboratory Experiments 132Since, for both cases, the bulk Richardson number is small, billowing is observed andincreases as J decreases.Although this set of experiments occurs at a lower Reynolds number, the Prandtlnumber and the ratio R are still too large to numerically model. The lower value of theReynolds number does, however, allow us to easily compute the evolution of the flowfor several values of R (and hence several values of Pr) to see if any trends are present.For this purpose, we use the flow parameters from the experimental results of Guez &Lawrence [20]. Since the value of e is not given for this experiment, we set e = 0.5. Thisvalue is given by Lawrence et al. [38]. Since experiments of Guez & Lawrence [20] wereperformed in the same facility as those of Lawrence et al. [38], it is not unreasonable toassume that the density interface displacements are approximately equal.When J = 0.03 and e = 0.5, the right moving instability is the only unstable modeand is the only mode contained in the initial perturbation. At low Reynolds numbers,the initial amplitude of the perturbation has a significant effect on how the flow develops.If the initial amplitude is too small, then diffusion acts to stabilize the flow before theperturbation has a chance to fully develop. This phenomenon was observed by Klaassen &Peltier [29] in their study of Kelvin-Helmholtz instabilities. In the experiments of Guez &Lawrence [20], the flow is forced at a frequency corresponding to a = 0.45, allowing us toinitialize the perturbations with a relatively large amount of kinetic energy (K’ = 0.005).To compare our results with the photographs in figure 5.3, we wish to output theresults at intervals equal to the time required for the wave to travel one wavelength.This time interval, t0 is given byA= U+Cr2irh/21— (U+c.6U)where c is the non-dimensional phase speed obtained from linear stability analysis. InChapter 5. Comparison with Laboratory Experiments 133non-dimensional time, we obtain— irh 5’U— &(U + c 6U) h/227r6U&(U + c 6U) (5.2)= 4.75Results from the numerical simulations with R = 3, 6, and 10 (and hence Pr = 9,36, and 100, respectively) are shown in figures 5.5, 5.6 and 5.7. In all cases, nx = 64,nz = 128, and At = 0.0625. Although the results are shown at intervals of 4.75, ascalculated above, the billows nevertheless move across the frame as time increases. Thismovement is explained by recalling that, as discussed in Chapter 4, the phase speed of aright moving wave increases with time when e 0.For all values of R, the flow behaves qualitatively the same. A finger of heavierfluid is drawn up into the upper layer and forms a billow. The finger becomes thinnerwith a sharper tip as R (and hence Pr) increases. For R = 6 and 10, the billow hasseparated by t = 57. As the billow wraps around itself, it is stretched and becomesthinner. Eventually, the upper portion of the billow becomes so thin that it diffuses intothe ambient fluid. At this point, the billow divides into two disconnected parts. WhenR = 3, the billow remains connected at all times since the initial thickness of the billowis much larger than that of the other two cases. It is observed that the diffusion rateincreases as R decreases. This change in diffusion rate can be related to the decreasingvalue of the Prandtl number and is consistent with the results of Klaassen & Peltier [28]who studied the effect of the Prandtl number on Kelvin-Helmholtz billows. For R = 10we observe a resolution problem at the selected grid size. Since the results for R = 6 arequalitatively the same as for the R = 10 case and do not change when the grid is refined,we use this case for comparison with the experimental results of Guez & Lawrence [20].Although the results from the numerical simulations depicted in figure 5.6 and thoseChapter 5. Comparison with Laboratory Experiments 134of the laboratory experiments of Guez & Lawrence [20] reproduced in figure 5.3 arequalitatively the same, there are several differences between them. In the numericalsimulations the tips of the figures in the billows are not as sharp as those of the laboratoryexperiments. This difference between experimental and numerical observations is a resultof the values of R and Pr used in the numerical simulations being much smaller than thoseof the experiments. This difference becomes less pronounced as R and Pr are increased.Another noticeable difference between the two results, is that there is not as muchwrapping around of the billows in the numerical simulations as in the laboratory experiments. There are two possible causes of this discrepancy. Firstly, Klaassen & Peltier [28]found that as the Prandtl number increases, there is more wrapping around of fingers ina Kelvin-Helmholtz instability before the displaced fluid diffuses into the ambient flow.Since we have used Pr = 36, which is much smaller than the actual value, to model theflow, perhaps the large amount of diffusion suppresses the wrapping around of the billow.Secondly, there is some uncertainty in determining the flow conditions for the laboratoryexperiments. It is possible that the actual value of the bulk Richardson number is slightlysmaller than 0.03. Since the amount of billowing is sensitive the the value of the bulkRichardson number, we expect more curling up of the billows as J decreases.Breakdown of the billows occurs at approximately the same time in both the numericalsimulations and the laboratory experiments. The cause of this breakdown is different,however, for the two cases. In the laboratory experiment the billows collapse when theflow becomes three-dimensional. Since the numerical simulations restrict the flow totwo-dimensions, transition to three-dimensional flow cannot possibly be the cause of thebreakdown of the billows. Instead, it is due to the top part of the billow diffusing intothe ambient fluid. Because of this difference, it does not make physical sense to comparethe two flows beyond this point. Nevertheless, it is of interest to note that the two flowsremain qualitatively the same.Chapter 5. Comparison with Laboratory Experiments 135Next we examine the experimental results of Lawrence et al. [38] when J = 0.06 (figure 5.4). Since the results of the nonlinear simulations with R = 6 compared favourablywith the experimental results of Guez & Lawrence [20] and both experiments were conducted in the same facility, we use the same value to compare with the results of Lawrenceet al. [38]. The only difference is that no forcing is present in the laboratory study ofLawrence et al. [38]. Nevertheless, we found that a fairly large initial amplitude of theperturbation was still required in order for the perturbation to become fully developed.Having a large initial amplitude poses no problems since only the right moving mode isan unstable mode when J = 0.06 with = 0.5. When both right and left modes are unstable, more care must be taken in choosing the initial amplitudes of the perturbations,as is discussed below.Results of the numerical simulations with J = 0.06 are shown in figure 5.8 (here= 4.95). As with the J = 0.03 case, the fingers of the billows are not as thin as in theexperimental results of Lawrence et al. [38] (figure 5.4). Also, the fingers do not curl upas much in the numerical simulations. The possible explanations for these differences arethe same as the J = 0.03 case. Although the results of the numerical simulations and thelaboratory experiments differ at later times, qualitative agreement as the perturbationinitially develops is excellent.Lawrence et al. [38] examined the development of perturbations for a range of Richardson numbers (figure 5.4). Comparing numerical results with those of Lawrence et al. [38]is difficult when J is large and both right and left moving modes are unstable. The problem lies in selecting the correct initial amplitude of the perturbation. If the amplitudeis not sufficiently large, then diffusion damps the perturbation before it has a chance togrow. If the initial amplitude is too large, the perturbation does not spend much time inthe linear regime, resulting in the left mode being of approximately the same amplitudeChapter 5. Comparison with Laboratory Experiments 136as the right mode. The experimental observations indicate, however, that the initial amplitude is such that the perturbation spends sufficient time in the linear regime for theright mode to dominate. Finding the correct initial amplitude is a matter of trial anderror. We have not attempted to numerically reproduce these experiments.Chapter 5. Comparison with Laboratory Experiments 137Sequence of photographs showing Holmboe waves (figure 7.lb of Zhu [72]). The upperlayer is from left to right, the lower layer from right to left. The grid markings are 5cm apart and the photographs were taken at 0.5 second intervals. See table 5.1 for flowconditions.Figure 5.1: Experiments of Zhu [72]CD CjD CIIC C).CD c_— I!-4. CD CD -4 Ci) C CD Ci) C Cl)p 0 p p bI I I.Chapter 5. Comparison with Laboratory Experiments 139Figure 5.3: Experiments of Guez & Lawrence [20]Sequence of photographs showing development of Kelvin-Helmholtz type billows fromGuez & Lawrence [20]. The flow is from left to right. The field view of each photographis approximately 25cm with an overlap of approximately one wavelength between eachsuccessive photograph. See table 5.1 for flow conditions.Chapter 5. Comparison with Laboratory Experiments 140II II II II IIFigure 5.4: Experiments of Lawrence et al. [38]Sequence of photographs showing development of shear instabilities at various Richardsonnumbers (figure 6 from Lawrence et al. [38]). The flow is from left to right. The field viewof each photograph is approximately 25cm with an overlap of approximately 5cm betweeneach successive photograph. As the bulk Richardson number increases, less billowing isobserved. See table 5.1 for flow conditions when J = 0.06.IIInNJ0NJNJ0NJ0p 0 p NJp 0 p NJ p pC CCD0CD CcJIICDII:J) CCCDcn.CDo->0 C))c±a-.o-—1CDCD IICD CDCCDCi) CC.CD C)1 CD Ci) Ci) C CD CD U)0 0 p aC,N) C, 0 00 p pIIp•NJ0NJ•0 p NJ p pNJ0NJ•NJ0NJ•C),0 0 C 0 N) CI, 0 0p 0 0 p p bN) 9) 0 0 0 CC) 0o•p NJ P P 0 0p 0 P NJID b C C CoSNJ0NJSp 0 S 0 0.i[...1K.,.P 0 SN) C CPChapter 5. Comparison with Laboratory Experiments 1429cqo•o t= 14.2500 t= 19.0000 .=.75004 4 4 4oo2o:4o:6 0:8 1.0= 28.5000 t 33.2500 t = 38.0000 t = 42.75004 4 4 4. .•.:04 08. •1.O47.5000 t = 52.2500 t = 57.0000 t = 61.7500= 66.5000 t = 71.2500 t = 76.0000 t = 80.75004 4 4 40.0.26081.0 o.o:20o.8l.o o.oo:24’o:6o81.oFigure 5.6: Results from numerical simulations for J = 0.03 with R = 6Density contours of results from numerical simulations with J = 0.03, Re = 25, R = 6,Pr = 36, = 0.5, with & 0.45. Contours are from -0.72 to 0.72 in equal increments of0.18. x has been scaled with respect to its period 27r/&.Chapter 5. Comparison with Laboratory ExperimentsFigure 5.7: Results from numerical simulations for J = 0.03 with R = 10143Density contours of results from numerical simulations with J = 0.03, Re = 25, R = 10,Pr = 100, e = 0.5, with & = 0.45. Contours are from -0.72 to 0.72 in equal increments of0.18. x has been scaled with respect to its period 27r/&.t = 950000 t = 14.2500 I = 190000 = 23.7500420—2—404204I = 38.00001.041.0= 42.750042—40 0 0.2 0.4 0.6 0.8= 33.250042 0.2 0.4 0.6 0.8= 28.5000—2—4 —40.0—40.2 0.4 0.6 2.842461.75001.00.0 0.2 0.4 0.6 0.8 1.0=__47.500040.0 0.2 0.4 0.6 0.8 1.0= 57.00000—2—44200.2 0.4 0.6 0.8 1.0= 66.50004= 71.250043 0.2 0.4 0.6 0.8 1.0t = 80.75000.2 0.4 0.6 0.8 1.0 C420—2—40.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0Chapter 5. Comparison with Laboratory Experiments 144= 0.00000 t = 4.93750 t = 9.87500 t = 14.81254 4 4 42 2 2 2_- — ---- z—4 —4 —4—40.00.2 0.4 0.60.8 1.0 0.0 0.2 0.60.81.0 0.0 0.2 0.4 0.6 1.0 0.00.20.40.60.81.019.7500 87 29.6250 t= 34.56254 4 4 40:6 1.0 0!839.5000 t = 44.4375 t = 49.3750 t = 54.31254 4 4 4O.OO:2 0:6 o!4 0:6= 59.2500 t 64.1875 t = 69.1250 t = 74.06254 4 4 41.0Figure 5.8: Results from numerical simulations for J 0.06Density contours of results from numerical simulations with J = 0.06, Re = 25, R = 6,Pr = 36, e = 0.5, with & = 0.45. Contours are from -0.72 to 0.72 in equal increments of0.18. x has been scaled with respect to its period 27r/&.Chapter 6Summary and ConclusionsWe have studied the effect of displacing a thin density interface with respect to thecenter of the shear layer on the stability of a parallel flow. Linear stability analysis wasused to determine how increasing the density interface displacement changed the overallstability of the flow. In particular, we concentrated on how this displacement affectsHolmboe instabilities. Results of both inviscid theory and viscous theory were considered.Numerical simulations were then performed to examine the nonlinear evolution of theseinstabilities. Finally, results of the nonlinear simulations were compared with those oflaboratory experiments. We only considered displacements of the density interface belowthe center of the shear layer.6.1 Linear Stability AnalysisWhen the density interface is displaced with respect to the center of the shear layer,results of linear stability analysis indicate that right and left moving modes, which haveequal growth rates when the density interface displacement is zero, split into a fastergrowrng mode in one direction and a more slowly growing mode moving in the otherdirection. When the effects of horizontal boundaries are negligible and the density interface is placed below the center of the shear layer, it is the right moving wave that hasthe larger growth rate.145Chapter 6. Summary and Conclusions 1466.1.1 Inviscid ResultsInviscid theory was used to study the linear stability of piecewise linear backgroundprofiles. This has the advantage of admitting an exact dispersion relation. Also, otherthan the density interface displacement, the only additional parameter of this flow isthe bulk Richardson number. For the non-symmetric case, decreasing the height of thevertical domain damps the growth rate of the right mode more quickly than the growthrate of the left mode. This difference in rate of decrease for the two modes results in theleft mode having the larger growth rate when the total height of the domain is less thanor equal to five times the shear layer thickness. When the density interface correspondsto the center of the shear layer (the symmetric case), the results of Smyth [51] werereproduced with transition from Kelvin-Helmholtz to Holmboe instabilities occurring atsmaller values of the bulk Richardson number as the domain height decreases.For pure flolmboe waves, there is a region where the maximum growth rate increaseswith increasing Richardson number, indicating that three-dimensional instabilities arepossibly more unstable than two-dimensional instabilities. It was found, however, that itis indeed the two-dimensional instabilities that dominate. For non-symmetric Holmboeinstabilities, the maximum growth rate of the right mode decreases monotonically asthe Richardson number increases and so there is no possibility of the three-dimensionalinstabilities dominating. The left mode still has a region for which the growth rateincreases with the Richardson number, but, as with the symmetric case, two-dimensionalinstabilities remain more unstable.6.1.2 Viscous ResultsSince the nonlinear numerical simulation requires that there are no discontinuities in theinitial data or its derivatives, results of linear stability analysis for a viscous flow withChapter 6. Summary and Conclusions 147hyperbolic tangent background profiles were examined. Three additional parameterswere introduced: the Reynolds number, the Prandtl number, and the ratio of shear layerthickness to density interface thickness, R. As the density interface displacement wasincreased, the precise behaviour of right and left modes depended on the value of R.For large finite values of R (R > 4), the behaviour was similar to the discontinuousdensity interface case which corresponds to R = cc. The most unstable mode of theright moving wave occurs at a smaller wave number whereas the maximum growth rateof the left mode occurs at a larger wave number. When R < 4, however, there is achange in this behaviour. The most unstable right wave now occurs at a larger value ofthe wave number. The limb on the stability diagram corresponding to the left movingwave continues to tilt towards the wave number axis, but also separates itself from theaxis. This behaviour becomes more pronounced as R decreases.As with the inviscid case, the effect of decreasing the domain height was examined.For the symmetric case, the behaviour was significantly different than the inviscid case.Instead of the transition from Kelvin-Helmholtz to Holmboe instabilities being shiftedto a smaller Richardson number as the domain height decreases, the unstable modesnear this transition are the first to be stabilized. This is a direct result of the modes inthis region being the least unstable. When the total domain height is less than or equalto five times the height of the shear layer thickness, there is a stable region separatingKelvin-Helmholtz instabilities from Holmboe instabilities. For the non-symmetric case,the growth rate of the right mode decreases more rapidly than that of the left mode asthe boundaries get closer to the shear layer. Unlike the inviscid case, however, the rightmode always remains more unstable than the left mode.For R = 3 with Pr = 9 and Re = 300, Smyth & Peltier [55] have shown that threedimensional instabilities are more unstable than two-dimensional instabilities when thedensity interface displacement is zero. As the density interface displacement increases,Chapter 6. Summary and Conclusions 148two-dimensional instabilities dominate for the faster growing right mode. The left mode,however, still exhibits characteristics necessary for three-dimensional instabilities to bedominant and thus may be a mechanism through which perturbations become three-dimensional. Further investigation is required to determine if three-dimensional instabilities are indeed dominant for the left moving wave. As R increases, three-dimensionalinstabilities are less likely to be primary instabilities as the results tend to those of thepiecewise linear case.Studying the shapes of the eigenfunctions predicted from linear stability analysis provided insight into the growth mechanism of the instabilities. For the symmetric case,Smyth & Peltier [54] showed that when both right and left modes are present, the instabilities oscillate between growth and decay phases of equal duration. During the growthphase, the critical level acts as a source, extracting energy from the mean flow and thedensity interface as a sink, converting perturbation kinetic energy to potential energy.The same behaviour is observed for the non-symmetric case. One important difference,however, is that, as the density inlerface displacement increases, the right mode is responsible for a larger portion of the energy exchange between the mean kinetic energy,the perturbation kinetic energy, and the potential energy reservoirs.When the density interface displacement is small (0 0.25), the growth anddecay of the perturbations are in phase with and governed by the transfer of energybetween the mean flow and the perturbations. The maximum in the perturbation kineticenergy occurs just after the waves have passed through each other. The minimum occursjust after the waves have reached the stage of maximal separation and are starting toapproach each other. When the density interface displacement is large (E = 0.5), thegrowth and decay of the perturbations are still governed by the exchange of energybetween the mean flow and the perturbations, but the two quantities are no longer inphase.Chapter 6. Summary and Conc1u.ions 1496.2 Nonlinear SimulationsNumerical simulations were used to study the nonlinear evolution of the perturbations. Inorder to examine the effect of initial conditions on the development of the perturbationsand to study nonlinear interactions between the right and left moving waves, we ranthree sets of simulations. We first perturbed the flow with the more unstable right mode.Next, we tracked the evolution of a flow initially perturbed with the more slowly growingleft mode. Finally, we perturbed the flow with both right and left modes.For the symmetric case, the same behaviour is observed when the flow is initiallyperturbed with either the right or left mode. The mode that is initially present dominatesthe flow at early stages of the evolution. Eventually, the other mode grows and the flowtends to a state where both modes are present. When the perturbation initially containsboth modes, they grow at the same rate and the flow remains symmetric throughoutits evolution. Thus, although the early behaviour depends on which modes are used toinitialize the flow, long term behaviour appears to be independent of initial conditions.When the flow is initially perturbed with the right mode only, the amount of timespent in the linear regime decreases as the density interface displacement increases. Thistrend is due to the increased growth rate of the right mode, as predicted by linear stabilityanalysis. The left mode does eventually develop, but takes longer to become noticeableas the density offset is increased. Once again, the behaviour can be attributed to theresults of linear stability analysis which predict a smaller growth rate for the left modeas the density offset is increased.When the perturbations initially contain only the left mode, the flow spends moretime in the linear regime as the density offset increases. Also, the right mode affects theflow more quickly. For very large density offsets, the right mode affects the flow whileit is still in the linear regime and dominates the flow at an early stage in its evolution.Chapter 6. Summary and Conclusions 150The behaviour of the flow when both modes are initially present is not surprising in viewof the above results. When the density offset is large, the right mode quickly dominatesthe flow. The left mode, however, continues its slow growth and becomes noticeable atlater times in the simulation. One problem encountered was our inability to resolve theseparation of a fluid parcel from the perturbation. Resolution problems occurred whenthe density offset was large and restricted the time for which the simulation was valid.Thus little can be said about the long term behaviour of such flows.Growth rates in the linear stage of the simulations were compared to the results oflinear theory. Excellent agreement was found for the symmetric case. For the non-symmetric case, the right mode grows more slowly and the left more quickly than thepredicted values. When the flow has both modes present, the overall linear growth ratelies between those of the right only and left only cases, but is dominated by the rightmode.Initial phase speeds agree with those predicted from linear stability analysis. As timeincreases, the phase speed of the right moving wave increases. That of the left movingwave increases for small density offsets and decreases for large density offsets. Thisbehaviour indicates that perhaps all waves are tending to the same phase speed in thenonlinear regime. When oniy one mode is present, there are no oscillations in the phasespeed. When both modes are present, however, the phase speeds of the waves change asthey pass through each other. In the nonlinear regime, the waves speed up as as theyapproach each other and slow down as they pass each other. Also, there is a decrease inthe average phase speed of the waves when both waves are present.In contrast to the linear regime where oscillations in the perturbation kinetic energyare governed by the exchange between the mean flow and the perturbations, in thenonlinear regime, oscillations in the perturbation kinetic energy are closely linked to theexchange of energy between the potential energy reservoir and the perturbation kineticChapter 6. Summary and Conclusions 151energy reservoir. Maxima in the perturbation kinetic energy occur when the waves arepassing through each other, whereas the minima occur when the waves are maximallyseparated.6.3 Comparison with Laboratory ExperimentsResults of the numerical simulations were compared with the experimental observationsof Zhu [72] (J 0.4, e 0), Guez & Lawrence [20] (J 0.03, e 0.5), and Lawrence etal. [38] (J 0.06, e 0.5). Due to numerical considerations, we were unable to match thehigh Prandtl number and thin density interface thickness of the salt stratified laboratoryexperiments. To determine the effect of the Prandtl number on the development ofinstabilities, we ran a series of numerical simulations with several values of the Prandtlnumber. The density interface thickness used in the simulations are related the Prandtlnumber by R = /N, where 1/R is a measure of the density interface thickness. As thePrandtl number increased and the density interface thickness decreased, the waves orbillows formed by the instabilities became more sharply defined.Good qualitative agreement was obtained between the results of the numerical simulations and those of the laboratory experiments. Main differences were due to the thickerdensity interface and the lower Prandtl number used in the numerical simulations. Sinceobtaining quantitative information about a flow from experimental observations is oftendifficult, simulations may be of use to quantify the instabilities studied here.6.4 Future WorkAlthough the above study has allowed us to gain better insight into the nonlinear development of non-symmetric Holmboe instabilities, many issues have been brought up andremain unresolved. Some of the more important ones are mentioned here.Chapter 6. Summary and Conclusions 152Due to numerical restrictions, we were unable to resolve the long term behaviour ofthe flow. This behaviour is of interest as it appears that the eventual state of the flow isindependent of the initial conditions. Also, computing the long term evolution of the flowwould allow one to verify the conjecture that the phase speed of the waves asymptotesto a single value. When the density offset was large, denser fluid from the lower layerseparated from the wave and was released into the less dense upper layer. Since we wereunable to consistently compute the evolution of this phenomenon, it is desirable to findout the time and manner in which it occurs. Resolving this mechanism would also permitthe study of long term behaviour of these instabilities.The maximum amplitude of the perturbation kinetic energy appears to be independent of the initial conditions. Since the perturbation kinetic energy contains the totalcontribution of both modes, determining the individual contributions of right and leftmodes would allow one to better study the interaction between the two waves.The nonlinear simulations examined here forced right and left modes to have the samewave length. For the non-symmetric case, linear stability analysis indicates that the mostunstable right and left modes occur at different wave numbers. It would be of interestto see how this would affect the evolution of the instabilities.Finally, we restricted our study to the two-dimensional development of instabilities.Although this is often valid for the initial development of the perturbations, it is welldocumented that the flow eventually becomes three-dimensional. Hence determining themanner in which it becomes three-dimensional is of importance.List of Symbols(j horizontally averaged (mean) qnantity, unless otherwise noted 74(.) vertically integrated quantity 74(“) discrete Fourier transform of a quantity 76() differentiation with respect to time 76a nondimensional wave number 13real part ofa 14& wave number used in numerical simulation 74wave number at which the maximum growth rate occurs 307 coefficient of expansion for stratifying agent 8two-dimensional Laplacian operator 8tIh discrete two-dimensional Laplacian operator 75Lip maximum density variation 9Lit time step size in numerical simulation 77size of horizontal space step in numerical simulation 75Liz size of vertical space step in numerical simulation 756 delta function 25density scale used for non-dimensionalizing 12SU velocity scale used for non-dimensionalizing 12amplitude of perturbations 23e nondimensional vertical density interface displacement 25‘7 density interface thickness 318 concentration of stratifying agent 8153List of Symbols 154concentration of stratifying agent at standard state 8angle with respect to x-axis of three-dimensional perturbations . .. .29it diffusivity of stratifying agent 8A dimensional wave length 130A amplification factor of semi-discretized scheme 77p molecular viscosity 8ii kinematic viscosity 11p density of fluid 8complex vertical amplitude of p’ from normal mode analysis 13p’ deviation of density from hydrostatic equilibrium 9Pa vertical density distribution of fluid in hydrostatic equilibrium 9Fourier coefficient for 3 32Po density at the standard state 8P1 mean density in upper layer 3P2 mean density in lower layer 3approximate value of p’(iAx, j/Xz, t) from numerical simulation... 76complex vertical amplitude of i/” from normal mode analysis 13Fourier coefficient for 32nth stream function amplitude function for wave number a 15x phase angle of perturbation kinetic energy 47initial mean stream function distribution 13approximate value of “(iAx, jz, t) from numerical simulation .. . 761’ stream function 11perturbation of stream function from steady state 13w vorticity 79List of Symbols 155complex vertical amplitude of vorticity ( — iaii) 45A matrix obtained from Galerkin’s method to solve linear problem .. 33A coefficients to complete solution to linear problem 15a th coefficient of dispersion relation for piecewise linear profiles .. .26G(R1,R2) conversion of energy from reservoir R1 to reservoir R2 81c complex nondimensional wave speed 13th wave speed associated with the wave number a 15Cg group velocity 14c, imaginary part of c 13c. real part of c 13complex wave speed associate with right moving mode 28complex wave speed associate with left moving mode 28D(R1) loss of energy from reservoir R1 due to diffusion 81discrete first derivative in x-direction 75discrete second derivative in x-direction 75discrete first derivative in z-direction 75discrete second derivative in z-direction 75d vertical density interface displacement 25g gravitational acceleration 8g’ modified gravitational acceleration 130h shear layer thickness 25H total depth of flow 25imaginary part of a complex value 15J bulk Richardson number 3J value of J below which Holmboe instabilities do not occur 39List of Symbols 156JT transition from Kelvin-Helmholtz to Holmboe waves occurs at JT . 27K kinetic energy of mean flow 80K’ perturbation (or wave) kinetic energy 74K’ filtered perturbation kinetic energy 80L length scale used for non-dimensionaliziug 12N number of Fourier modes used to solve linear problem 32nx number of grid points in x-direction for numerical simulation 83nz number of grid points in z-direction for numerical simulation 83F potential energy 80Pr Prandtl number 12p pressure 8deviation of pressure from hydrostatic equilibrium 9complex vertical amplitude of p’ from normal mode analysis 29Pa vertical pressure distribution of fluid in hydrostatic equilibrium .... 9R ratio of shear layer thickness to density interface thickness 31real part of a complex value 13Re Reynolds number 12Sc Schmidt number 12Ri gradient Richardson number 19T total energy 81t time 8time interval between outputs for nonlinear numerical simulation . 130U initial mean horizontal velocity distribution 24(I average horizontal velocity between upper and lower layers 12U1 meau velocity in upper layer 3List of Symbols 157U2 mean velocity in lower layer 3u velocity vector 29ü complex vertical amplitude of u’ from normal mode analysis 29u’ perturbation velocity vector 29u horizontal component of fluid velocity 8ü complex vertical amplitude of u’ from normal mode analysis 45streamwise component of velocity perturbation 23instantaneous perturbation of horizontal velocity from mean flow .74V vector containing Fourier coefficients q. and pn 33spanwise component of velocity perturbation 29w vertical component of fluid velocity 8complex vertical amplitude of w’ from normal mode analysis 29vertical component of velocity perturbation 23x horizontal coordinate in streamwise direction 8y horizontal coordinate in spanwise direction 29z vertical coordinate 8Bibliography[1] P. BALDwIN AND P. H. ROBERTS, The critical layer in stratified shear flow, Mathematika, 17 (1970), pp. 102—119.[2] W. E. BOYCE AND R. C. DIPRIMA, Elementary Differential Equations and Boundary Value Problems, John Wiley & Sons, New York, 3rd ed., 1977.[3] P. BRADSHAw, An Introduction to Turbulence and its Measurement, PergamonPress, Oxford, 1971.[4] F. K. BR0WAND AND Y. H. WANG, An experiment on the growth of small disturbances at the interface between two streams of different densities and velocities, inProceedings of the International Symposium on Stratified Flows, Novosibirsk, SovietUnion, August 29- 311972.[5] F. K. BR0wAND AND C. D. WINANT, Laboratory observations of shear-layer instability in a stratified fluid, Boundary-Layer Meteorology, 5 (1973), pp. 67—77.[6] C. CANUTO, M. Y. HuS5AINI, A. QuATER0NI, AND T. A. ZANG, Spectral Methods in Fluid Mechanics, Springer-Verlag, Berlin, 1988.[7] C. P. CAuLFIELD, Stratification and Buoyancy in Geophysical Flows, PhD thesis,Cambridge University, Cambridge, England, 1991.[8] C. P. CAuLFIELD AND W. R. PELTIER, Three dimensionalization of the stratifiedmixing layer, Physics of Fluids, 6 (1994), pp. 3803—3805.[9] S. CHANDRASEKHAR, Hydrodynamic Hydromagnetic Stability, Oxford UniversityPress, Oxford, 1961.[10] B. W. CHAR, B. LE0NG, K. 0. GEDDEs, M. B. MONAGAN, G. H. GONNET,AND S. M. WATT, Maple V Language Reference Manual, Springer-Verlag, NewYork, 1991.[11] D. CoLsoN, Wave-cloud formation at Denver, Weatherwise, 7 (1954), pp. 34—35.[12] P. G. DRAzIN AND W. H. REID, Hydrodynamic Stability, Cambridge UniversityPress, Cambridge, 1981.158Bibliography 159[13] C. B. EMMANUEL, B. R. BEAN, L. G. MCALLISTER, AND J. R. POLLARD,Observations of Helmholtz waves in the lower atmosphere with acoustic sounder,Journal of the Atmospheric Sciences, 29 (1972), pp. 886—892.[14] H. J. S. FERNANDO, Turbulent mixing in stratified fluids, Annual Review of FluidMechanics, 23 (1991), pp. 455—493.[15] H. G. FIsHER, E. J. LIsT, R. C. Y. KoR, J. IMBERGER, AND N. H. BROOKS,Mixing in Inland and Coastal Waters, Academic Press Inc., San Diego, 1979.[16] K. S. GAGE, Linear viscous stability theory for stably stratified shear flow: A review,Boundary-Layer Meteorology, 5 (1973), pp. 3—17.[17] M. GAsTER, A note on the relation between temporally-increasing and spatially-increasing disturbances in hydrodynamic stability, Journal of Fluid Mechanics, 14(1962), pp. 222—224.[18] E. E. GOSSARD AND J. H. RICHTER, The shape of internal waves of finite amplitude from high-resolution radar sounding of the lower atmosphere, Journal of theAtmospheric Sciences, 27 (1970), pp. 971—973.[19] M. C. GREGG, Diapycnal mixing in the thermocline: A review, Journal of Geophysical Research, 92 (1987), pp. 5249—5286.[20] L. G. GuEz AND G. A. LAwRENCE, Hydrogen bubble tracking and the determination of velocity and vorticity fields in a stratified mixing layer. Submitted toExperiments in Fluids.[21] P. HAzEL, Numerical studies of the stability of inviscid stratified shear flows, Journalof Fluid Mechanics, 51(1972), pp. 39—61.[22] H. VON HELMHOLTz, On discontinuous movements of fluids, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 36 (1868),pp. 337—346.[23] J. H0LMB0E, On the behavior of symmetric waves in stratified shear layers, Geofys.Publ., 24 (1962), pp. 67—113.[24] L. N. HOWARD, Note on a paper of John W. Miles, Journal of Fluid Mechanics, 10(1961), pp. 509—512.[25] W. KELVIN, Hydrokinetic solutions and observations, The London, Edinburgh, andDublin Philosophical Magazine and Journal of Science, 42 (1871), pp. 362—377.Bibliography 160[26] G. H. KEuLEGAN, Interfacial instability and mixing in stratified flows, Journal ofResearch of the National Bureau of Standards, 43 (1949), pp. 487-500.E27] G. P. KLAAssEN, The Transition to Turbulence in Stably Stratified Parallel Flows,PhD thesis, University of Toronto, Toronto, Canada, 1982.[28] G. P. KLAAssEN AND W. R. PELTIER, The effect of Prandtl number on the evolution and stability of Kelvin-Helmholtz billows, Geophysical and Astrophysical FluidDynamics, 32 (1985), pp. 23—60.[29] , Evolution of finite amplitude Kelvin-Helmholtz billows in two spatial dimensions, Journal of Atmospheric Sciences, 42 (1985), pp. 1321—1339.[30] , The onset of turbulence in finite-amplitude Kelvin-Helmholtz billows, Journalof Fluid Mechanics, 155 (1985), pp. 1—35.[31] , The role of transverse secondary instabilities in the evolution of free shearlayers, Journal of Fluid Mechanics, 202 (1989), pp. 367—402.[32] C. G. Koop, Instability and turbulence in a stratified shear layer, Tech. ReportRep. USCAE 134, Univ. Southern Calif., Dept. Aerospace Ellg., 1976.[33] C. G. KooP AND P. K. BROWAND, Instability and turbulence in a stratified layerwith shear, Journal of Fluid Mechanics, 93 (1979), pp. 135—159.[34] C. G. Koop, C. D. WINANT, AND F. K. BR0wAND, 1995. Personal communication.[35] D. K0PPEL, On the stability of flow of a thermally stratified fluid under the actionof gravity, Journal of Mathematical Physics, 5 (1964), pp. 963—982.[36] H. LAMB, Hydrodynamics, Cambridge University Press, Cambridge, 6th ed., 1932.[37] G. A. LAwRENcE, 1995. Personal communication.[38] G. A. LAWRENCE, F. K. BR0WAND, AND L. G. REDEK0PP, The stability of asheared density interface, Physics of Fluids A, 3 (1991), pp. 2360—2370.[39] S. A. MAsL0wE, Shear flow instabilities and transition, in Hydrodynamic Instabilities and the Transition to Turbulence, H. L. Swinney and J. P. Gollub, eds.,Springer-Verlag, Berlin, 2nd ed., 1985, pp. 181—228.[40] T. MAxw0RTHY AND F. K. BR0WAND, Experiments in rotating and stratifiedflows: oceanographic application, Annual Review of Fluid Mechanics, 7 (1975),pp. 273—305.Bibliography 161[41] J. W. MILES, On the stability of heterogeneous shear layers, Journal of Fluid Mechanics, 10 (1961), PP. 496—508.[42] A. S. M0NIN AND A. M. YAGL0M, Statistical Fluid Mechanics, The MIT Press,Cambridge, Massachusetts, 1971.[43] S. NIsHIDA AND S. YOSHIDA, Stability and eigenfunctions of disturbances in stratified two-layer shear flows, in Proceedings of the Third International Symposium onStratified Flows, Pasadena, California, 3—5 February 1987 1987.[44] P. C. PATNAIK, F. S. SHERMAN, AND G. M. CoRcos, A numerical simulation ofKelvin-Helmholtz waves of finite amplitude, Journal of Fluid Mechanics, 73 (1976),Pp. 215—240.[45] J. PEDL0SKY, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1979.[46] C. PELLAcANI, Shear instabilities in stratified fluids; linear theory, La Rivista DelNuovo Cimento, 6 (1983), pp. 1—26.[47] S. POND AND G. L. PICKARD, Introductory Dynamical Oceanography, PergamonPress, New York, 2nd ed., 1983.[48] 0. P0uLIQuEN, J. M. CH0MAz, AND P. HuERRE, Propagating Holmboe waves atthe interface between two immiscible fluids, Journal of Fluid Mechanics, 266 (1994),pp. 277—302.[49] LORD RAYLEIGH, On the stability, or instability, of certain fluid motions, Proc.London Math. Soc., 11(1880), pp. 57—70.[50] F. S. SHERMAN, J. IMBERGER, AND G. M. CoRcos, Turbulence and mixing instably stratified waters, Annual Review of Fluid Mechanics, 10 (1978), PP. 267—288.[51] W. D. SMYTH, Hb1mboe waves, Master’s thesis, University of Toronto, Toronto,Canada, 1986.[52] , Instability and Turbulence in Density-Stratified Mixing Layers, PhD thesis,University of Toronto, Toronto, Canada, 1990.[53] W. D. SMYTH, G. P. KLAASSEN, AND W. R. PELTIER, Finite amplitude Holmboewaves, Geophysical and Astrophysical Fluid Dynamics, 43 (1988), Pp. 181—222.[54] W. D. SMYTH AND W. R. PELTIER, The transition between Kelviri-Helmholtz andHolmboe instability: An investigation of the overreflection hypothesis, Journal ofAtmospheric Sciences, 46 (1989), pp. 3698—3720.Bibliography 162[55] , Three-dimensional primary instabilities of a stratified, dissipative, parallel flow,Geophysical and Astrophysical Fluid Dynamics, 52 (1990), pp. 249—261.[56] , Instability and transition in finite-amplitude Kelvin-He lmholtz and Holmboewaves, Journal of Fluid Mechanics, 228 (1991), Pp. 387—415.[57] H. B. SQUIRE, On the stability of three-dimensional disturbances of viscous flowbetween parallel walls, Proceedings of the Royal Society A, 142 (1933), pp. 621—628.[58] C. STAQUET, Two-dimensional secondary instabilities in a stably-stratified shearlayer, in Proceedings of the 4th International Symposium on Stratified Flows, Grenoble, France, 29 June — 2 July 1994.[59] G. STRANG, Introduction to Applied Mathematics, Wellesley-Cambridge Press,Cambridge, MA, 1986. pp. 453—456.[60] H. U. SvERDRUP, M. W. JOHNsON, AND R. H. FLEMING, The Oceans, TheirPhysics, Chemistry, and General Biology, Prentic-Hall, Inc., Englewoods Cliffs, N.J.,1942.[61] H. TANAKA, Quasilinear and nonlinear interactions offinite amplitude perturbationsin a stably stratified fluid with hyperbolic tangent shear, Journal of the MeteorologicalSociety of Japan, 53 (1975), pp. 1—31.[62] G. I. TAYLOR, Effect of variation in density on the stability of superposed streamsof fluid, Proceedings of the Royal Society A, 132 (1931), pp. 499—523.[63] S. A. THORPE, A method of producing a shear flow in a stratified fluid, Journal ofFluid Mechanics, 32 (1968), pp. 25—48.[64] D. J. TRITT0N AND P. A. DAvIEs, Instabilities in geophysical fluid dynamics, inHydrodynamic Instabilities and the Transition to Turbulence, H. L. Swinney andJ. P. Gollub, eds., Springer-Verlag, Berlin, 2nd ed., 1985.[65] J. S. TURNER, Buoyancy Effects in Fluids, Cambridge University Press, Cambridge,1973.[66] Y. H. WANG, An experimental study of the instability of a stably stratified free shearlayer, Journal of Fluid Mechanics, 71(1975), pp. 563—575.[67] R. C. WEAST, ed., CRC Handbook of Chemistry and Physics, CRC Press, Inc.,Boca Ratoii, Florida, 1st student ed., 1988.[68] F. M. WHITE, Viscous Fluid Flow, McGraw Hill, New York, 2nd ed., 1991.Bibliography 163[69] J. D. WOODS, Wave-induced shear instability in the summer thermocline, Journalof Fluid Mechanics, 32 (1968), pp. 791—800.[70] N. Y0NEMITsu, The Stability and Interfacial Wave Phenomena of a Salt WedgeFlow, PhD thesis, University of Alberta, Edmonton, Canada, 1991.[71] S. YO5HIDA, On a mechanism for breaking of interfacial waves, Coastal Engineeringin Japan, 20 (1973), pp. 7—15.[72] Z. Zilu, Exchange flow through a channel with an underwater sill, PhD thesis,University of British Columbia, Vancouver, Canada, 1995.Appendix ADispersion Relation for Piecewise Linear ProfilesIn section 3.2 we considered an inviscid flow with the following nondimensional background velocity aild density profiles.1—1 z>—ePa =(1It was shown that the normal modes approach to linear stability produces the nondimensional Taylor-Goldstein equation.(U— c)( — 2) — U + 2J(z + = 0, (A.2)For the background flow described by (A.1), (A.2) can be solved exactly. Except atz = —1, —e, and 1, the Taylor-Goldstein equation reduces to—= 0. (A.3)Thus we can write the solution to (A.2) asA1e+ B1e’— A2e+ B2eA3e+ B3eA4e+ B4e’1U= z—ll<z<H/2—l<z<l—H/2<z<—1(A.1)1<z<H/2—E < Z < 1—l <z < —E—H/2<z<—1(A.4)164Appendix A. Dispersion Relation for Piecewise Linear Profiles 165In order to determine the coefficients A1 through B4, we impose the following boundary conditions and matching conditions. We require q = 0 at z +H/2 givingAie°72+Bie”2 = 0,(A.5)A4e”2+B4e’2 = 0.Next, we require that q be continuous, givingAie+Bie = A2e+B,A2e6+ B2e = A3e+ B3e, (A.6)A3e+ B3ea = A4e+ B4ea.To derive the matching conditions, we rewrite the Taylor-Goldstein equation (A.2) as[(U— c)— Uz] — 2(U — c) + 2J6(z + e)u0 (A.7)On integrating (A.7) across a discontinuity z from z — to z + and letting —* 0, weobtain the following jump conditions.Atz=+1Z\ [(U — c)qz — Uzq] = 0, (A.8)and at z = —e[(U— c)—U] + 2JU = 0. (A.9)Imposing the jump conditions (A.8) and (A.9) we obtain the following equations for A1through B4.(—1 + c) [a(Ai — A2)e — a(Bi — B2)e] = A1e + B1e,(e + c) [a(A3 — A2)e + a(B2 — B3)e] = (A2e+B2e), (A.10)(—1— c) [a(A3 — A4)e — a(B3— B4)ej = A3e+ B3ea.Appendix A. Dispersion Relation for Piecewise Linear Profiles 166Equations (A.5), (A.6), and (A.10) form a complete set of equations for A1 throughB4 and can be written in the formA1A2A3A4M = 0 (A.11)B1B2B3B4The dispersion relation is obtained by setting det M = 0. Using Maple [10] (a symbolic programming language) to simplify calculations, the following dispersion relationis obtained.C4 + a1c + a2c + a3c + a4 = 0 (A.12)wherea1 28, (A.13)Ja2 = 82c(e2— 1) (1 + e2— 2e° cosh 2ae)1+4a2(e— 1) {(i + 2a) — e4 — e2 [(i — 2a) — e4j+4e(sinh2a — 2ccosh2c)}, (A.14)—J sinh(2) ( e2 H_fl + e2 — 2e)a3 =— 1)2(e — 1) {_(l + 2a) + e4 + e2 [(1— 2)2— e4j+4e(— sinh 2c + 2a cosh 2)}, (A.15)andAppendix A. Dispersion Relation for Piecewise Linear Profiles 167_2a4 =42(e2H— 1) {_(1 + 2a) + e4 + e2 [(1— 2a)— e4j+4e(—sinh2c + 2ccosh2a)}—J43(e2— 1) {_(i +2)2— e4 + e2 — 2a) — e4j+2 [2e(2a2 — 1) +e2_1)(1 — 2a) +e2(2 + 1)] cosh2e+4e(cosh 2c — 2c sinh 2a)}. (AJ6)Appendix BMaximum Growth Rates for Holmboe’s EquationIn chapter 3, we determine the value of the bulk Richardson number, J, at which transition from Kelvin-Helmholtz to Holmboe instabilities occurs for the two-layered floworiginally considered by flolmboe [23]. In addition, we examine the possibility that instabilities are initially three-dimensional. To this end, we must determine the maximumgrowth rate ac for a fixed value of the bulk Richardson number. Note that Smyth etal. [53] found the maximum growth rate for a fixed wave-number.For the case studied by Holmboe we set 0 and let 11 —* cx in the dispersionrelation obtained in appendix A for the piecewise linear profile described by (3.9). Weobtain the following dispersion relation,c4+a2=0, (B.1)where—J a+a_a2= 2 (B.2)a4 = (B.3)1 — 2o + e2a= 2 (B.4)Since we are interested in finding the most amplified growth rate crc, for a fixed valueof J, we define o- = —icrc to be the complex growth rate. The above dispersion relationbecomesoA + b2o + b4 = 0, (B.5)168Appendix B. Maximum Growth Rates for Holmboe’s Equation 169with= aJ+aa_, (B.6)aJa2. (B.7)Equation (B.5) yields two types of instabilities when J > 0: Holmboe instabilitieswhen b — 4b < 0, and Kelvin Helmholtz instabilities when b — 4b > 0 and b2 < 0. Wewish to address the following question: given a bulk Richardson number, J, at what wavenumber a does the most unstable mode occur (i.e., the largest value of ar)? Due to thenature of the different instabilities, it is necessary to consider the two types separately.B.1 Holmboe InstabilitiesWhen the flow is unstable to Holmboe instabilities b— 4b < 0 and equation (B.5) hasfour roots of the form a = +a,. + ia. Solving for a2 in equation (B.5), one can write—b2 + i/4b4 —(a + ia)2 = (B.8)2Equation B.8 can be used to derive the following equation for ar:a + + -(b — 4b) = 0 (B.9)In order to determine the maximum growth rate and the wave number at which thisgrowth rate occurs, we differentiate (B.9) with respect to a and set 3ar/9 0. Theresulting equation defines the curve in the a-J plane along which r is maximized:—— 0 (B.10)Rewriting the above in terms of a and J we obtain—2a”2J+(1 — 6a — e + 4ae_2) Ju/2 +ah/2(2 —4a — 2e4a) 0 (B.11)A1 A2 A3Appendix B. Maximum Growth Rates for Holmboe’s Equation 170Since it is a quadratic in J1/2, we find that, solving for J, there are two possible solutions.The solution determining the curve along which the growth rate is maximized for a fixedJ is given by_________— (—A2 —— 4A1A3—2A ) (B.12)It should be noted that this solution is only valid in the region of Holmboe instabilities.To determine the minimum value of (cr, J) for which (B.12) is valid, one must solve— 4b = 0 with J given by (B.12). The minimum value (, J) = (0.488, 0.007) isobtained.B.2 Kelvin-Helmholtz InstabilitiesFor Kelvin-Helmholtz instabilities b — 4b > 0 and b2 < 0. In this case, solutions toequation (B.5) are of the form u = ar. Thus we can use the above methods on (B.5)directly. The curve along which the growth rate is maximized is defined by:—B + — 4BB3 (B.13)2B1whereB1 16a2(7 — 1), (B.14)B2 8a — 32a2 + + (2 — 32a + 96o2 — 48)7 +(—4 + 32a — 88c2 + 32a)72+ (2 — 8a + 32a2)7, (B.15)B3 —1 + 8o — 24a + 32€ — 16a4 + (1 — 6a + 20a — 40€ + 32a4)’y +(2— 10cr + 8a2 + 83)72 + (—2 + 6c — 12a2 + 16a3)-y +(—1 + 2 + 8a2)y4+-y5, (B.16)and = e2a.Appendix B. Maximum Growth Rates for Holmboe’s EquationFigure B.1: Holmboe’s [23] stability diagramContours are of constant growth rate, act. Dashed line indicate values of a for which thegrowth rate is maximum for a given J.171By examining the sign ofô2ur/öa it is possible to determine the range of a for whichthe growth rate along the curve given by equation (B.13) is indeed a maximum. Thecurves defined by (B.12) and (B.13) are shown in figure B.1.0.500.400.30-)0.200.100.000.0 0.5 1.0 1.5cc2.0Appendix CEnergy EquationsHere we derive the energy equations used for diagnostics of the results from linear stabilityanalysis and from the nonlinear simulations. The equations derived here are only validfor flows with the same boundary conditions as the flow under present consideration.These are essentially the same equations as those derived by Klaassen [27] in his studyof Kelvin-Helmholtz instabilities. For a more general description of energy equations forturbulent flows see IVlonin & Yaglom [42].C.1 Kinetic EnergyWe start with the nondimensional form of the momentum equations with the Boussinesqapproximation (2.12) and (2.13).uu + wu = —p + (C.1)Wt + nw + ww —p — Jp’ + Aw. (C.2)Assuming that p constant (a valid approximation when 6p << p0), the equation governing the evolution of the kinetic energy is found by multiplying (C.1) by u and (C.2)by w and then adding the two together. We obtainUUt + WWt + u(uu + ww) + w(uu + ww)=—(up + wp) — Jwp’ + (uu + wAw). (C.3)Equation (C.3) is the general form for the evolution of kinetic energy in a flow. The172Appendix C. Energy Equations 173last term in (C.3) is usually divided into a viscous diffusion term and viscous dissipationterm:nu + w’w = + w2) _[(u)2 + (u)2 + (w)2 + (w)2] (C4)viscous diffusion viscous dissipationFor simplicity in the following derivations, however, we leave these lumped together anddeal with this issue as needed.Since we wish to separate the kinetic energy of the mean flow from that of the perturbations, we writeu(x, z, t) = ü(z, t) + u’(x, z, t), (C.5)w(x,z,t) w’(x,z,t), (C.6)where, as defined in chapter 4, (j= (j dx)/L, where L is the length of the domain inthe x-direction, indicates that the quantity has been horizontally averaged. Substitutingthe above expressions for u and w into (C.3), and taking the horizontal average, we obtainthe following.1 (2) + () + () + u (2+ w’w ++ + ‘‘u + + + ui(C7)If we define (u’2 + w’2)/2 to be the perturbation kinetic energy and ü2/2 the meankinetic energy, then (C.7) describes the evolution of the sum of the two kinetic energies.In order to separate the mean and perturbation kinetic energies, we consider the equation governing the conservation of momentum in the x-direction. Taking the horizontalaverage of (C.1) we get(C.8)Appendix C. Energy Equations 174which when multiplied by ü gives the evolution equation for the mean kinetic energy:(u) + u+ uw’ — -uu. (C.9)To determine the equation for the horizontally averaged perturbation kinetic energy,we subtract (C.9) from (C.7), which gives(u!2 + + a + + u’2 + u’w’w + u’w’u + w’2 += ———JZ77 + -- (u’Au’ + w’.’). (C.1O)An excellent description of the terms in (C.1O) for the homogeneous case can be foundin Bradshaw [3].Equations (C.9) and (C.1O) can be simplified by using the periodicity in the xdirection. It is easily shown that= w’w = 0. (C.11)Also,u’p + w’p = u’p + (u + w’jp’ + w’p, by incompressibility,= (u’p’) H- (w’p’), (C.12)= (w’p’), by periodicity in x.Using the above simplifications, the equation for the evolution of the mean kinetic energycan be written as(c.13)and the equation for the evolution of the horizontally averaged perturbation kineticenergy becomes1(I2 + !2) + + + w’2 +=— (w’p’) —J+ (u’u’ + w’w’). (C.14)Appendix C. Energy Equations 175Equation (C.14) is used in chapter 3 to help understand the linear evolution of the flow.In order to understand how the mean flow and perturbations interact in the simillations described in chapter 4, we study the energy budgets for the horizontally averagedand vertically integrated mean kinetic energy, perturbation kinetic energy, and potentialenergy. Vertically integrating (C.13) and (C.14) we obtainK + (ü7) = —(uu), (C.15)andK + (u’w’w) + (u’w’u) + (w’w’w) + (a)=— ((w’p’)) — J() + (u’u’ + w’w’), (C.16)where K = (u2), K’ = (u’2 + w’2), and (.) = f—H/2 •dz.Equations (C.15) and (C.16) can be simplified using the boundary conditions ‘ == p’ = 0 at z = +H/2 (or equivalently w’ = 0 at z = +H/2) as well as theperiodicity in the x-direction. It is evident that(w’w’wØ = ((w’p’)) = 0. (C.17)Also,K u’w’w + u’w’u) = —(w’2u+ u’2w), integration by parts,= (w’2+ u’2L’) by incompressibility, (C 18)= (w’2w), by periodicity in x,= 0 w’=Oatz=+H/2.Appendix C. Energy Equations 176and____1 L H/2(üw’u)= 71 1 üw’udzdx,L H/2LL IH/2Z, integration by parts,0 at z +H/21 L H/2L IH/2 (u’üw’ — u’üu) dzdx, by incompressibility, (C.19)1 L H/2= L I H/2 üu’w’dxdz by periodicity in x,Finally,Ku’u’ + w’w’) = + w!2)) — ((u)2 + (u)2 + (w)2 + (w)2)= _((u)2 + (u)2 + (w)2 + (w)2) (C.20)and(aa) = - (()2) (()2) (C.21)Thus (C.15) and (C.16) become— (C.22)andK=-(u) - J() - + (u)2 + (w)2 + (w)2). (C.23)C.2 Potential EnergyAt a given point in the fluid domain, the potential energy per unit volume is given bygpz, or, in non-dimensional form, Jpz. In order to compute the potential energy of theAppendix C. Energy Equations 177flow, we consider the nondimensional form of the heat equation (2.9),Pt + UPx + WPz RePr (C.24)Separating the density into a mean and perturbation,p(x, z, t) — (z, t) + p’(x, z, t) (C.25)and using similar expansions for u and w, (C.24) becomesPt + Pt + UP + + WPZ + Wp= RePr + + pj. (C.26)Multiplying by Jz to obtain the non-dimensional potential energy, taking the horizontalaverage, and integrating vertically gives(Jzt) + J(z(u’p + w’p))= RePr (C.27)We can simplify (C.27) using the following:(z(u’p + w’p)) = (zu’p + (u + w)p’ + w’p), by incompressibility,= (z(u’p’) + (w’p’)),= (z(w’p’)), by periodicity in x, (C.28)= —(7), integration by parts,w’ = 0 at z = ±H/2.Equation (C.27) simplifies toJ(wp) + RePr’ (C.29)where P = J(z) is the total potential energy of the flow. Unlike the kinetic energy, thepotential energy cannot be separated into mean and perturbation parts. Instead, thepotential energy combines the effects of both the mean flow and the perturbations.Equations (C.22), (C.23), and (C.29) are used in chapter 4 to compute the energybudget of the flow as the perturbations evolve.Appendix DPeriod of Perturbation Kinetic EnergyWhen 0, linear stability analysis predicts that the phase speed of right and left modesare different. Using the normal mode approach we can write the velocity perturbationas the sum of the right and left modes, i.e.,z, t) = u)(x, z, t) + u’(x, z, t), (D.l)= G(r)(z)cosx(r) + S(T)(z)sinx + C(z)cosx +S(1)(z)sinx(, (D.2)where (‘1 and (l) are the components corresponding to the right and left modes, respectively and (r) = — c)t), = — c1)t). For simplicity we consider only thehorizontal velocity component u’, but the results are valid when the vertical velocitycomponent w’ is also present. The periods of the right and left modes are given byT(r)=27r/o4r) and T1 27r/oc!). When the two modes are superimposed, the periodof the resulting wave is given byT = lcm(T,T), (D.3)the lowest common multiple of the two periods.In chapter 3, we silperimpose the two modes and examine how the horizontally averaged perturbation kinetic energy changes with time, giving a better understanding of thelinear evolution of the flow. Here, we compute the period of the horizontally averagedperturbation kinetic energy.178Appendix D. Period of Perturbation Kinetic Energy 179= 1JLU,2dx= 1+ u’(1)2dx,= 1 jL ) + s(r) ) + c( + x]dx,= [c( + S(r)2 + +2o+2 (C(l)G(r) + cos a(c —+ 2 (C(r)S(l) + C(1)S(r)) sin (e — c1))t], (D.4)where L = 2ir/a is the wavelength of the flow (i.e., the period in the x-direction). From(D.4), the period of the horizontally averaged perturbation kinetic energy is 2ir/a(cT)—
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Non-symmetric Holmboe waves
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Non-symmetric Holmboe waves Haigh, Susan Patricia 1994
pdf
Page Metadata
Item Metadata
Title | Non-symmetric Holmboe waves |
Creator |
Haigh, Susan Patricia |
Date Issued | 1994 |
Description | When two flows of different velocity and density meet, a shear layer with a density gradient is formed. Under certain conditions this flow can be unstable. A statically stable stratified shear flow in which the density interface is much thinner than the shear layer thickness can be linearly unstable to two modes with equal growth rates and equal and opposite phase speeds. The superposition of these two modes is called a Holmboe instability. This instability is only possible when the flow is symmetric about the center of the shear layer. We examine the effect of breaking this symmetry by allowing the center of the density interface to be displaced with respect to the center of the shear layer. There are three major components to this study: linear stability analysis, nonlinear numerical simulations, and comparison with laboratory experiments. Linear stability analysis is used to examine the effect of the density interface offset on the overall stability of the flow. Both inviscid theory with piecewise linear background velocity and density profiles, and viscous theory with smooth background profiles are used. As in previous studies, it was found that the growth rate of one mode increases and that of the other mode decreases as the density interface displacement is increased. The precise behaviour depends on the relative thickness of the density interface with respect to the shear layer thickness. For inviscid theory with piecewise linear background profiles, we show that the initial perturbations must be two-dimensional. When the effects of viscosity and diffusion are included, however, it may be possible that the weaker mode is initially three-dimensional. Detailed analysis of the energy transfer in the linear regime indicates that when the background flow loses its symmetry, it is the mode with the larger growth rate that is primarily responsible for the extraction of energy from the mean flow. Two-dimensional numerical simulations are used to examine the nonlinear development of “non-symmetric” Holmboe instabilities. We start by perturbing the flow with the stronger mode predicted by linear theory. We then examine the response of the flow to weaker mode. Finally, we impose both modes. By comparing the development of these three flows we are able to study the interaction between the two modes and the effect of initial conditions on the development of instabilities. Although the initial development of instabilities depends on the initial conditions, this dependence weakens as the density interface offset is increased. Also, preliminary results indicate that long term behaviour of the perturbations are independent of initial conditions. Results of the numerical simulations are compared to both symmetric and nonsymmetric Holmboe instabilities that have been observed in laboratory experiments. Since we are unable to compute the flow at the high Prandtl number of the salt stratified experimental flows, we run the simulation for a range of increasing Prandtl numbers to determine how instabilities of thermally stratified flows (with low Prandtl number) and salt stratified flows differ. Results of these simulations indicate that differences between the experimental and numerical flows can be attributed to the thicker density interface and lower Prandtl number used in the numerical simulations. When the density interface is thicker and the Prandtl number lower, the waves or billows formed by the instabilities are not as sharply defined as in the experiments. |
Extent | 3261191 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080013 |
URI | http://hdl.handle.net/2429/7225 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-059707.pdf [ 3.11MB ]
- Metadata
- JSON: 831-1.0080013.json
- JSON-LD: 831-1.0080013-ld.json
- RDF/XML (Pretty): 831-1.0080013-rdf.xml
- RDF/JSON: 831-1.0080013-rdf.json
- Turtle: 831-1.0080013-turtle.txt
- N-Triples: 831-1.0080013-rdf-ntriples.txt
- Original Record: 831-1.0080013-source.json
- Full Text
- 831-1.0080013-fulltext.txt
- Citation
- 831-1.0080013.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080013/manifest