Vortex pairing and mixing in stratifiedshear flowsbyWenjing DongB.Eng., Tsinghua University, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Civil Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2017c© Wenjing Dong 2017AbstractKelvin-Helmholtz (KH) instabilities are an important source of mixing in oceans, lakes andthe atmosphere. The process of vortex pairing can increase the amount of mixing. First,the effects of initial conditions on vortex pairing and mixing are studied by running DirectNumerical Simulations with a variety of initial perturbations. It is shown that when thesubharmonic component of the perturbation is out of phase relative to the KH mode, vortexpairing is delayed or even eliminated. The amount of mixing in the simulations where thesubharmonic mode is out of phase is approximately half of that in the simulations wherethe subharmonic mode is in phase. The time of pairing is also found to be sensitive to thephase of the subharmonic mode. A slight change of the phase can change time of pairingsignificantly when the subharmonic mode is close to being out of phase. Second, the effectsof Prandtl number on KH instabilities, vortex pairing and mixing are studied. It is foundthat KH instabilities and vortex pairing are suppressed at higher Prandtl numbers, whichtends to reduce the amount of mixing. This effect is counteracted by enhanced three-dimensional motions in higher Prandtl number flows. However, the general trend is formixing and mixing efficiency to decrease as Prandtl number is increased.iiLay AbstractClimate and ocean circulation models rely on parameterization of turbulence. Understand-ing small scale turbulent events in the atmosphere, oceans and lakes helps optimize the pa-rameterization. The Kelvin-Helmholtz instability, a common instability causing small scaleturbulence in nature, is studied using numerical simulations. The objective is to investigatemixing of scalars, e.g. temperature and salinity, occurring in Kelvin-Helmholtz instabilitiesunder different circumstances. The effect of initial conditions is studied by running simula-tions with a variety of initial perturbations. Then the effect of molecular properties of thescalar, for example whether the scalar is temperature of salinity, on mixing is investigated.iiiPrefaceA version of chapter 2 has been submitted for publication. The authors are myself, E.W. Tedford, M Rahmani, and G. A. Lawrence. I developed the ideas under the guidanceof E. W. Tedford, M Rahmani, and G. A. Lawrence. I was responsible for preparing themanuscript and E. W. Tedford, M Rahmani, and G. A. Lawrence revised the manucript.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Linear modal stability theory . . . . . . . . . . . . . . . . . . . . . . 21.2.2 Linear non-modal stability theory . . . . . . . . . . . . . . . . . . . 51.3 Nonlinear evolution of KH instabilities . . . . . . . . . . . . . . . . . . . . 61.4 Mixing efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.5 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Sensitivity of vortex pairing and mixing to initial conditions . . . . . . 102.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.1 Mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.2 Direct Numerical Simulations . . . . . . . . . . . . . . . . . . . . . 142.2.3 Diagnostic tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 Pairing mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.4 Two-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4.1 Degree of modality . . . . . . . . . . . . . . . . . . . . . . . . . . . 19vTable of Contents2.4.2 Phase evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4.3 Lagrangian trajectories . . . . . . . . . . . . . . . . . . . . . . . . . 202.4.4 Growth rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4.5 Visualization of vertical velocity of the subharmonic component . . 232.5 Three-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.5.1 Kinetic energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.5.2 Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293 The effect of Prandtl number on KH instabilities . . . . . . . . . . . . . 313.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3 Two-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3.1 KH instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3.2 Vortex pairing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.4 Three-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.4.1 Vortex pairing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4.2 Kinetic energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4.3 Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.4.4 Mixing efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.5 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46viList of TablesTable 2.1 Numerical parameters for all the simulations. . . . . . . . . . . . . . 15Table 2.2 Table of all simulations and diagnostic variables. . . . . . . . . . . . . 17Table 3.1 Numerical parameters for all the simulations. . . . . . . . . . . . . . 33viiList of FiguresFigure 1.1 Growth rate as a function of wavenumber. . . . . . . . . . . . . . . . 6Figure 1.2 Snapshots of density field . . . . . . . . . . . . . . . . . . . . . . . . 7Figure 2.1 Vorticity snapshots of two-dimensional simulations. . . . . . . . . . . 13Figure 2.2 Conceptual drawing of vortex pairing mechanism. . . . . . . . . . . 18Figure 2.3 Evolution of r for random perturbation simulations. . . . . . . . . . 20Figure 2.4 Evolution of phase, particle trajectories, and evolution of growth rate. 21Figure 2.5 Normalized vertical velocity of the subharmonic component. . . . . . 24Figure 2.6 Kinetic energy of the subharmonic mode and three-dimensional motions. 26Figure 2.7 The increase of total and background potential energy. . . . . . . . . 27Figure 2.8 Saturation time of the subharmonic mode as a function of the phase. 29Figure 3.1 Kinetic energy of the KH and subharmonic modes. . . . . . . . . . . 34Figure 3.2 Thickness of the mean velocity and density profiles. . . . . . . . . . 35Figure 3.3 Density field for two-dimensional simulations. . . . . . . . . . . . . . 36Figure 3.4 Kinetic enregy of the subharmonic mode in 2D and 3D simulations. 37Figure 3.5 The evolution of the kinetic energy components. . . . . . . . . . . . 39Figure 3.6 The increase of total and background potential energy. . . . . . . . . 41Figure 3.7 Instantaneous and cumulative mixing efficiency. . . . . . . . . . . . . 42viiiAcknowledgementsFirst and foremost, I would like to thank my supervisor Greg Lawrence for his supportand guidance during my studies at UBC. Thanks him for being so kind to me and giveme much freedom about taking courses and doing research. Second, I would like to thankEdmund Tedford for always being so helpful. It is always a pleasure to talk and discussmy research with him. I also would like to thank Mona Rahmani for her guidance on myresearch. Special thanks goes to her help in using the numerical code for this thesis andrevising my paper. I also would like to thanks the members of EFM groups. Lastly butmost importantly, thanks to my parents.ixChapter 1Introduction1.1 MotivationNumerical modeling of environmental flows has been a valuable tool to understand complexflows in nature, for example atmospheric boundary flows over buildings or internal seichesin lakes. Subgrid-scale parameterizations are used in these models to describe small scaledynamic process that cannot be resolved due to computational constrains. Many closureschemes have been proposed to parametrize small scale turbulence, however the resultsgiven by different schemes can have non-negligible differences especially for inland waterbodies or near-shore regions (Ivey et al., 2008). A better understanding of the small scaleturbulence events in the atmosphere and ocean is therefore helpful to improve the accuracyof numerical models.In the atmosphere, ocean, and lakes, fluids are stratified due to temperature or salin-ity. Shear may exist in these stratified fluids due to wind, inflows, or topography. Whenshear exists, shear instabilities may arise, causing the transition to turbulence and enhanc-ing transport of scalars. There are three primary instabilities in stratified free shear flows,namely Kelvin-Helmholtz (KH) instabilities, Holmboe instabilities, Taylor-Caulfield insta-bilities. The velocity and density profiles determine which type of instabilities occurs. KHinstabilities are the most widely known among these three instabilities. For this reason, thisthesis is focused on KH instabilities, characterized by the roll-up of the density interface(Thorpe, 1971) and streamlines of cat’s eye shape (Drazin & Reid, 2004). After KH in-stabilities reach maximum amplitude, secondary instabilities grow causing the transition toturbulence. Most mixing occurs during this turbulent stage due to highly three-dimensionalsmall scale motions (Caulfield & Peltier, 2000). However, flow in the pre-turbulent stage canaffect mixing and mixing efficiency. When vortex pairing exists, it increases the Reynoldsnumber by merging KH billows and increasing the length scale of large scale motion (Moser& Rogers, 1993). Three-dimensionality which causes more mixing (Rahmani et al., 2014)is increased as a result. If turbulence is in the transition regime, mixing efficiency is alsoincreased (Shih et al., 2005). Therefore, we consider KH instabilities when vortex pairingexists.11.2. Theoretical background1.2 Theoretical backgroundMuch theoretical work on KH instabiities focuses on the linear stability theory which deter-mines whether an infinitesimal perturbation of specific wavenumber is stable by solving anequation for the perturbation. Linear modal theory introduces normal mode perturbationthat assumes growth rate is a constant to simplify the problems. Modal theory predicts themost unstable mode successfully for many hydrodynamic problems. However, the modaltheory only considers part of the possible perturbations and therefore non-modal stabilitytheory has been developed to find new unstable modes or understand non-modal growth(e.g. Arratia et al., 2013; Kaminski et al., 2014; Guha & Lawrence, 2014).1.2.1 Linear modal stability theoryLinear modal stability theory for inviscid and non-diffusive stratified shear flows was estab-lished by Taylor (1931) and Goldstein (1931) independently. The resulting governing equa-tion is called the Taylor-Goldstein (TG) equation. Miles (1961) and Howard (1961) provedthe gradient Richardson number Ri = N2/(∂U/∂z)2, where N defined as√−g/ρ0∂ρ¯/∂zis the Brunt-Va¨isa¨la¨ frequency, is less than 0.25 somewhere is a necessary condition forinstabilities to occur. In this section, we derive the viscous-diffusive TG equations, fromwhich TG equation can be obtained by neglecting diffusivity and viscosity.If we assume the fluid is incompressible and subject to the Bousinessq approximation,the governing equations for the total field variables u, p, and ρ are∇ · u = 0, (1.1)∂u∂t+ (u · ∇)u = − 1ρ0∇p− g ρρ0kˆ + ν∇2u, (1.2)∂ρ∂t+ (u · ∇)ρ = κ∇2ρ, (1.3)where kˆ is unit vertical vector, ρ0 is a reference density, ν is viscosity, and κ is diffusivity.Consider a horizontal parallel flow U(z, t) and an undisturbed stratification ρ¯(z, t) as thebackground flow. This flow is in equilibrium with the hydrostatic pressure p¯(z). SubstitutingU , ρ¯, and p¯ into equation (1.1), (1.2), and (1.3), we obtain the governing equations forU(z, t), ρ¯(z, t), and p¯(z, t), i.e.,∂U∂t= ν∂2U∂z2, (1.4)0 = − 1ρ0∂p¯∂z− g ρ¯ρ0. (1.5)21.2. Theoretical background∂ρ¯∂t= κ∂2ρ¯∂z2, . (1.6)Therefore, background shear U(z, t) and density ρ¯(z, t) are diffusing independently. Back-ground pressure p¯(z, t) can be calculated by integrating equation (1.5) vertically. Squires’theorem demonstrates that two-dimensional perturbation is more unstable than three-dimensional perturbation for shear instabilities in homogeneous shear flows. Koppel (1964)derived the governing equations for three-dimensional small perturbations in viscous anddiffusive fluids and proved that two-dimensional perturbations are sufficient to consider thestability of stratified shear flows. Hence we consider the stability of a two-dimensionalperturbation. The perturbed velocity, density, and pressure fields areu = U(z, t)i + u′(x, z, t), p = p¯(z, t) + p′(x, z, t), ρ = ρ¯(z, t) + ρ′(x, z, t), (1.7)where u′, p′, and ρ′ are perturbations. The equations for perturbations u′, p′, and ρ′ areobtained by substituting equation (1.7) into equation (1.1), (1.2), and (1.3) and subtractingthe equation (1.4), (1.5), and (1.6) for background flow. Dropping non-linear terms, weobtain the equations for infinitesimal perturbations, i.e.,∂u′∂x+∂w′∂z= 0, (1.8)∂u′∂t+ w′dUdz+ U∂u′∂x= − 1ρ0∂p′∂x+ ν∇2u′, (1.9)∂w′∂t+ U∂w′∂x= − 1ρ0∂p′∂z− g ρ′ρ0+ ν∇2w′, (1.10)∂ρ′∂t+ U∂ρ′∂x+ w′dρ¯dz= κ∇2ρ′. (1.11)Introducing stream function ψ′ defined byu′ =∂ψ′∂zand w′ = −∂ψ′∂x, (1.12)equation (1.9), (1.10), and (1.11) become∂2ψ′∂t∂z− ∂ψ′∂x∂U∂z+ U∂2ψ′∂x∂z= − 1ρ0∂p′∂x+ ν∇2∂ψ′∂z, (1.13)− ∂2ψ∂t∂x− U ∂2ψ′∂x2= − 1ρ0∂p′∂z− g ρ′ρ0− ν∇2∂ψ′∂x, (1.14)∂ρ′∂t+ U∂ρ′∂x− ∂ψ′∂x∂ρ¯∂z= κ∇2ρ′. (1.15)31.2. Theoretical backgroundIf we assume normal mode solutions of the following form,ψ′ = φˆ(z)eiα(x−ct), ρ′ = ρˆ(z)eiα(x−ct), and p′ = pˆ(z)eiα(x−ct), (1.16)where α is the wavenumber and is real because the flow is unbounded in x direction. c canbe complex and we denote it by cr + ici. The real part of c, cr, is phase speed, while theproduct of its imaginary part ci with α, αci, is growth rate. If ci is negative, the perturbationis stable and decreases exponentially. However, if ci is positive, the perturbation is unstbleand grows exponentialy. Substituting in normal mode solutions, equation (1.13), (1.14),and (1.15) become(U − c)dφˆdz− ∂U∂zφˆ = − 1ρ0pˆ+νiα(d2dz2− α2)dφˆdz, (1.17)α2(U − c)φˆ = −g ρˆρ0− 1ρ0dpˆdz− νiα( d2dz2− α2)φˆ, (1.18)(U − c)ρˆ+ ρ0gN2φˆ = κ1iα(d2dz2− α2)ρˆ. (1.19)Eliminating pressure pˆ in equation (1.17) , we obtain the equation for φˆ, i.e.,(U − c)( d2dz2− α2)φˆ− ∂2U∂z2φˆ =gρ0ρˆ+νiα(d2dz2− α2)2φˆ. (1.20)Equation (1.20) and equation (1.19) are viscous diffusive TG equations. They determinethe growth rate of each wavenumber mode. Note that the growth rate is in fact not aconstant as the shear layer and background density profile are diffusing. During a shorttime, we can assume viscous diffusive TG equations are correct for a constant c sincediffusion occurs slowly. Also, these two equations cannot be reduced into one equation forφˆ due to diffusivity.If we assume the fluid is inviscid and non-diffusive, equation (1.20) and (1.19) could bereduced to one equation for φˆ only, i.e.,(U − c)( d2dz2− α2)φˆ− ∂2U∂z2φˆ− N2U − cφˆ = 0. (1.21)This is the original TG equation.We non-dimensionalize velocity by a velocity scale q, density by ρ1, and length by l,then the viscous diffusive TG equations become(U∗ − c∗)( d2dz∗2− α∗2)φˆ∗ − ∂2U∗∂z∗2φˆ∗ = Jρˆ∗ +1Re1iα∗(d2dz∗2− α∗2)2φˆ∗, (1.22)41.2. Theoretical background(U∗ − c∗)ρˆ∗ + ρ∗0g∗N∗2φˆ∗ =1RePr1iα∗(d2dz∗2 − α∗2)ρˆ∗, , (1.23)where variables with ∗ are dimensionless variables, Re = ql/ν, Pr = ν/κ, and J =ρ1gl/ρ0q2. J is the bulk Richardson number and measures the importance of stratifica-tion to the shear.The growth rate and phase speed of each mode can be predicted by the TG equationsonce background velocity and density profiles are known. For this thesis, hyperbolic tangentvelocity and density profiles are used as the initial background flow as in Hazel (1972), i.e.,U =∆U2tanh(2zh0)and ρ¯ = −∆ρ2tanh(2zδ0), (1.24)where h0 is the thickness of the shear and δ0 is the thickness of the density stratification.In this thesis, we only consider the case where δ0 = h0. We use ∆U , ∆ρ, and h0 as velocity,density, and length scales to non-dimensionalize field variables, thereforeRe =∆Uh0νand J =∆ρgh0ρ0∆U2. (1.25)We fix Reynolds number Re at 1200, bulk Richardson number J at 0.07, and vary Prvarying from 1 to 64. This combination of parameters allows KH instabilities and vortexpairing. The growth rate as a function of wavenumber is almost independent of Pr andit is shown in figure 1.1 for Re = 1200, Pr = 1, and J = 0.07. When αh0/2 is between0 and 1, the perturbation is unstable. The most unstable mode occurs at αh0/2 = 0.435with growth rate 0.14. The phase speed for this wavenumber is zero and therefore the mostunstable mode is stationary.Although linear modal stability theory gives the correct prediction of the most unstablemode’s wavenumber, it only considers the perturbations where time t and vertical coordinatez can be separated. In general the stability problem is an initial value problem and cannotbe reduced to an eigenvalue problem. Also, the growth rate is a function of time. Whetherthe most unstable mode predicted by modal theory is the most unstable one for all possibleperturbations needs to be answered by non-modal stability theory.1.2.2 Linear non-modal stability theoryThe drawbacks of linear modal stability theory are it only considers a special class ofperturbation and does not explain the physical mechanism of instabilities. Farrell & Ioannou(1996) and Schmid & Henningson (2012) established a generalized stability theory for non-modal perturbations. Constantinoul & Ioannou (2011) applied the theory to show there ispotential for large transient growth even if the background flow is stable for normal mode51.3. Nonlinear evolution of KH instabilities0 0.2 0.4 0.6 0.8 1αh0200.0250.050.0750.10.1250.15αcih0∆UFigure 1.1: Non-dimensionalized growth rate αcih0/∆U as a function of half non-dimensionalized wavenumber αh0/2. The two vertical dashed lines correspond to the mostunstable mode and its subharmonic mode.perturbations. On the other hand, Holmboe (1962) first proposed the resonant interactionbetween progressive interfacial waves at discontinuities of vorticity profile or density profileleads to exponentially growing instaiblities in shear flows. Acoording to his theory, KHinstability is caused by the resonant interaction between two stable vorticity waves located attwo discontinuity points. Holmboe instability is caused by the resonant interaction betweena vorticity wave and a gravity wave. Holmboe’s theory has been extended by Baines (1994),Caulfield (1994), and Carpenter et al. (2011). Carpenter et al. (2011) showed that resonantinteraction occurs when stable progagating waves “phase lock”. Guha & Lawrence (2014)showed that two stable waves phase lock and reach the growth rate predicted by linearmodal stability theory under certain conditions in idealized piecewise velocity and densitybackground flows. It has also been shown that the transient growth rate can be eithersmaller or larger than the linear modal growth rate.1.3 Nonlinear evolution of KH instabilitiesLinear theory is successful when perturbation is small. However, it soon breaks up whennon-linear effects become non-negligible. To study the non-linear evolution of KH insta-bilities, researchers have carried out laboratory experiments (e.g. Thorpe, 1971; Winant &Browand, 1974; Koop & Brownad, 1979; Ho & Huang, 1982) and numerical simulations(e.g. Corcos & Sherman, 1976; Moser & Rogers, 1993; Caulfield & Peltier, 2000; Smyth61.3. Nonlinear evolution of KH instabilitieset al., 2001; Mashayek & Peltier, 2012; Salehiour & Peltier, 2015b) in both homogeneousand stratified shear flows. In this section, we give a brief introduction to the non-linearevolution of KH instabilities in stratified fluids.Figure 1.2: Snapshots of density field at non-dimensional time (a) t∗ = 75, (b) t∗ = 107,(c) t∗ = 200, (d) t∗ = 400 for a simulation with Re = 1200, Pr = 16, and J = 0.07.KH instabilities typically evolve through four stages. During the first stage, the primaryinstability, i.e. KH instability, grows and saturates. It has been shown that this stage isalmost two-dimensional (Caulfield & Peltier, 2000). The kinetic energy of the perturbationquickly increases as mean kinetic energy transfers energy to it through shear production.The potential energy of the system also increases through buoyancy flux. The density fieldduring this stage is shown in figure 1.2 (a).During the second stage, secondary instabilities grow and lead to the transition to tur-bulence. Different secondary instabilities may appear or dominate depending on Reynoldsnumber Re, Prandtl number Pr, or bulk Richardson number J . Most secondary instabilitiesare three-dimensional, e.g. the convective core instability (e.g. Caulfield & Peltier, 2000;Klaassen & Peltier, 1985) and braid instabilities (Mashayek & Peltier, 2012). The convectivecore instability occurs for the parameters chosen in this thesis. Besides three-dimensionalsecondary instabilities, there is a widely known two-dimensional secondary instability, i.e.vortex pairing. Whether vortex pairing occurs depends on the competition between this71.4. Mixing efficiencysecondary instability and three-dimensional secondary instabilities. Vortex pairing has beenobserved in both experiments (see Thorpe, 1971; Winant & Browand, 1974; Koop & Brow-nad, 1979; Ho & Huang, 1982) and three-dimensional simulations (see Moser & Rogers,1993; Smyth & Moum, 2000; Smyth et al., 2001; Rahmani et al., 2014). Figure 1.2 (b)shows the density field of a simulation during this stage when vortex pairing occurs. Ifthree-dimensional secondary instabilities have a higher growth rate than vortex pairing anddestroy the KH billows, vortex pairing does not occur. However, if vortex pairing growsmuch quicker than other three-dimensional secondary instabilities, two KH billows mergeand Reynolds number is increased due to the increasing length scale of the large scale motion(Moser & Rogers, 1993). During this stage, three-dimensionality becomes observable.During the third stage, the flow becomes highly turbulent and KH billows no longer exist.A snapshot of the density field during this stage is shown in figure 1.2 (c). Most of mixingoccurs during this stage (Caulfield & Peltier, 2000; Rahmani et al., 2014) because smallscale motions are prevalent in the middle of the shear layer. As a result, potential energyincreases significantly. However, turbulence cannot be sustained and gradually decays.When the turbulence intensity, measured by buoyancy Reynolds number Reb, becomessmall, the active turbulent stage ends.During the fourth stage, three-dimensional motions are weak and the flow graduallyreturns to a laminar state (see figure 1.2 (d)). Potential energy and mean kinetic energybecome approximately constant. Both the thickness of velocity and density profiles arehigher compared to the initial state. Also, the flow becomes more stable than the initialstate.1.4 Mixing efficiencyThe quantification of turbulent diffusivity, κρ, plays an important role in turbulent modelingof stratified flows. Osborn (1980) derived a formula relating turbulent diffusivity to dissipa-tion rate of turbulent kinetic energy by assuming homogeneous and stationary turbulence,i.e.,κρ = Γε′N2, Γ =Rf1−Rf (1.26)where Rf is flux Richardson number (a measure of mixing efficiency) and defined as theratio of buoyancy flux b to shear production P (Kundu & Cohen, 2008). Γ is the fluxcoefficient and is usually assumed to be 0.2 in physical oceanography. However, turbulencein nature is always inhomogeneous and non-stationary (Ivey et al., 2008). Also, studies havefound that Rf depends on stratification (Rehmann & Koseff, 2004), turbulence intensity(Gargett, 1988; Caulfield & Peltier, 2000; Barry et al., 2001; Salehiour & Peltier, 2015b;81.5. ObjectivesSalehiour et al., 2016), and Prandtl number (Barry et al., 2001; Smyth et al., 2001; Rahmaniet al., 2016; Salehiour & Peltier, 2015b).Shih et al. (2005) extended the the definition of flux Richardson number to non-stationaryturbulence by defining Rf as,Rf =bb+ ε′, (1.27)and found the relation between flux coefficient Γ and buoyancy Reynolds number Reb inuniformly stratified homogeneous turbulent flows. However as Smyth & Moum (2000)pointed out, mixing efficiency in uniformly stratified turbulence is different from that innon-uniformly stratified turbulence and Venaille et al. (2017) verified mixing efficiency doesrely on the stratification profile. Also, the definition in Osborn (1980) and Shih et al.(2005) does not distinguish the irreversible mixing from the reversible mixing which canbe returned to kinetic energy (Winters et al., 1995). To avoid this problem, Salehiour &Peltier (2015a) defined turbulent diffusivity using the irreversible flux defined in Winterset al. (1996) and background density profile. By doing this, Salehiour & Peltier (2015b)connects turbulent diffusivity to irreversible mixing efficiency commonly used in numericalstudies. Moreover, Salehiour et al. (2016) created a contour plot of the irreversible mixingefficiency as a function of buoyancy Reynolds number and stratification. However, thedifficulty of this new diffusivity is that it is not easy to compute. Also, the effect of Prandtlnumber has not been included.1.5 ObjectivesThe objective of this thesis is to investigate the effects of initial conditions and Prandtlnumber on vortex pairing and mixing in KH instabilities. Re = 1200 and J = 0.07 arechosen to allow vortex pairing and the length of the computation domain is two wavelengthsof the KH instabilities. Chapter 2 focuses on the effect of initial condition on vortex pairingand mixing. Random perturbations are used to simulate unforced KH instabilities. Vortexpairing and mixing initialized by the same structure and amplitude random perturbationswith different phase of the subharmonic mode are compared. A few additional simulationsare initialized by eigenfunctions to viscous diffusive TG equations as a comparison withrandom perturbed simulations. Chapter 3 investigates the effect of Prandtl number onvortex pairing and mixing properties.9Chapter 2Sensitivity of vortex pairing andmixing to initial conditions2.1 IntroductionFluids are often stratified in nature due to temperature or salinity, or both. The existenceof shear may give rise to instabilities in these otherwise stably stratified flows. Kelvin-Helmholtz (KH) instabilities, also called Rayleigh instabilities in homogeneous fluids, are oneof the most widely known shear instabilities. KH instabilities have been studied extensivelyin both homogeneous and stratified fluids using laboratory experiments (e.g. Thorpe, 1973;Browand & Winant, 1973; Winant & Browand, 1974), field observations (Seim & Gregg,1994; Moum & Farmer, 2003; Geyer et al., 2010), and numerical simulations (e.g. Patnaiket al., 1976; Klaassen & Peltier, 1985; Caulfield & Peltier, 2000; Mashayek & Peltier, 2012,2013; Rahmani et al., 2014). They are characterized by stationary two-dimensional periodicelliptic vortices called KH billows, which are connected by thin tilted braids of high strainrate (Corcos & Sherman, 1976).When KH instabilities reach maximum amplitude (saturate), they are susceptible toseveral secondary instabilities, e.g. vortex pairing (Browand & Winant, 1973; Winant &Browand, 1974; Koop & Brownad, 1979; Ho & Huang, 1982), convective core instabilitydue to the overturn of fluid caused by the roll-up (e.g. Klaassen & Peltier, 1985; Caulfield &Peltier, 2000), and instabilities that are located in braid regions and extract energy from themean shear or strain (see Mashayek & Peltier, 2012). Which secondary instabilities exist ordominate depends on non-dimensional parameters governing the flows, i.e. Reynolds num-ber, Richardson number, and Prandtl number. In low to intermediate Reynolds numberflows, vortex pairing is a dominant two-dimensional secondary instability. It can increasethe vertical scale of motion and thickness of the shear layer (Corcos & Sherman, 1984;Smyth & Peltier, 1993). As a result, the effective Reynolds number is also increased. Sincethe amount of mixing and mixing efficiency are higher for higher Reynolds numbers in themixing transition regime (Rahmani et al., 2014), vortex pairing is an efficient way to en-hance mixing and mixing efficiency. The dominant three-dimensional secondary instabiltyin this Reynolds number regime is the convective core instability (Klaassen & Peltier, 1985;102.1. IntroductionCaulfield & Peltier, 2000). Caulfield & Peltier (2000) show that the growth rate of the con-vective core instability mainly comes from the mean shear, while the two-dimensional KHinstability acts as a catalyst in the sense that it provides the flow on which the secondaryinstability grows. This competition of vortex pairing and three-dimensional secondary in-stabilities determines whether vortex pairing occurs or not. The competition is dependenton the initial non-dimensional parameters, and also on the details of the initial pertur-bations (Caulfield & Peltier, 2000; Metcalfe et al., 1987), e.g. the amplitudes of KH, thesubharmonic components, and three-dimensional motions.Some researchers have studied the dependence of secondary instabilities on initial condi-tions in shear layers without density stratification, for example Patnaik et al. (1976), Ho &Huang (1982), Ho & Huerre (1984), and Metcalfe et al. (1987). Patnaik et al. (1976) showthat shredding replaces pairing when the KH and the subharmonic modes are completelyout of phase. One vortex is strengthened and the other is weakened in that case. However,shredding is seldom observed in experiments due to the existence of ambient noise otherthan pure eigenfunctions of the Orr-Sommerfeld equation. Ho & Huang (1982) study thespreading rate of a spatially varied shear layer under different forcing. They show that per-turbing the flow with different frequencies can change the the number of vortices mergingtogether and the spreading rate of the shear layer. They also show that the subharmonicmode of the most unstable mode is necessary for pairing, otherwise pairing is significantlydelayed. Metcalfe et al. (1987) study the effect of the most unstable mode, the subharmonicmode, and one three-dimensional mode in initial perturbations on primary and secondaryinstabilities of homogeneous shear flows using numerical simulations. Their results showthat vortex pairing can suppress the modal growth rate of a three-dimensional mode whenthe subharmonic mode reaches finite amplitude and the three-dimensional mode is small.However, one should be careful that this may be only valid for flows initialized by eigen-functions of certain amplitude.Numerical investigations of shear instabilities in stratified flows have also found thatvortex pairing depends on initial conditions, e.g. Klaassen & Peltier (1989) and Smyth& Peltier (1993). Klaassen & Peltier (1989) obtain the amplitude ratios of the first threeharmonics (with wavenumber 12αkh, αkh, and32αkh where αkh is the wavenumber of the mostunstable mode predicted by the Taylor-Goldstein equation) in a two-wavelength domainfrom a numerical simulation perturbed by white noise. They employ these three modesto study the impact of the phase of the subharmonic and third harmonic mode on vortexpairing for two-dimensional flows. Their results demonstrate that pairing is delayed andthe growth rate of the subharmonic mode is decreased if the subharmonic and the thirdmodes are out of phase relative to KH instabilities. However, the time of vortex pairingmay be sensitive to the 32αkh mode if the subharmonic mode is out of phase. Smyth &112.1. IntroductionPeltier (1993) show similar results in two-dimensional simulations.The effects of initial perturbations on vortex paring and mixing in stratified shear flowswarrants further investigation. In the present study, we investigate how vortex pairingis influenced by the phase difference of the subharmonic component relative to KH com-ponent and extend the study of pairing to three-dimensional flows, by carrying out eightnumerical simulations with a variety of initial perturbations. We also investigate whetherinitial non-modal growth of these two modes, or the existence of other modes in the initialperturbations, affects vortex pairing and mixing properties.There are four types of perturbations commonly used to excite primary and secondaryinstabilities in numerical simulations of stratified shear flows. The first type is pure two-dimensional eigenfunctions (e.g. Patnaik et al., 1976; Klaassen & Peltier, 1989). This typeis usually used for two-dimensional simulations. The second type of perturbation is a com-bination of the eigenfunction of the most unstable mode and a smaller amplitude randomnoise to excite secondary instabilities (e.g. Smyth & Peltier, 1993; Caulfield & Peltier, 2000;Mashayek & Peltier, 2012; Rahmani et al., 2014; Salehiour & Peltier, 2015b; Rahmani et al.,2016). An implicit assumption of this type of perturbation is that the KH mode dominatesinitially and three-dimensional perturbations grow linearly while KH instabilities saturate.The third type of perturbations is a combination of sinusoidal waves of wavenumbers αkhand 12αkh with vertical exponential decay, and random noise (e.g. Smyth & Moum, 2000;Smyth & Winters, 2003; Alexakis, 2009). The fourth type is pure random noise (e.g. Car-penter et al., 2010; Mashayek & Peltier, 2013). The present study was motivated by theobservation that two sets of small random perturbations, that are generated by the samerandom number generator with different seeds, can trigger very different flows, i.e. vortexpairing occurs in one simulation but not in the other simulation. To verify the conjecturethat if the subharmonic mode is out of phase relative to KH mode vortex pairing may notoccur, we conduct numerical simulations where the background flow is perturbed by thesame pure random number perturbations with differing phase of the subharmonic mode.We also compare flows perturbed by pure random perturbations with flows perturbed byeigenfunctions.The numerical methods and diagnostic tools are described in section 2.2. A simplifiedpairing mechanism is described in section 2.3. Section 2.4 describes the process of vortexpairing and the growth rate of the subharmonic mode in two-dimensional simulations. Insection 2.5, three-dimensional results are compared with two-dimensional results and mixingis compared in different simulations.122.1. IntroductionFigure 2.1: Vorticity snapshots of two-dimensional simulations R02D (left panel, phase ofthe subharmonic mode is θMsub ≈ 0) and Rpi2 2D (right panel, phase of the subharmonic modeis θMsub ≈ −pi2). The snapshots labeling with 3D are vorticity field of simulations R03D (left)and Rpi2 3D (right) at y =Ly2. Pairing is delayed in two-dimensional simulation Rpi2 2D butcompletely eliminated in the three-dimensional simulation Rpi2 3D. Black circles are fluidparticles located at vortex centers at t = 30.132.2. Methodology2.2 Methodology2.2.1 Mathematical modelHyperbolic tangent functions are used for the background velocity and density profiles, asfirst introduced by Hazel (1972),ρ¯ = −∆ρ2tanh(2zδ0)and U =∆U2tanh(2zh0), (2.1a, b)where δ0 is the thickness of the density interface and h0 is the thickness of the velocityinterface. Four non-dimensional parameters characterize the flows, i.e. bulk Richardsonnumber J , Reynolds number Re, Prandtl number Pr, and the scale ratio R which aredefined asJ =∆ρgh0ρ0(∆U)2, Re =∆Uh0ν, Pr =νκ, R =h0δ0. (2.2)In this study, J = 0.07, Re = 1200, Pr = 16, R = 1. The flow is susceptible to Kelvin-Helmholtz instabilities for this combination of J and R (see Smyth & Winters, 2003, for areview of instability types ).We assume the fluid is incompressible and apply the Boussinessq approximation, so thegoverning equations for the system are∇ · u = 0, (2.3)DuDt= − 1ρ0∇p− ρρ0gkˆ + ν∇2u, (2.4)DρDt= κ∇2ρ, (2.5)where κ is molecular diffusivity, ν is kinetic viscosity, ρ0 is a reference density, and kˆ is theunit vertical vector.2.2.2 Direct Numerical SimulationsThe governing equation (2.3), (2.4), and (2.5) are solved by a pseudo-spectral code developedby Winters et al. (2004) and later improved by Smyth et al. (2005). The code employs thethird order Adams-Bashforth time stepping scheme. Boundary conditions are horizontallyperiodic and vertically free slip and no flux for our simulations.The resolution of DNS is typically determined by the Kolmogrov scale, Lk = (ν3/ε′)1/4,in homogeneous fluids where ε′ is viscous dissipation rate of turbulent kinetic energy. Moin& Mahesh (1998) suggest that the grid spacing in DNS should be O(Lk). In stratified flowswith Pr > 1, the smallest scale that needs to be resolved is the Batchelor scale (Batchelor,142.2. MethodologyRey Pr J Lx/h0 Ly/h0 Lz/h0 Nx Ny Nz1200 16 0.07 14.43 7.22 15 320 160 320Table 2.1: Numerical parameters for all the simulations. The number of grid points is forthe velocity field and is half that of the density field.1959), LB = Lk/√Pr. In our simulations, ∆z/LB is always less than 4.0 and ∆z/LK isalways less than 2.0 (∆z of the density field is half of the velocity field.). The dissipationrate ε′ used to calculate LK is averaged within the initial shear layer −h0/2 < z < h0/2where turbulence is the most energetic. The domain length Lx is set to two wavelengths ofthe most unstable mode to allow vortex pairing. The spanwise width of the domain Ly forthe three-dimensional simulations is one wavelength of the primary KH instability, whichis at least six wavelengths of most unstable spanwise mode (for the most unstable spanwisewavenumber see Klaassen & Peltier, 1989). The height is 15h0 which is sufficient to removethe effects of the horizontal boundaries. The numerical details are summarized in table 2.1.We ran two sets of simulations to study the effect of initial perturbations on vortexpairing and mixing. One set is perturbed by random perturbations, three of them aretwo-dimensional and three of them are three-dimensional simulations. The other set isperturbed by eigenfunctions of KH and the subharmonic modes with the same amount ofkinetic energy (defined in equation (2.13)) as in random perturbation simulations. Theeigenfunction simulations are performed in two dimensions only. All simulations are listedin table 2.2.For random perturbation simulations, inherently three dimensional random pertur-bations of u′ and w′ are added to the background flow to excite instabilities in three-dimensional simulations. They are given by the following equations,u′ = cru(x, y, z)∆U2(1−∣∣∣∣tanh(2zh0)∣∣∣∣) , (2.6)w′ = crw(x, y, z)∆U2(1−∣∣∣∣tanh(2zh0)∣∣∣∣) , (2.7)where ru and rw are random numbers between -1 and 1, and c sets the maximum amplitudeof perturbations. In the present study, c = 0.1, which is the same magnitude as in thesimulations of Smyth & Winters (2003) and Carpenter et al. (2010), and small enough forperturbations to grow linearly initially. To study the effects of the phase of the subharmoniccomponent, these simulations are initialized by the same perturbations except that the phaseof the subharmonic mode is changed. The initial conditions in two-dimensional simulations152.2. Methodologyare spanwise averaged values of those in corresponding three-dimensional simulations.The phase of each wavenumber component is defined in terms of two-dimensional verticalvelocity w2d , i.e.,A sin(k2piLxx+ θk)= <{eikαxwˆ2d,k(z =Lz2)}, (2.8)where A is a real number, wˆ2d,k is the kth Fourier component of w2d. Note that k = 2corresponds to the KH component and k = 1 corresponds to the subharmonic component,i.e. a wavelength equal to the domain length, Lx. Initially each Fourier component ex-periences non-modal growth and the phase of every component defined in equation (2.8)changes. When the subharmonic component becomes approximately modal, i.e. identicalto the eigenfunction of the Taylor-Goldstein equation (this occurs around t = 24 for threerandom perturbation simulations), the phase becomes almost constant for some time untilnon-linear effects become important. We define θsub as the phase of the subharmonic com-ponent relative to the KH component and θMsub as the phase of the subharmonic componentwhen it becomes modal. We use three different phases in the random perturbation simu-lations. They are designated R followed by the approximate modal phase θMsub, and 2D or3D depending on the number of dimensions used in the simulation. For example, θMsub isapproximately 0, −pi4 , and −pi2 respectively in three-dimensional simulations R03D, Rpi4 3D,and Rpi2 3D (exact phase values are listed in table 2.2).Simulations perturbed by eigenfunctions are named by the same rule, but the firstletter is E, indicating that they are perturbed by eigenfunctions. For these eigenfunctionperturbed simulations, the phase of the subharmonic mode does not change initially andθMsub is the initial value. θMsub = −pi2 for E pi2 2D and θMsub = 0 for E02D.2.2.3 Diagnostic toolsHereafter, all results are presented in non-dimensional form. Specifically, velocity is non-dimensionalized by ∆U , density is non-dimensionalized by ∆ρ, length by h0, time by h0/∆U ,and pressure by ρ0(∆U)2. Following Caulfield & Peltier (2000), the velocity is decomposedinto three parts, i.e.,u = 〈u〉xy, u2d = 〈u〉y − 〈u〉xy, u3d = u− u− u2d. (2.9a–c)Given these definitions, the total kinetic energy K is defined asK =〈u · u〉xyz2, (2.10)162.2. MethodologyRandom perturbations EigenfunctionsRun R02D Rpi4 2D Rpi2 2D R03D Rpi4 3D Rpi2 3D E02D Epi2 2DtM 24 24 24 24 24 24 0 0θMsub/pi 0.04 - 0.21 -0.46 0.05 -0.21 -0.48 0 0.5tkh 80 82 85 80 81 84 81 87tsub 106 110 146 107 113 129 104 249tp 108 112 148 109 113 – 107 249t3d – – – 143 146 123 – –tf – – – 221 214 187 – –M(t = 400)/K0 – – – 0.072 0.068 0.036 – –Ec(t = tf ) – – – 0.229 0.223 0.198 – –Table 2.2: The simulations beginning with E are perturbed by eigenfunctions, while thosebeginning with R are perturbed by the same random perturbations except different phasesof the subharmonic mode. The kinetic energy of KH and the subharmonic modes in eigen-function simulations are the same as in random perturbation simulations. tkh, tsub, and tpare defined in section 2.4. t3d, tf and Ec are defined in section 2.5.and can be partitioned into three parts K, K2d, K3d, i.e.,K = K +K2d +K3d, (2.11)whereK =〈u · u〉z2, K2d =〈u2d · u2d〉xz2, K3d =〈u3d · u3d〉xyz2, (2.12a–c)and the subscripts indicate averaging over that direction.Fourier transforms are applied to u2d and w2d in order to identify the contribution of eachwavenumber component to the total K2d, so that the kinetic energy of the kth componentisK2d,k = 〈|uˆ2d,k|2 + |wˆ2d,k|2〉xz, k > 1 (2.13)where uˆ2d,k and wˆ2d,k are the Fourier components of u2d and w2d of wavenumber 2pik/Lx.Hence,K2d =Nx2∑k=1K2d,k. (2.14)Note that k = 1 corresponds to the subharmonic component and we denote it as Ksub.172.3. Pairing mechanismαkhx/2w'(a)0 pi2 pi3pi2 2piαkhx/2w'(b)0 pi2 pi3pi2 2piKHSubharmonicFigure 2.2: Conceptual drawing of vortex pairing demonstrating the effect of phase of thesubharmonic mode. (a) θsub = 0 (b) θsub = −pi2 . The blue solid line is KH mode and the reddash-dotted line is the subharmonic mode. Blue circles denote the location of KH vortexcenters and red squares are vortex centers associated with the subharmonic mode. The twoarrows show the directions of the mean flow.k = 2 corresponds to KH instability and we denote it as Kkh. These two components of thekinetic energy characterize the kinetic energy of the primary and subharmonic components.2.3 Pairing mechanismBefore discussing results, we first illustrate the vortex pairing mechanism in figure 2.2 usingthe vertical velocity of the KH instability (blue) and subharmonic mode (red). This figureshows two cases with the phase equal to 0 and −pi2 . The circles represent core centersassociated with KH instabilities and the squares are potential vortex cores associated withthe subharmonic mode. The arrows show the direction of the mean flow. In figure 2.2 (a),the phase of the subharmonic mode θsub = 0 and its associated core center is at the centreof two KH billows. The subharmonic mode displaces the left KH billow upward and theright KH billow downward. The two KH billows are then advected toward each other by themean flow, cross each other, and merge into one larger billow. This is the optimal phase forpairing. In figure 2 (b), the phase of the subharmonic mode is θsub = −pi2 and two KH corecenters are at the nodes of the subharmonic mode. The phase of the subharmonic modeneeds to shift towards 0 or pi to impose a vertical velocity at the two KH core centers forvortex pairing to occur. As a result, pairing is delayed until this shift occurs.182.4. Two-dimensional results2.4 Two-dimensional results2.4.1 Degree of modalityIn the randomly perturbed simulations, initially the subharmonic component experiencesnon-modal growth. We use the Hermitian angle (Scharnhorst, 2001) between the subhar-monic component wˆ2d,sub and the initial eigenfunction of the subharmonic mode wˆeig,sub tomeasure the degree of modality, i.e.,r(t)=∣∣∣∣∣Nz∑k=1wˆ∗2d,subwˆeig,sub∣∣∣∣∣√∑k|wˆ2d,sub|2∑k|wˆeig,sub|2(2.15)where ∗ denotes the complex conjugate and || denotes the amplitude of a complex number.This Hermitian angle r(t) is the ratio between the length of the orthogonal projectionof the subharmonic component onto the the eigenfunction to the length of itself. It isalways between 0 and 1, and r = 1 only when the subharmonic component of the randomperturbation is identical to the eigenfunction.Figure 2.3 shows the evolution of r for the three random perturbation simulations R02D,Rpi4 2D, and Rpi2 2D. Initially, r is close to zero because the subharmonic component of theinitial random perturbations is significantly different from the eigenfunction. However, asthe subharmonic component evolves to the eigenfunction, r increases to 1. We define thetime required for the subharmonic component to become modal, tM , as the time when rfirst exceeds 0.99, which is when t = 24 for these three simulations. Note that before t = 24,the three simulations are evolving linearly and therefore appear identical in figure 2.3.2.4.2 Phase evolutionThe phases for the five two-dimensional simulations perturbed by either random pertur-bations or the eigenfunctions are plotted in figure 2.4 (a). Initially, in the eigenfunctionperturbation simulations, the phase does not change. In the three random perturbationsimulations, before t = 25, the phases change because of non-modal growth. Since thesubharmonic component is growing linearly, the phase differences between any two of thethree random perturbation simulations are constant. The phase is not meaningful duringthis time as the subharmonic component has not become modal and the phase is the phaseat the vertical center (z = Lz2 ) of the domain (see equation 2.8). The phase before tM istherefore shown as thin lines. Between t = 25 and t = 50, the phases stay almost constant.During this period, the phases are approximately 0,−pi4 ,−pi2 for R02D, Rpi4 2D, and Rpi2 2D192.4. Two-dimensional results0 25 50 75 100 125t00.20.40.60.81rR02DRpi42DRpi22DFigure 2.3: Evolution of r for random perturbation simulations R02D, Rpi4 2D, and Rpi2 2D.(see table 2.2). After t = 50, the phases of the subharmonic mode in simulations Rpi4 2D andRpi2 2D shift toward 0 (figure 2.4 (a)). For the eigenfunction simulation Epi2 2D, the phasebegins to shift after t = 200 and is not shown in this figure.2.4.3 Lagrangian trajectoriesThe Lagrangian trajectories of two fluid particles located at the density inflection points att = 30 for R02D, Rpi4 2D, and Rpi2 2D are tracked as representatives of KH billow centres.Time t = 30 is the earliest time that the vortex centres are identifiable. The trajectories oftwo fluid particles initially located at Lx/4, 3Lx/4 are also tracked for the two eigenfunctionsimulations, E02D and E pi2 2D. The Lagrangian trajectories show good agreement with thepaths of inflection points on the ρ0 isopycnal until the two KH billows get close to eachother. Figure 2.1 shows the evolution of the vorticity field with the two fluid particles shownas black circles. These two fluid particles approximately represent the vortex centres untilsmall scale motions prevail, e.g. at t = 146 for simulation R02D.Figure 2.4 (b) shows the horizontal coordinates of the two fluid particles. Initiallythe horizontal distance between the two fluid particles for the eigenfunction simulationsis almost constant since the KH instabilities are stationary. After approximately t = 75,the two fluid particles quickly converge for the optimal phase simulation E02D. When thehorizontal distance becomes zero, one fluid particle (billow) is on top of the other. Thistime is defined as the time of pairing, tp, and listed in table 2.2 (for three-dimensionalsimulations, tp is obtained by averaging two sets of trajectories each composed of 21 fluidparticles spread over the spanwise direction). The two vortices become closer and mergein simulation E pi2 2D, but only at t = 249 (see table 2.2) because the phase needs to shift202.4. Two-dimensional results0 25 50 75 100 125 150t-1-0.500.51θsub/pi0 25 50 75 100 125 150t0.250.50.7511.251.51.75x/λ(b)R02DE02DR pi2 2DE pi2 2DR pi4 2D0 25 50 75 100 125 150t-0.0500.050.10.151 2d(lnK2d,sub)/dt(c)→tkhTGFigure 2.4: (a) Evolution of the phase of the subharmonic mode. The phase before tM isshown as thin lines. The phase shift toward 0 in simulation E pi2 2D begins after t = 150 andis not shown in this figure. (b) x coordinates of two fluid particles located at two vortexcenters at t = 30 for random perturbation simulations or located at Lx/4 and 3Lx/4 initiallyfor eigenfunction perturbation simulations. Pairing occurs at t = 249 for simulation E pi2 2Dand is not shown in this figure (see table 2.2). (c) Growth rate of the subharmonic mode.The vertical dashed line indicates the saturation time of KH instabilities in simulationRpi4 2D.212.4. Two-dimensional resultstoward the optimal value and pairing is delayed. Unlike the results of the eigenfunctionsimulations, in the random perturbation simulations there is an oscillation of the fluidparticles. For these three simulations, the distance between two fluid particles is alwaysthe smallest for R02D and largest for Rpi2 2D. Also, pairing occurs first in R02D and lastin Rpi2 2D, so tR02Dp < tRpi42Dp < tRpi22Dp . Relating figure 2.4 (b) with the phase evolutionin figure 2.4 (a), we find that two fluid particles begin to move together when the phaseapproximately stops changing. Therefore, pairing occurs earliest if the subharmonic is inphase and latest if it is out of phase. tRpi42Dp is just slightly larger than tR02Dp , which suggeststhat if the subharmonic mode is not close to being out of phase, the difference caused bythe phase is small. Also, the time of pairing for simulation R02D is close to E02D, butpairing happens much earlier for simulation Rpi2 2D than Epi2 2D. Hence, the existence ofthe other components in initial perturbations does not affect the time of pairing if thephase is optimal. However, if the subharmonic mode is close to ±pi2 , the time of pairing issensitive to initial perturbations, e.g. the existence of the other components or whether thesubharmonic component is exactly the eigenfunction.2.4.4 Growth rateIn this section, we discuss temporal changes in the kinetic energy of the subharmonic com-ponent. Figure 2.4 (c) shows the growth rate for three random perturbation simulations andtwo eigenfunction simulations. tkh is defined as the time when Kkh (the 2D kinetic energyof the KH mode) reaches its first maximum, also called saturation time. This occurs atapproximately t = 81 with slight variations for all simulations (see table 2.2). The verticaldashed line in this figure indicates tkh for Rpi4 2D.Initially, the growth rates of the subharmonic mode for the two eigenfunction simulations(E02D and E pi2 2D) are the same and decline as the shear layer diffuses. The estimatedgrowth rate using the viscous diffusive Taylor-Goldstein equation and the mean flow hasthe same decreasing trend as the growth rates based on K2d,sub (the 2D kinetic energy of thesubharmonic mode). For the random perturbation simulations (R02D, Rpi4 2D, and Rpi2 2D),the growth rate during the non-modal stage of growth can be either smaller or larger thanthe modal growth rate as shown by Guha & Lawrence (2014). For all five simulations, thegrowth rate is independent of the phase in the initial linear stage of growth, i.e. beforearound t = 45.After around t = 45, non-linear effects and the phase become important. After t = 45and before tkh, the growth rates in the late pairing simulations (Epi2 2D and Rpi2 2D) signif-icantly drop compared to the other simulations, which is consistent with the simulationsof Klaassen & Peltier (1989). The growth rate in the Rpi4 2D simulation stays closer to theearlier pairing simulations R02D and E02D. Therefore the subharmonic mode is only sup-222.4. Two-dimensional resultspressed if the phase is close to ±pi2 . The growth rates in simulations R02D and E02D arealmost the same and the growth rates in simulations Rpi2 2D and Epi2 2D are almost the same.This suggests that if the phase and amplitude of the subharmonic and KH are the same foran eigenfunction perturbation and a random perturbation, the growth of the subharmonicmode is the same before saturation of KH instabilities. In other words, the existence ofthe other components in initial perturbations and initial non-modal growth have negligibleeffects on the growth rate during this non-linear growth stage.At saturation of KH, in R02D, Rpi4 2D, and E02D, the phase is close to optimal. Thegrowth rates then quickly decrease to zero. In these three simulations, the first zero crossingof the growth rate is close to tp and denotes the saturation of the subharmonic mode, i.e.the global maximum of Ksub.In Rpi2 2D, after tkh the growth rate begins to increase. This coincides with the phaseshift toward the optimal value (see figure 2.4 (a)). In this simulation, the saturation ofthe subharmonic mode occurs much later at t = 146. In E pi2 2D, the phase remains at−pi2 and the growth rate continues to decrease. Unlike the other simulations, the growthrate crosses zero before saturation of the subharmonic mode (at t = 249). Comparisonbetween simulations R02D, Rpi4 2D, and E02D shows that the growth rate is almost thesame even until pairing as long as the phase is not close to −pi2 . However, comparisonbetween simulations Rpi2 2D and Epi2 2D shows that the growth rate is sensitive to the exactform of the subharmonic itself or the existence of the other components in initial conditions.In table 2.2, the saturation times of KH and the subharmonic mode and time of pairing(tkh, tsub, and tp) are summarized. Pairing occurs first in simulation E02D, second inR02D, third in Rpi4 2D, fourth in Rpi2 2D, and last in Epi2 2D. In all simulations tsub is closeto the time of pairing, tp. Ho & Huang (1982) obtain a similar result qualitatively in theirexperiments.2.4.5 Visualization of vertical velocity of the subharmonic componentWe show the evolution of normalized vertical velocity of the subharmonic component, i.e.<{eiαkhx/2wˆ2d,sub/〈|wˆ2d,sub|2〉z} at five discrete times in simulations R02D and Rpi2 2D infigure 2.5 to examine its structure change. At t = 10, the velocity component still showssome randomness because it has not become modal as indicated in the low value of r infigure 2.3. The structure is the same for two simulations except a phase difference, as theflow is evolving linearly. At t = 25, the velocity component looks regular and is almostidentical to the eigenfunction (not shown in this figure). Again, it is consistent with thefact that r is more than 0.99. The asymmetry in the velocity component comes from theinitial randomness. At t = 75 (before tkh), the structure does not change much and stillstays close to the structure at t = 25 in simulation R02D. In simulation Rpi2 2D, however,232.4. Two-dimensional resultsR02Dt =100 5 10-505t =250 5 10-505t =750 5 10-505t =1000 5 10-505t =1250 5 10-505(z−0.5Lz)/h0Rpi22Dt =100 5 10-505t =250 5 10-505t =750 5 10-505t =1000 5 10-505t =1250 5 10x/h0-505Figure 2.5: The normalized vertical velocity of the subharmonic component, i.e.<{eiαkhx/2wˆ2d,sub/〈|wˆ2d,sub|2〉z} in two-dimensional simulation R02D (left panel) and simu-lation Rpi2 2D (right panel).242.5. Three-dimensional resultsthe velocity component clearly changes its structure and this change is caused by non-linearinteraction between different modes. This is also shown by the fact that r hardly changesin simulation R02D, but drops in simulation Rpi2 2D. Smyth & Peltier (1993) show that thegrowth rate of the subharmonic component mainly comes from the mean flow even in thenonlinear growth stage before saturation of KH instabilities. Hence, the structure changeof the subharmonic velocity component in late pairing simulation Rpi2 2D accounts for thedecrease in growth rate shown in figure 2.4 (c), given that the mean flow is almost the samebefore tkh for two simulations. At t = 100, two KH billows start to pair in simulation R02D(see figure 2.4 (b)). The corresponding velocity component becomes slightly deviated fromthe eigenfunction. This deviation is more significant for t = 125.In simulation Rpi2 2D, att = 100 the velocity component is significantly different from that at t = 25 and so is truefor r. Correspondingly, the growth rate is still low. At t = 125, the structure becomessimilar to the eigenfunction again, more accurately the structure at t = 80 in simulationR02D. This explains the high growth rate at this time shown in figure 2.4 (c).2.5 Three-dimensional results2.5.1 Kinetic energyKinetic energy of the subharmonic component in two- and three-dimensional random pertur-bation simulations are compared in figure 2.6 to illustrate the effects of three-dimensionalmotions on vortex pairing. The peak energy of the subharmonic mode is reduced in allthree-dimensional simulations. The kinetic energy of the subharmonic mode, Ksub, is leastaffected by three-dimensional motions in the optimal phase simulation, i.e. 0 phase. Thisis because vortex pairing happens before the three-dimensional motions develop. The satu-ration time of the subharmonic mode for R03D is almost identical to the R02D simulation(see figure 2.6 (a) and table 2.2). The comparison between Rpi4 2D and Rpi4 3D is similar (seefigure 2.6 (b)). Thus when the phase is close to the optimal value, vortex pairing is delayedslightly in three-dimensional simulations.Vortex pairing does not occur in the three-dimensional simulation Rpi2 3D, which is shownin the significant difference between the two lines in figure 2.6 (c). The peak in the three-dimensional kinetic energy in Rpi2 3D is much smaller and occurs earlier (see table 2.2) thanthe peak in the three-dimensional kinetic energy in R03D and Rpi4 3D. During the extratime needed in Rpi2 3D for the phase shift from −pi2 to ∼ 0, the three-dimensional motionsgrow and suppress vortex pairing before it can occur.252.5. Three-dimensional results0 25 50 75 100 125 150 175t00.020.040.060.08(a)Ksub, R02DKsub, R03DK3d, R03D0 25 50 75 100 125 150 175t00.020.040.060.08(b)Ksub, Rpi42DKsub, Rpi43DK3d, Rpi43D0 25 50 75 100 125 150 175t00.020.040.060.08(c)Ksub, Rpi22DKsub, Rpi23DK3d, Rpi23DFigure 2.6: Kinetic energy of the subharmonic mode Ksub and three-dimensional kineticenergy K3d in two-dimensional and three-dimensional simulations.2.5.2 MixingWe follow the framework in Winters et al. (1995) to compare the mixing between simula-tions. The potential energy is then defined as,P = g〈ρz〉v. (2.16)Potential energy P is partitioned into background potential energy Pb and available potentialenergy Pa defined asPb = g〈ρzb〉v, Pa = P − Pb, (2.17a, b)where zb is the location of fluid parcels after being re-arranged into a statically stable state(see Winters et al., 1995; Caulfield & Peltier, 2000). Available potential energy characterizesthe energy that can be exchanged between potential energy and kinetic energy, while theincrease in background potential energy quantifies irreversible mixing in a closed system.The amount of mixing caused by fluid’s motion isM = ∆Pb − φit =∫φMdt, (2.18)where φM is the rate of mixing and φi is the mixing caused by molecular diffusion withoutfluid motion, i.e.,φi =−κg(ρ¯top − ρ¯bottom)Lz. (2.19)262.5. Three-dimensional results0 100 200 300 400t00.020.040.060.08∆P/K0(a)0 100 200 300 400t00.020.040.060.08(∆Pb−φit)/K0(b)R03DR pi4 3DR pi2 3DFigure 2.7: (a) The increase in total potential energy ∆P divided by initial kinetic energyK0, (b) the increase in mixing M divided by initial kinetic energy K0.Cumulative mixing efficiency (Caulfield & Peltier, 2000) is used as a measure of overallmixing properties in this study. It is defined asEc =∫ tft3dφMdt∫ tft3dφMdt+∫ tft3dεdt, (2.20)where tf is defined as the time when buoyancy Reynolds number Reb = ε′/ν〈N2〉z firstdrops below 20 after t3d. This period is chosen as studies indicate that turbulence is activeas long as Reb > 20 (Smyth & Moum, 2000).Figure 2.7 (a) shows the increase in the total potential energy, ∆P , with time. The firstpeak in R03D and Rpi4 3D represents vortex pairing and the first peak in Rpi2 3D representssaturation of KH . Time variation of ∆P is the same for the three simulations before t = 80,indicating that the KH instability is independent of the subharmonic mode. The peak dueto vortex pairing does not exist for Rpi2 3D since pairing never occurs. Overall, the increasein total potential energy with vortex pairing is much higher than without vortex pairingbecause vortex pairing efficiently increases the vertical scale of fluid’s motion. The increasein total potential energy in Rpi4 3D is slightly lower than R03D.Figure 2.7 (b) shows the amount of mixing divided by initial kinetic energy for three-dimensional simulations, M/K0. For all simulations, the amount of mixing is negligiblebefore saturation of KH instabilities as the flow is well-organized and mixing is only slightlygreater than the mixing due to molecular diffusion. After about t = 130, the amount ofmixing significantly increases as small scale motions reach sufficient amplitude. As turbu-lence subsides, the amount of mixing gradually approaches a constant. The difference in272.6. Discussionmixing among the three simulations becomes noticeable after about t = 150. The amountof mixing in simulation R03D with vortex pairing is higher than simulation Rpi2 3D with-out pairing, since vortex pairing increases the amount of stirring. At t = 400 the amountof mixing M/K0 in simulation R03D (0.072) is double that of simulation Rpi2 3D (0.036).Therefore the density interface after re-laminarization is sharper in Rpi2 3D. Mixing in simu-lation Rpi4 3D stays close to and is slightly lower than simulation R02D. Hence, the amountof mixing is significantly reduced if the subharmonic mode is out of phase and vortex pairingdoes not occur. If the phase of subharmonic mode relative to KH mode is not close to ±pi2 ,the effect of phase on mixing is small.Finally the cumulative mixing efficiency defined in equation (2.20) as a measure ofmixing efficiency during the active turbulence stage is compared in table 2.2 among the threedifferent initial conditions to obtain the range of mixing efficiency. We find the maximumand minimum cumulative mixing efficiency are 0.229 (with pairing) and 0.198 (withoutpairing). Eigenfunction simulations with the same phase and amount of three-dimensionalmotions give the same mixing efficiency. The cumulative mixing efficiency for turbulencecaused by shear instabilities is around 0.2 as indicated in many numerical studies (Peltier& Caulfield, 2003). Therefore the strong and weak pairing simulations set the bounds forcumulative mixing efficiency for flows with the same initial non-dimensional parameters andthe same amount of perturbation energy.2.6 DiscussionWe find that vortex pairing is sensitive to initial conditions when the phase of the subhar-monic mode is close to ±pi2 . For simplicity in this discussion, we use tsub to characterize thetime of pairing. In general, tsub is a function of all modes in initial conditions, not only thesubharmonic mode.We consider the sensitivity of tsub to the phase of the subharmonic mode by runningtwo-dimensional simulations perturbed by KH and the subharmonic mode eigenfunctions.Note that tsub is minimal in simulations with θMsub = kpi, and is maximal when θMsub =(k + 12)pi where k is an integer. Also, tsub is an even function of θMsub and symmetric aboutθMsub = (k +12)pi. We calculate tsub for 12 discrete phases in the range of −pi2 to 0 and plotthe results and the symmetric part between 0 and pi2 in figure 2.8. The figure shows thatthe time of pairing is not sensitive to the phase when the phase is close to 0. However, thetime of pairing varies significantly near θMsub = ±pi2 , which indicates that tsub is sensitive tothe phase when the phase is close to ±pi2 . Also, slight deviation from the eigenfunction ofthe subharmonic mode in initial conditions can also change the time of pairing. Hence, thedelay of vortex pairing is sensitive to the functional form of the subharmonic component,282.7. Conclusions100125150175200225250t sub-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5θMsub/piEigenfunctionsRandom perturbationsFigure 2.8: Saturation time of the subharmonic mode, tsub, as a function of the phase θMsubfor two-dimensional simulations perturbed by eigenfunctions and random perturbations.Note tsub is approximately equal to the time of pairing (see table 2.2).phase, and other modes in two-dimensional simulations when the phase is close to ±pi2 .Figure 2.8 also shows the time of pairing tsub for two-dimensional random perturbationsimulations. The results show the same trend as the eigenfunction simulations. The initialnon-modal growth and existence of other modes in the initial conditions cause the slightdifference between random perturbation and eigenfunction results.2.7 ConclusionsTwo-dimensional and three-dimensional DNSs were performed to investigate the effect ofinitial conditions on vortex pairing and mixing that are subject to KH instabilities. Wefocused on the impact of phase of the subharmonic mode on vortex pairing and mixing. Intwo-dimensional simulations, vortex pairing is delayed if there is a phase difference betweenKH and the subharmonic mode. This delay is small when the subharmonic mode is slightlyout of phase. However, when the phase difference becomes close to ±pi2 , the delay increasessignificantly and is also sensitive to the other details of the initial conditions. Compar-ison between eigenfunction simulations and random perturbation simulations shows thatthe response of the flow is similar provided the phase difference is small and not close to±pi2 . If the subharmonic mode is close to being out of phase, small deviations from theeigenfunctions of the KH and subharmonic mode will alter the time of pairing.In three-dimensional simulations, vortex pairing is always suppressed by three-dimensional292.7. Conclusionsmotions and the suppression is greater when the phase difference is larger. Three-dimensionalmotions can even grow to sufficient amplitude and eliminate pairing when the phase dif-ference is sufficiently large. Weaker pairing leads to less mixing and a reduced mixingefficiency. The decrease in mixing and mixing efficiency is greater as the subharmonic modebecomes close to being out of phase, similar to the behavior of the delay of pairing.30Chapter 3The effect of Prandtl number onKH instabilities3.1 IntroductionKH instabilities are an important mechanism for mixing in statically stable stratified fluids,e.g., in both the atmosphere and ocean. In the atmosphere, the air is stratified due totemperature and the Prandtl number Pr is about 0.7. In the ocean, the Prandtl numberis about 7 if the stratifying agent is heat and is about 700 if the stratifying agent is salt.In this chapter, the main objective is to understand the dependence of mixing on Prandtlnumbers.There have been many studies in KH instabilities, but most numerical studies havefocused on flows with Prandtl number close to 1 (e.g. Caulfield & Peltier, 2000; Mashayek& Peltier, 2013) as high Prandtl number flows require a high resolution for the densityfield. In experimental studies, the stratifying agent is usually salt and Prandtl numberis about 700 (e.g. Thorpe, 1973; Winant & Browand, 1974). However, high resolutionmeasurements were lacking in these experimental studies. Moreover, the dependence ofmixing on Prandtl number in KH instabilities in stratified flows is complicated. Klaassen& Peltier (1985), Mashayek & Peltier (2013), Salehiour & Peltier (2015b), and Rahmaniet al. (2016) have shown that turbulence intensity is stronger at higher Prandtl numberand therefore the potential for mixing is higher. However, lower molecular diffusivity candecrease the amount of mixing that occurs. Therefore, the effect of Prandtl number onmixing is still not well understood and needs further investigation.Some researchers have investigated the effect of Prandtl number on KH instabilities,secondary instabilities, and mixing properties, for example Klaassen & Peltier (1985), Smythet al. (2001), Rahmani et al. (2016), Mashayek & Peltier (2012), Salehiour & Peltier (2015b).Klaassen & Peltier (1985) show that Prandtl number can suppress primary KH instabilitiesand increase the growth rate of the secondary convective core instabilities in intermediateReynolds number stratified flows. Mashayek & Peltier (2012) and Salehiour & Peltier(2015b) studied the effect of Prandtl number on several secondary instabilities at highReynolds number. Smyth et al. (2001) show that mixing efficiency for Pr = 7 is significantly313.2. Methodslower than that for Pr = 1 in a two wavelength domain, where vortex pairing occurs.Rahmani et al. (2016) show that while mixing and mixing efficiency are slightly increasedin low Reynolds number flows because high Prandtl number enhances three-dimensionalmotions, mixing and mixing efficiency are reduced due to smaller diffusivity in intermediateReynolds number flows. Salehiour & Peltier (2015b) studied the effect of Prandtl numberin high Reynolds number flows for Prandtl number up to 16. They show that mixingefficiency is reduced in higher Prandtl number flows and the dependence is weaker than inlow Reynolds number flows. In high Reynolds number flows, three-dimensional secondaryinstabilities prevent vortex pairing. However, in low to intermediate Reynolds numberflows, vortex pairing can increase the vertical scale (Corcos & Sherman, 1984; Smyth &Peltier, 1993; Moser & Rogers, 1993) of the motions and the amount of mixing (Rahmaniet al., 2014). In the present study, we investigate the dependence of vortex pairing, three-dimensional motions, and mixing on Prandtl number.The numerical methods and simulations are described in section 3.2. Section 3.3 dis-cusses the effect of Prandtl number on primary KH and vortex pairing in two-dimensionalflows to remove the effect of three-dimensional motions. Three-dimensional results aredescribed and compared with two-dimensional results in section 3.4.3.2 MethodsAs in Chapter 2, we use hyperbolic tangent functions as background density and velocityprofiles, i.e.,ρ¯ = −∆ρ2tanh(2zδ0)and U =∆U2tanh(2zh0). (3.1a, b).In all simulations, the initial density interface thickness, δ0, is set to be equal to the thicknessof the velocity profile, h0. We keep Reynolds number Re, bulk Richardson number Junchanged and vary Prandtl number Pr from 1 to 64 to examine the effect of Prandtlnumber. In this chapter, Re = 1200 and J = 0.07. For this combination of non-dimensionalnumbers, the background flow is susceptible to KH instabilities.Since density variation is small, the Boussinessq approximation is assumed. The govern-ing equations, the Navier-Stokes equations with the Boussinessq approximation, are solvedby the numerical code developed by Winters et al. (2004) and later improved by Smyth et al.(2005). The grid spacing of density field is half of that of velocity field in the code. The codeuses the third order Adams-Bashforth time stepping scheme. The boundary conditions arehorizontally periodic and vertically free slip and no flux.We run both two and three dimensional numerical simulations to investigate the effect ofPrandtl number on the flows. Inherently three-dimensional random perturbations are added323.2. MethodsRun Pr Nx Ny Nz tkh tsub t3d tfP1(2D) 1 640 - 640 92 126 – –P1(3D) 1 640 320 640 92 126 188 249P9(2D) 9 640 - 640 93 129 – –P9(3D) 9 640 320 640 92 130 180 238P16(2D) 16 640 - 640 92 130 – –P16(3D) 16 640 320 640 92 130 178 243P64(2D) 64 1024 - 1024 92 134 – –P64(3D) 64 1024 512 1024 92 134 182 248Table 3.1: Numerical parameters for all the simulations. The number of grid points is forthe density field and is double of the velocity field.to the background flow to excite primary and secondary instabilities for three-dimensionalsimulations. The perturbations are added to velocity fields and are given by equation (2.6)and (2.7). The maximum amplitude of perturbations is 0.08∆U , which is small enoughto ensure the flow evolves linearly initially. The same perturbations are used for differentPrandtl number simulations to remove the effect of initial conditions. The initial condi-tions in two-dimensional simulations are the spanwise averaged values in three-dimensionalsimulations.The horizontal domain length Lx is two wavelengths of the most unstable mode pre-dicted by Taylor-Goldstein equation to allow vortex pairing. The spanwise width for threedimensional simulations is one wavelength of the KH instability, which is at least six wave-lengths of the most unstable spanwise mode (for the most unstable spanwise wavenumbersee Klaassen & Peltier, 1985, 1989). The domain height is 15h0, which is sufficient to makethe flows unaffected by vertical boundaries. The smallest scale needed to be resolved inDNS for homogeneous turbulence is Kolmogrov scale, Lk = (ν/ε′)1/4, where ν is the molec-ular viscosity and ε′ is dissipation rate of turbulent kinetic energy. Moin & Mahesh (1998)suggest that grid spacing should be O(Lk) in DNS. The resolution for stratified flows withPr > 1 is however determined by the Bachelor scale, LB = Lk/√Pr. In this chapter,dissipation rate ε′ averaged within −h0/2 < z < h0/2 where turbulence is most energetic isused to calculate the Bachelor scale. For all simulations, ∆z/LB < 4.3. Grid numbers (fordensity field) and other numerical parameters are listed in table 3.1. All simulations arenamed Px(nD) where P indicates the simulations are investigating Prandtl number effects,x is the value of the Prandtl number, and n indicates whether the simulation is 2D or 3D.333.3. Two-dimensional results3.3 Two-dimensional resultsIn this section, we examine the effect of Prandtl number on KH instability and vortexpairing in two-dimensional simulations. The focus is the effect of Prandtl number on vortexpairing.3.3.1 KH instabilityWe first compare KH instability in simulations with different Prandtl number. Figure 3.1(a) shows the kinetic energy of the KH mode, Kkh. The first maximum of Kkh correspondsto the saturation of KH instabilities and this time is defined as the saturation time ofKH instabilities, tkh. The saturation time of KH instabilities is almost the same for allsimulations and shown as stars in the figure.0 25 50 75 100 125 150 175t00.511.522.53Kkh×10-3(a)P1(2D)P9(2D)P16(2D)P64(2D)0 25 50 75 100 125 150 175t00.0020.0040.0060.0080.01Ksub(b)Figure 3.1: (a) Kinetic energy of the KH mode, Kkh and (b) subharmonic mode, Ksub intwo-dimensional simulations. Saturation time of KH instability and subharmonic mode, tkhand tsub, are shown as stars and squares.Figure 3.1 shows that, before t = 50, the amplitude of Kkh is negligible for all thesimualtions. After t = 50, the difference between different Prandtl number simulationsappears. KH instabilities are always weaker for higher Prandtl number flows before thefirst maximum. To measure the thickness of the mean velocity profile and mean densityprofile, two integral scales are introduced by Smyth & Moum (2000), i.e. Iu and Iρ:Iu =∫ Lz/2−Lz/21−(2u¯∆U)2dz, and Iρ =∫ Lz/2−Lz/21−(2ρ¯∆ρ)2dz. (3.2)Bulk Richardson number based on these two scales, J(t) = ∆ρgI2u/ρ0Iρ, is used to measurethe strength of stratification. Figure 3.2 shows the evolution of Iu, Iρ, and J(t). Until343.3. Two-dimensional results0 25 50 75 100t11.251.51.7522.252.5I u/h0Pr = 1Pr = 9Pr = 16Pr = 640 25 50 75 100t11.251.51.7522.252.5I ρ/δ00 25 50 75 100t0.070.080.090.10.11J(t)Figure 3.2: The integral thickness of the mean velocity profile lu (a), thickness of the meandensity profile lρ (b), and bulk Richardson number J(t) based on lu and lρ.t = 50, the thickness of velocity profile is the same for all the simulations as KH mode isstill very weak and molecular diffusion is the main process expanding the mean velocityprofile. However, the thickness of density profile is higher in the Pr = 1 than the otherPr cases due to higher molecular diffusivity. For the other Pr cases, diffusion of the meandensity is negligible. Hence, before t = 50, bulk Richardson number J(t) is greater andthe stratification is stronger in higher Pr flows. Also, bulk Richardson number is muchlower in Pr = 1 than the other cases. Therefore, the reduction in Kkh is much larger fromPr = 1 to 9 than from Pr = 9 to 64. Figure 3.3 shows the density field at tkh for differentPrandtl number flows. It shows that from Pr = 1 to 9 the height of KH billows significantlydecreases, and from Pr = 9 to 64 the height decreases slightly.3.3.2 Vortex pairingFigure 3.1 (b) shows the kinetic energy of the subharmonic component Ksub. Before sat-uration of KH (shown as stars), Ksub is approximately the same in all the simulations.Therefore, the growth of the subharmonic mode is almost independent of Prandtl numberbefore the subharmonic mode becomes dominant.After saturation of KH instabilities, the difference between different Prandtl numbersimulations becomes noticeable. From Pr = 1 to 64, the growth of the subharmonic com-ponent decreases and the first maximum is reduced. The first (also global) maximum isdefined as the saturation of the subharmonic component and indicates the merging of twoKH billows. The reduction is more significant from Pr = 1 to 9 than from Pr = 9 to 64.Therefore, high Prandtl number has a suppressing effect on vortex pairing and the sup-pressing effect is more significant from Pr = 1 to 9 than from Pr = 9 to 64. As explainedin chapter 2, the mechanism of vortex pairing is that one KH vortex is displaced upwardand the other downward. Figure 3.2 shows that at tkh (around t = 92), bulk Richardsonnumber J increases as Pr increases because of slower diffusion of density. Hence, in high353.3. Two-dimensional resultsPrandtl number flows, KH vortices need to overcome stronger stratification to pair.Figure 3.3: Density field at four instances for two-dimensional simulations. From left toright, Pr = 1, 9, 16, and 64. Row (a) shows density field at t = 92, around saturationtime of KH instabilities for all cases. Row (b) shows the density field at t = 130, aroundsaturation of subharmonic mode. Row (c) shows the density field at t = 150 when themerged billow has not been destroyed by small scale motions. Row (d) shows the densityfield at t = 170 when the merged billow has been destroyed.The saturation time of the subharmonic mode is defined as tsub. Saturation time ofthe subharmonic mode slightly increases from Pr = 1 to 64 as shown in figure 3.1 (b) andtable 3.1. Hence, vortex pairing is slightly delayed in high Prandtl number flows. Figure3.3 (b) shows the density field at t = 130 (around tsub). In Pr = 1 simulation, the mergedbillow core is almost isopycnal due to diffusion. However, in higher Prandtl number flows,the structure of two KH billows are still identifiable due to lower diffusivity and the shape363.4. Three-dimensional resultsof the merged KH billows is less similar to an ellipse. The height of the merged billowdecreases as Prandtl number increases, indicating that stirring, i.e. the increase of totalpotential energy, decreases as Prandtl number increases.After two KH billows pair, Ksub oscillates with a high amplitude in Pr = 1. Thisoscillation is caused by the nutation of merged billow (Guha & Lawrence, 2014). However,in the other flows, the amplitude of oscillation is very small because small scale motionsdevelop and destroy the merged billow. Figure 3.3 (c) shows that at t = 150 there is still nosmall scale motions in Pr = 1 simulation and the merged billow is regular. As Pr increases,small scale motions become more dominant. After one oscillation, Ksub quickly drops andstops oscillating because the merged billow is destroyed in all simulations as shown in figure3.3 (d).0 25 50 75 100 125 150 17500.0020.0040.0060.0080.010 25 50 75 100 125 150 17500.0020.0040.0060.0080.010 25 50 75 100 125 150 17500.0020.0040.0060.0080.010 25 50 75 100 125 150 17500.0020.0040.0060.0080.01Figure 3.4: Kinetic energy of the subharmonic component, Ksub in two-dimensional sim-ulations and corresponding three-dimensional simulations: (a) Pr = 1; (b) Pr = 9; (c)Pr = 16; (d) Pr = 64.3.4 Three-dimensional resultsIn this section, we examine effects of Prandtl number on vortex pairing and mixing whenthree-dimensional motions exist. As Caulfield & Peltier (2000) and Rahmani et al. (2014)373.4. Three-dimensional resultshave shown, before saturation of KH instabilities, three-dimensional motions can be ne-glected if initial perturbations are small enough. We first compare vortex pairing in thetwo-dimensional results with that in the three-dimensional results to illustrate the effectof three-dimensional motions. Then we investigate the effect of Prandtl number on three-dimensional motions. We finally examine how mixing properties are influenced by themolecular property, i.e., Prandtl number.3.4.1 Vortex pairingWe compare the kinetic energy of the subharmonic component, Ksub, in two-dimensionalsimulations with that in three-dimensional simulations in figure 3.4. For all Prandtl numbersimulations, the maximum of Ksub is reduced by three-dimensional motions. Nevertheless,the saturation time of the subharmonic component is not affected by three-dimensionalmotions in all simulations. Hence, three-dimensional motions can suppress pairing, butdo not delay pairing. Also, the reduction of maximum Ksub caused by three-dimensionalmotions does not show a clear trend with Prandtl number.3.4.2 Kinetic energyWe examine the three kinetic energy components, i.e. K¯, K2d, and K3d, in all three-dimensional simulations with different Prandtl number in figure 3.5. For all simulations, be-fore saturation of KH instabilities, three-dimensional motions are negligible. Two-dimensionalmotions are dominated by the KH mode. As shown in the two-dimensional results, KH in-stabilities are stronger in lower Prandtl number flows. Moreover, more energy is transferredto two-dimensional kinetic energy K2d from the mean kinetic energy K in lower Prandtlnumber flows. Also as explained in the two-dimensional results, the effects of Prandtl num-ber on the flows during this stage are mainly caused by the differing diffusion rate of densityfield before KH instabilities reach finite amplitude. The three kinetic energy componentsare similar for flows of Pr = 9 to Pr = 64 as diffusion is negligible before KH instabilitiesreach finite amplitude.After the saturation of KH instabilities and before vortex pairing, three-dimensionalmotions gradually become important, though two-dimensional motions are still strongerthan three-dimensional motions. At the time of pairing, K3d is slightly higher for higherPr number flows. This is consistent with Klaassen & Peltier (1985), who find that thegrowth rate of the dominant three-dimensional secondary instability is higher for higherPr. However, two-dimensional kinetic energy K2d is higher in lower Pr flows. Hence, highPr enhances three-dimensional motions, however KH instabilities and vortex pairing aresuppressed by high Pr. Also, it should be noted that the mean kinetic energy is mainly383.4. Three-dimensional results0 50 100 150 200 250 300 350 4000.60.70.80.91K/K0(a)Pr = 1Pr = 9Pr = 16Pr = 640 50 100 150 200 250 300 350 40000.020.040.060.080.1K2d/K0(b)0 50 100 150 200 250 300 350 400t00.020.040.060.080.1K3d/K0(c)Figure 3.5: The three kinetic energy components, (a) K, (b) K2d, and (c) K3d, in three-dimensional simulations. Three key times, tkh, tsub, t3d, are shown as stars, squares, andcircles.393.4. Three-dimensional resultslost due to vortex pairing rather than KH instabilities.After the time of pairing, tsub, for flows where Pr is higher than 1, three-dimensionalkinetic energy K3d continues to increase and reaches its maximum value. As Pr increases,the maximum of K3d also increases. Note that maximum K3d is less than maximum of K2dfor Pr = 1 flow and the other way around for higher Pr flows. Therefore, three-dimensionalmotions are more dominant in higher Pr flows. At t3d, more mean kinetic energy is lost forhigher Pr flows because of stronger three-dimensional motions, though less mean kineticenergy is lost at the time of pairing for higher Pr flows.After three-dimensional kinetic energy reaches its maximum, turbulence decays and theflows re-larminarize. In the end, more kinetic energy is lost and velocity field is diffusedmore in higher Pr flows. The kinetic energy is either lost to heat or transferred to potentialenergy as discussed in the next section.3.4.3 MixingFigure 3.6 (a) shows the increase of potential energy normalized by initial kinetic energy and(b) shows the normalized background potential energy. Before saturation of KH instabilities,the amount of stirring, ∆P/K0, decreases as Pr increases, as the KH wave kinetic energyKkh does. The decreasing trend of stirring with Pr can be explained by the lower height ofKH billows in higher Pr simulations as shown in figure 3.3 (a). During this stage, mixing ismainly caused by molecular diffusion as the flow is still well organized. In addition to loweramount of stirring in higher Pr flows, lower molecular diffusivity decreases the amount ofmixing further. The amount of mixing during this stage is negligible for all simulations.After saturation of KH instabilities, potential energy continues to increase due to vortexpairing. The first maximum of potential energy occurs after time of pairing because therotation of merged billow brings more heavy fluid particles upward as shown in figure 3.3 (b)and (c) and increases potential energy. The first maximum of potential energy is lower forhigher Pr flows since the size of the merged KH billow is smaller for higher Pr flows. Theamount of mixing occurred during this stage is dominant by molecular diffusion as three-dimensional motions are still small. Hence, less stirring and lower molecular diffusivity leadto less mixing in higher Pr flows. Also, the amount of mixing is much lower for Pr > 1than Pr = 1.After the first maximum of potential energy, potential energy decreases and reaches alocal minimum at t3d. During this stage, two-dimensional kinetic energy K2d also decreases,but three-dimensional kinetic energy increases. Hence, the energy stored in available poten-tial energy is transferred to three-dimensional kinetic energy K3d through buoyancy flux.During this stage, the merged KH billow is destroyed and small scale motions become preva-lent. Also, mixing occurs at a higher rate due to turbulent mixing than the previous stage.403.4. Three-dimensional resultsThe amount of mixing still decreases as Pr increases.After three-dimensional kinetic energy reach its maximum, part of kinetic energy istransferred to potential energy and potential energy increases again. After some oscillation,potential energy gradually reaches a constant as the flow returns to an approximate laminarstate. The energy stored in available potential energy is transferred to background potentialenergy by mixing and background potential energy increases to total amount of potentialenergy. From Pr = 1 to 9, the amount of mixing significantly decreases because vortexpairing is suppressed by higher Pr and less stirring occurs in Pr = 9. Although three-dimensional motions are stronger in Pr = 9, the effect of weaker pairing on mixing is moresignificant. From Pr = 9 to Pr = 16, mixing slightly increases due to stronger three-dimensional motions, although pairing is weaker in the latter case. From Pr = 16 to 64,mixing is almost the same.0 50 100 150 200 250 300 350 400 450t00.020.040.060.080.1∆P/K0(a)Pr = 1Pr = 9Pr = 16Pr = 640 50 100 150 200 250 300 350 400 450t00.020.040.060.080.1∆Pb/K0(b)Figure 3.6: (a) The increase of potential energy normalized by initial kinetic energy,∆P/K0. (b) The increase of background potential energy normalized by initial kineticenergy, ∆Pb/K0. Three key times, tkh, tsub, t3d are shown as stars, squares, and circles.413.4. Three-dimensional results3.4.4 Mixing efficiencyIn this section, we first compare the evolution of mixing efficiency in different Pr flows infigure 3.7 (a). Following Caulfield & Peltier (2000), the instantaneous mixing Ei is definedasEi =φMφM + ε. (3.3)Mixing efficiency Ei is negligible before t = 50 as KH mode has not reached finite amplitude.After that, mixing efficiency begins to rise for all Pr flows due to the rolling up of thedensity interface. During this stage, mixing efficiency Ei decreases as Pr increases, becausemolecular diffusion is the dominant mechanism of mixing. Even until around t = 180, thisrelation between Ei and Pr is valid. After that, the flow becomes fully turbulent and mixingefficiency is around 0.2 as observations and numerical simulations have shown. Mixingefficiency for Pr = 1 is larger than that for the other Pr. However, there is not a clear trendfrom Pr = 9 to 64. Therefore, we compare cumulative mixing efficiency Ec during the activeturbulent stage t ∈ [t3d, tf ] in figure 3.7 (b). From Pr = 1 to 9, cumulative mixing efficiencysignificantly decreases. Nevertheless, from Pr = 9 to Pr = 16, Ec slightly increases.Although this non-monotonic trend is inconsistent with existed studies, for example Smyth& Moum (2000), Salehiour & Peltier (2015b), and Rahmani et al. (2016) , it is still possiblebecause mixing efficiency also depends on Reb and Ri. From Pr = 9 to 64, the change ofmixing efficiency is much less than that from Pr = 1 to 9.0 50 100 150 200 250 300t00.10.20.30.40.50.60.7EiPr = 1Pr = 9Pr = 16Pr = 640 10 20 30 40 50 60 70Pr0.20.2050.210.2150.220.2250.230.2350.24EcFigure 3.7: (a) Instantaneous mixing efficiency Ei and (b) cumulative mixing efficiency Ecduring the period t ∈ [t3d, tf ] .423.5. Conclusion and discussion3.5 Conclusion and discussionWe ran two-dimensional and three-dimensional simulations to examine the effect of Prandtlnumber on KH instabilities, vortex pairing, and mixing. Two-dimensional results show thatKH instabilities and vortex pairing are suppressed in high Prandtl number flows and thesuppression is more significant from Pr = 1 to 9 than from 9 to 64. The suppressing effectof higher Prandtl number on KH instabilities is mainly due to diffusion occurring beforeKH instabilities reach finite amplitude as Klaassen & Peltier (1985) propose. This effectdoes not appear in Salehiour & Peltier (2015b) and Rahmani et al. (2016) because theeigenfunction of KH mode is added to background flow in these two studies rather than apure random perturbation.As shown in the previous chapter, three-dimensional motions can suppress vortex pair-ing. Nevertheless, the suppressing effect of three-dimensional motions on vortex pairing doesnot show a clear trend with Pr. Since KH instabilities and vortex pairing are suppressed inhigher Prandtl number flows, the amount of stirring and mixing caused by these two eventsare also lower in higher Prandtl number flows. However, three-dimensional motions increasemonotonically and tend to increase mixing as Pr increases. Because these two contradictingeffects on mixing, mixing for Pr = 1 is much higher than the other Pr flows where mixing isalmost the same. It should be noted that three-dimensional motions are enhanced by higherPr is only valid in certain non-dimensional parameter space. Rahmani et al. (2016) findthat three-dimensional motions are enhanced in low Reynolds number flows and Salehiour& Peltier (2015b) shows that in high Reynolds number flows three-dimensional motions donot vary monotonically as Pr increases. In the present study, the growth rate of the convec-tive core instability (also the dominant three-dimensional secondary instabilities), increasesas Pr increases (Klaassen & Peltier, 1985). For other initial non-dimensional parameters,the growth rate of other secondary instabilities might be lower for higher Pr (see Mashayek& Peltier, 2012). Similar to mixing, cumulative mixing efficiency in Pr = 1 during theactive turbulent stage is much higher than the other Pr flows. Overall, cumulative mixingefficiency shows a decreasing trend with Pr.43Chapter 4ConclusionsKelvin-Helmholtz instabilities in stratified flows have been studied by running Direct Nu-merical Simulations. The effects of initial conditions and Prandtl number on vortex pairingand mixing have been investigated. In chapter 2, we compare the numerical results witha variety of initial perturbations. Klaassen & Peltier (1985) and Smyth & Peltier (1993)show that vortex pairing is delayed if the subharmonic mode is out of phase relative to theKH mode. We find that the delay increases as the phase difference increases by runningsimulations with eigenfunctions of the KH and subharmonic modes. Moreover, the delay issensitive to the phase and other details of the initial perturbations as the phase differencebecomes close to ±pi2 . In three-dimensional simulations, vortex pairing is always suppressedby three-dimensional motions and the suppression is greater if the phase difference is larger.Vortex pairing may even be eliminated by three-dimensional motions if the subharmonicmode is out of phase. In the case with the subharmonic mode being out of phase, mixing issignificantly reduced because of no pairing. The reduction of mixing also increases as thephase difference increases.In chapter 3, we examine the effect of Prandtl number on vortex pairing and mixing.Two-dimensional simulation results show that KH instabilities are suppressed in higherPrandtl number flows due to less diffusion of the density field before KH mode reachesfinite amplitude as proposed by Klaassen & Peltier (1985). The suppressing effects be-come weaker as Prandtl number increase from 9 to 64 than from 1 to 9. Comparison ofthree-dimensional and two-dimensional results shows that three-dimensional motions cansuppress vortex pairing,but not delay pairing. The suppression of three-dimensional mo-tions on vortex pairing does not show a clear trend with Prandtl number. However, three-dimensional motions are stronger in higher Prandtl number flows. This may be only validfor certain non-dimensional numbers as different secondary three-dimensional instabilitiesmay vary differently with Prandtl number. Whether mixing is decreased or increased asPrandtl number increases depend on the two contradicting effects: weaken pairing; enhancethree-dimensional motions. As Prandtl number increases from 1 to 9, vortex pairing issignificantly weakened, which results in less stirring. Hence, mixing decreases as Prandtlnumber increases from 1 to 9. From Pr = 9 to 64, mixing is almost the same. We finallyexamine the effect of Prandtl number on mixing efficiency. Mixing efficiency in Pr = 1 is44Chapter 4. Conclusionsmuch higher than that in the other Pr flows. Overall, mixing efficiency shows a decreasingtrend with Pr.45BibliographyAlexakis, A. (2009). Stratified shear flow instabilities at large Richardson numbers.Phys. Fluids., 21:054108.Arratia, C., Caulfield, C. P., & Chomaz, J. M. (2013). Transient perturbation growth intime-dependent mixing layers. J. Fluid Mech, 717(90-133).Baines, P. G. (1994). On the mechanism of shear flow instabilities. J. Fluid Mech., 276:327–342.Barry, M., Ivey, G. N., Winters, K. B., & Imberger, J. (2001). Measurements of diapycnaldiffusivities in stratified fluids. J. Fluid Mech, 442:267–291.Batchelor, G. K. (1959). Small-scale variation of convected quantities like temperature inturbulent fluid Part 1. General discussion and the case of small conductivity. J. FluidMech., 5:113–133.Browand, F. K. & Winant, C. D. (1973). Laboratory observations of shear-layer instabilityin a stratified fluid. Boundary-Layer Meteorology., 5:67–77.Carpenter, J. R., Tedford, E. W., Heifetz, E., & Lawrence, G. A. (2011). Instability instratified shear flow: Review of a physical interpretation based on interacting waves.Applied Mechanics Reviews, 64(6):060801.Carpenter, J. R., Tedford, E. W., Rahmani, M., & Lawrence, G. A. (2010). Holmboe wavefields in simulation and experiment. J. Fluid Mech., 648:205–223.Caulfield, C. P. (1994). Multiple linear instability of layered stratified shear flow. J. FluidMech., 258:255–285.Caulfield, C. P. & Peltier, W. R. (2000). The anatomy of the mixing transition in homoge-neous and stratified free shear layers. J. Fluid Mech., 413:1–47.Constantinoul, N. C. & Ioannou, P. J. (2011). Optimal excitation of two dimensionalholmboe instabilities. Physics of Fluids, 23(7):074102.46BibliographyCorcos, G. M. & Sherman, F. S. (1976). Vorticity concentration and the dynamics ofunstable free shear layers. J. Fluid Mech, 73:241–264.Corcos, G. M. & Sherman, F. S. (1984). The mixing layer: deterministic models of aturbulent flow. Part 1. introduction and the two-dimensional flow. J. Fluid Mech., 139:29–65.Drazin, P. G. & Reid, W. H. (2004). Hydrodynamic stability. Cambridge university Press.Farrell, B. F. & Ioannou, P. J. (1996). Generalized stability theory. Part I: Autonomousoperators. J. Atmos Sci., 53(14):2025–2040.Gargett, A. E. (1988). The scalling of turbulence in the presence of stable stratification. J.Phys. Oceanogr., 93:5021–5036.Geyer, W. R., Lavery, A. C., Scully, M. E., & Trowbridge, J. H. (2010). Mixing by shearinstability at high Reynolds number. J. Geophys. Res., 37:L22–607.Goldstein, S. (1931). On the stability of superposed streams of fluids of different densities.Proc. R. Soc. Lond., 132:524–548.Guha, A. & Lawrence, G. A. (2014). A wave interaction approach to studying non-modalhomogeneous and stratified shear instabilities. J. Fluid Mech., 755:336–364.Hazel, P. (1972). Numerical studies of the stability of inviscid stratified shear flows. J. FluidMech, 51:39–61.Ho, C.-M. & Huang, L.-S. (1982). Subharmonics and vortex merging in mixing layers.J. Fluid Mech, 119:443–473.Ho, C.-M. & Huerre, P. (1984). Perturbed free shear layers. Annu. Rev. Fluid Mech.,16:365–424.Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geophys.Publ, 24:67–113.Howard, L. (1961). Note on a paper of John W. Miles. J. Fluid Mech., 10:509–512.Ivey, G. N., Winters, K. B., & Koseff, J. R. (2008). Densifty stratification, turbulence, buthow much mixing? Annu. Rev. Fluid Mech, 40:169–184.Kaminski, A. K., Caulfield, C. P., & Taylor, J. R. (2014). Transient growth in stronglystratified shear layers. J. Fluid Mech, 758(R4).47BibliographyKlaassen, G. P. & Peltier, W. R. (1985). The effect of Prandtl number on the evolution andstability of Kelvin-Helmholtz billows. Geophys. Astrophys. Fluid Dynamics, 32:23–60.Klaassen, G. P. & Peltier, W. R. (1989). The role of transverse secondary instabilities inthe evolution of free shear layers. J. Fluid Mech., 202:367–402.Koop, C. G. & Brownad, F. K. (1979). Instability and turbulence in a stratified fluid withshear. J. Fluid Mech, 93:135–159.Koppel, D. (1964). On the stability of flow of a thermally stratified fluid under the actionof gravity. J. Math. Phys., 5:963–982.Kundu, P. K. & Cohen, I. M. (2008). Fluid Mechanics. Academic Press.Mashayek, A. & Peltier, W. R. (2012). The ‘zoo’ of secondary instabilities percursoryto stratified shear flow transition. Part 1 shear aligned convection, pairing, and braidinstabilities. J. Fluid Mech, 708:5–44.Mashayek, A. & Peltier, W. R. (2013). Shear-induced mixing in geophysical flows: does theroute to turbulence matter to its efficiency? J. Fluid Mech., 725:216–261.Metcalfe, R. W., Orszag, S. A., Brachet, M. E., & Riley, J. J. (1987). Secondary instabilityof a temporally growing mixing layer. J. Fluid Mech., 184:207–243.Miles, J. (1961). On the stability of heterogeneous shear flows. J. Fluid Mech., 10:496–508.Moin, P. & Mahesh, K. (1998). Direct numerication simulation: a tool in turbulenceresearch. Annu. Rev. Fluid Mech., 30:539–578.Moser, R. D. & Rogers, M. M. (1993). The three-dimensional evolution of a plane mixinglayer: pairing and transition to turbulence. J. Fluid Mech., 247:275–320.Moum, J. N. & Farmer, D. M. (2003). Structure and generation of turbulence at interfacesstrained by internal solitary waves propogating shoreward over the continental shelf.J. Phys. Oceanogr., 33:2093–2112.Osborn, T. R. (1980). Estimates of the local rate of vertical diffusion from dissipationmeasurements. J. Phys. Oceanogr., 10(83–89).Patnaik, P. C., Sherman, F. S., & Corcos, G. M. (1976). A numerical simulation of Kelvin-Helmholtz waves of finite amplitude. J. Fluid Mech., 73:215–240.Peltier, W. R. & Caulfield, C. P. (2003). Mixing efficiency in stratified shear flows.Annu. Rev. Fluid Mech., 35:135–167.48BibliographyRahmani, M., Lawrence, G. A., & Seymour, B. R. (2014). The effect of Reynolds numberon mixing in Kelvin-Helmholtz billows. J. Fluid Mech., 759:612–641.Rahmani, M., Seymour, B. R., & Lawrence, G. A. (2016). The effect of Prandtl number onmixing in low Reynolds number Kelvin-Helmholtz billows. Phys. Fluids., 28:054107.Rehmann, C. R. & Koseff, J. R. (2004). Mean potential energy change in stratified gridturbulence. Dyn. Atmos. Oceans., 37(271–294).Salehiour, H. & Peltier, W. R. (2015a). Diapycnal diffusivity, turbulent Prandtl numberand mixing efficiency in Boussinessq stratified turbulence. J. Fluid Mech, 775(464–500.).Salehiour, H. & Peltier, W. R. (2015b). Turbulent dispycnal mixing in stratified shear flows:the influence of Prandtl number on mixing efficiency and transition at high Reynoldsnumber. J. Fluid Mech, 773:187–223.Salehiour, H., Peltier, W. R., Whalen, C. B., & MacKinnon, J. A. (2016). A new charac-terization of the turbulent diapycnal diffusivities of mass and momentum in the ocean.Geophys. Res. Lett, 43(3370–3379).Scharnhorst, K. (2001). Angles in complex vector spaces. Acta Applicandae Mathematica.,69(1):95–103.Schmid, P. J. & Henningson, D. S. (2012). Stability and transition in shear flows. SpringerScience & Business Media.Seim, H. E. & Gregg, M. C. (1994). Detailed observations of a naturally occuring shearinstability. J. Geophys. Res., 99(C5):10 049–10 073.Shih, L. H., Koseff, J. R., Ivey, G. N., & Ferziger, J. H. (2005). Parameterization ofturbulence fluxes and scales using homogeneous sheared stably stratified turbulence sim-ulations. J. Fluid Mech., 525:193–214.Smyth, W. D. & Moum, J. N. (2000). Length scales of turbulence in stably stratified mixinglayers. Physics of Fluids, 12:1327.Smyth, W. D., Moum, J. N., & Caldwell, D. R. (2001). The efficiency of mixing in turbulentpatches: Inferences from direct simulations and microstructure observations. Journal ofPhysical Oceanography, 31(8):1969–1992.Smyth, W. D., Nash, J. D., & Moum, J. N. (2005). Differential diffusion in breakingKelvin-Helmholtz billows. J. Phy. Oceanogr., 35:1004–1022.49BibliographySmyth, W. D. & Peltier, W. R. (1993). Two-dimensional turbulence in homogenous andstratified shear layers. Geophs. Astrophys. Fluid Dynamics., 69:1–32.Smyth, W. D. & Winters, K. B. (2003). Turbulence and mixing in Holmboe waves.J. Phy. Oceanogr., 33:694–711.Taylor, G. I. (1931). Effect of variation in density on the stability of superposed streams offluid. Proc. R. Soc. Lond., 132:499–523.Thorpe, S. A. (1971). Experiments on the instability of stratified shear flows: misciblefluids. J. Fluid Mech., 46:299–319.Thorpe, S. A. (1973). Experiments on instability and turbulence in stratified shear flow.J. Fluid Mech., 61:731–751.Venaille, A., Gostiaux, L., & Sommeria, J. (2017). A statistical mechanics approach tomixing in stratified fluids. J. Fluid Mech., 810:554–583.Winant, C. D. & Browand, F. K. (1974). Vortex pairing: the mechanism of turbulentmixing-layer growth at moderate Reynolds number. J. Fluid Mech, 63:237–255.Winters, K. B., , & D’asaro, E. A. (1996). Diascalar flux and the rate of fluid mixing.J. Fluid Mech., 317:179–193.Winters, K. B., Lombard, P. N., Riley, J. J., & D’asaro, E. A. (1995). Availabel potentialenergy and mixing in density-stratified fluids. J. Fluid Mech., 289:115–128.Winters, K. B., MacKinnon, J. A., & Mills, B. (2004). A spectral model for process studiesof rotating density-stratified flows. J. Atmos. Ocean. Technol., 21:69–94.50