UNSTABLE WAVES ON A SHEARED DENSITY INTERFACE by Jeffrey Richard Carpenter M.A.Sc., University of British Columbia, 2005 B.Sc.Eng., University of Guelph, 2003 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2009 c© Jeffrey Richard Carpenter 2009 ii Abstract The Holmboe instability is known to occur in stratified shear layers that exhibit a relatively thin density interface compared to the shear layer thickness. At finite amplitude the instability appears as cusp-like propagating internal waves. The evolution and identification of these unstable waves is the subject of this thesis. The results are presented in four parts. First, the basic wave field resulting from Holmboe’s instability is studied both numerically and experimentally. In comparing basic descriptors of the wave fields, a number of processes are identified that are responsible for differences between the simulations and experiments. These are related to variations in the mean flow that arise due to the different boundary conditions. Holmboe waves are known to produce vertical ejections of interfacial fluid from the wave crests. This ‘ejection process,’ in which stratified fluid is transported against buoyancy forces, is caused by the formation of a vortex couple (i.e. two vorticies of opposite sign that travel as a pair). Results obtained by means of direct numerical simulations also show that the process is primarily two-dimensional and does not require the presence of both upper and lower Holmboe modes. Shear instability is then studied in the highly stratified Fraser River estuary. The observations are found to be in good agreement with the predictions of linear theory. When instability occurs, it is largely as a result of asymmetry between regions of strong shear and density stratification. The structure of the salinity intrusion is found to depend on the strength of the freshwater discharge, in addition to the phase of the tidal cycle. This has implications for whether estuarine mixing takes place through shear instability or boundary layer turbulence. Finally, the asymmetric stratified shear layer, which exhibits a vertical shift between the den- sity interface and the shear layer centre, is examined by the formulation of a diagnostic that is based on the ‘wave interaction’ mechanism of instability growth. This allows for a quantitative assessment of Kelvin-Helmhotz and Holmboe-type growth mechanisms in stratified shear layers. The predictions of the diagnostic are compared to results of nonlinear simulations and observa- tions in the Fraser River estuary. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Co-authorship statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction and overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Holmboe wave fields in simulation and experiment . . . . . . . . . . . . . . . . . . 13 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Linear stability of stratified shear layers . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Description of the numerical simulations . . . . . . . . . . . . . . . . . 15 2.3.2 Description of the laboratory experiment . . . . . . . . . . . . . . . . . 17 2.4 Wave structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Phase speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Spectral evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6.1 Frequency shifting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6.2 Wave energy spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7 Wave growth and amplitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.7.1 Wave growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.7.2 Comparison of saturated amplitudes . . . . . . . . . . . . . . . . . . . . 28 2.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3 Ejection process in Holmboe waves . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Ejection mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 iv Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Shear instability in the Fraser River estuary . . . . . . . . . . . . . . . . . . . . . 43 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Site Description and Data Collection . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3 General Description of the Salinity Intrusion . . . . . . . . . . . . . . . . . . . . 46 4.3.1 High Discharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3.2 Low Discharge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.4 Stability of Stratified Shear Flows . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4.1 Taylor-Goldstein Equation . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4.2 Miles-Howard criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4.3 Mixing Layer Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4.4 Solution of the TG equation for observed profiles . . . . . . . . . . . . . 52 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.6 Small scale overturns and bottom stress . . . . . . . . . . . . . . . . . . . . . . 59 4.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5 Unstable modes in asymmetric stratified shear layers . . . . . . . . . . . . . . . . 69 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5.2 Wave interaction interpretation of instability . . . . . . . . . . . . . . . . . . . . 73 5.3 Formulation of a diagnostic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3.1 Taylor–Goldstein equation . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3.2 Partitioning into kinematic and baroclinic fields . . . . . . . . . . . . . . 77 5.3.3 Piecewise profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3.4 Smooth profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3.5 Classification of modes . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4.1 Symmetric profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.4.2 Asymmetric profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.5 Application to field profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.6 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 vList of Tables 2.1 Important parameters in simulations and experiment . . . . . . . . . . . . . . . . 17 4.1 Transect details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 vi List of Figures 1.1 Profiles of a Kelvin-Helmholtz unstable stratified shear layer . . . . . . . . . . . 4 1.2 Profiles of a stratified shear layer with R = 3 . . . . . . . . . . . . . . . . . . . 5 1.3 Profiles of an asymmetric stratified shear layer (a = 0.5, R = 3) . . . . . . . . . 6 2.1 Growth rate and dispersion relations for the DNS . . . . . . . . . . . . . . . . . 15 2.2 Evolution of the background profiles . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Spatial changes in U(z) and layer depths in the laboratory . . . . . . . . . . . . 18 2.4 Sample density fields from the laboratory experiments and the simulations . . . . 19 2.5 Wave characteristics in the laboratory experiments and the simulations . . . . . . 21 2.6 Traced wave characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Holmboe wave pairing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.8 Spectral evolution in the simulation and experiment . . . . . . . . . . . . . . . . 25 2.9 Growth of kinetic energy in the simulations . . . . . . . . . . . . . . . . . . . . 28 3.1 Laboratory photograph of a Holmboe wave ejection . . . . . . . . . . . . . . . . 36 3.2 Illustration of Holmboe wave formation . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Time sequence of ejection process . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Wave interactions in the ejection process . . . . . . . . . . . . . . . . . . . . . . 39 4.1 Map of the lower 27 km of the Fraser River . . . . . . . . . . . . . . . . . . . . 44 4.2 Observed tides at Point Atkinson and New Westminster . . . . . . . . . . . . . . 45 4.3 Echo soundings observed during high discharge . . . . . . . . . . . . . . . . . . 47 4.4 Echo soundings observed during low discharge . . . . . . . . . . . . . . . . . . 49 4.5 Transect 1 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.6 Transect 2 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.7 Transect 3 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.8 Transect 4 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.9 Transect 5 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.10 Transect 6 stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.11 Selected density profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.12 KH and one-sided instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1 Stratified shear layer profiles considered . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Stability of the symmetric case . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.3 Stability of the asymmetric case . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4 KH and H stability diagrams . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.5 Partial growth rates of piecewise symmetric profiles . . . . . . . . . . . . . . . . 81 5.6 Partial growth rates of smooth symmetric profiles . . . . . . . . . . . . . . . . . 82 5.7 Partial growths along maximum growth curve . . . . . . . . . . . . . . . . . . . 83 5.8 Partial growth rates of piecewise asymmetric profiles . . . . . . . . . . . . . . . 84 5.9 Partial growth rates of smooth asymmetric profiles . . . . . . . . . . . . . . . . 85 5.10 Partial growths along maximum growth curve . . . . . . . . . . . . . . . . . . . 85 5.11 Profiles from field measurements . . . . . . . . . . . . . . . . . . . . . . . . . . 87 vii List of Symbols Roman Characters: a . . . . . . . . . . . Dimensionless asymmetry, wave amplitude arms . . . . . . . . Root-mean-square wave amplitude c . . . . . . . . . . . Complex phase speed cd . . . . . . . . . . Speed of internal wave on density interface cg . . . . . . . . . . Group speed ci . . . . . . . . . . . Imaginary component of the complex phase speed cr . . . . . . . . . . Real component of the complex phase speed cv . . . . . . . . . . Vorticity wave speed D . . . . . . . . . . Interfacial region domain d . . . . . . . . . . . Dimensional asymmetry offset E . . . . . . . . . . .Wave energy density F . . . . . . . . . . .Weight function G . . . . . . . . . . .Green’s function g . . . . . . . . . . . Gravitational acceleration g′ . . . . . . . . . . .Reduced gravitational acceleration h . . . . . . . . . . . Shear layer thickness i . . . . . . . . . . . . Imaginary unit, √−1 J . . . . . . . . . . . Bulk Richardson number K . . . . . . . . . . .Wave kinetic energy K0 . . . . . . . . . . Initial kinetic energy K2d . . . . . . . . . Two-dimensional kinetic energy K3d . . . . . . . . . Three-dimensional kinetic energy k . . . . . . . . . . . Dimensional wavenumber kmax . . . . . . . .Dimensional wavenumber of maximum growth Lx . . . . . . . . . . Domain length in the streamwise direction Ly . . . . . . . . . . Domain length in the spanwise direction Lz . . . . . . . . . . Domain length in the vertical direction m . . . . . . . . . . Density interface index N2 . . . . . . . . . Squared buoyancy frequency N2∗ . . . . . . . . . Modified square buoyancy frequency n . . . . . . . . . . . Vorticity interface index Pr . . . . . . . . . Prandtl number p . . . . . . . . . . . Pressure, interface index p̃ . . . . . . . . . . . Perturbation pressure p̄ . . . . . . . . . . . Hydrostatic background pressure q . . . . . . . . . . . Perturbation vorticity qK . . . . . . . . . .Kinematic portion of perturbation vorticity qB . . . . . . . . . . Baroclinic portion of perturbation vorticity R . . . . . . . . . . .Scale ratio viii Re . . . . . . . . . .Reynolds number Ri . . . . . . . . . . Gradient Richardson number r . . . . . . . . . . . Shear layer spreading rate S . . . . . . . . . . . Wave stretching t . . . . . . . . . . . .Time U(z) . . . . . . . .Background velocity profile U∗(z) . . . . . . . Modified background velocity profile U1, U2 . . . . . . Layer velocities ∆U . . . . . . . . . Velocity difference between layers Ū . . . . . . . . . . .Velocity of the mean flow u . . . . . . . . . . . Velocity component in the streamwise direction u . . . . . . . . . . . Velocity field vector ũ . . . . . . . . . . . Perturbation horizontal velocity u1d . . . . . . . . . One-dimensional velocity u2d . . . . . . . . . Two-dimensional velocity u3d . . . . . . . . . Three-dimensional velocity vb . . . . . . . . . . Velocity of the boat w . . . . . . . . . . . Velocity component in the vertical direction w̃ . . . . . . . . . . . Vertical perturbation velocity x . . . . . . . . . . . Streamwise component of Cartesian coordinate system y . . . . . . . . . . . Spanwise component of Cartesian coordinate system z . . . . . . . . . . . Vertical component of Cartesian coordinate system Greek Characters: α . . . . . . . . . . . Dimensionless wavenumber αmax . . . . . . . Dimensionless wavenumber of maximum growth δ . . . . . . . . . . . Density interface thickness η . . . . . . . . . . . Interface elevation, displacement field η̂ . . . . . . . . . . . Displacement eigenfunction η̂∗ . . . . . . . . . . Modified displacement eigenfunction θ . . . . . . . . . . . Phase function κ . . . . . . . . . . . Diffusion coefficient of density λ . . . . . . . . . . . Wavelength λ∗ . . . . . . . . . . Apparent wavelength ν . . . . . . . . . . . Coefficient of kinematic viscosity ρ . . . . . . . . . . . Density ρ0 . . . . . . . . . . Reference density ∆ρ . . . . . . . . . Density difference between layers ρ1, ρ2 . . . . . . . Layer densities ρ̄(z) . . . . . . . . Density profile σ . . . . . . . . . . . Intrinsic frequency, dimensionless growth rate σT . . . . . . . . . . Density less 1000 kg m−3 σKH . . . . . . . . KH partial growth rate σH . . . . . . . . . H partial growth rate ix σself . . . . . . . . Self interaction partial growth rate ψ . . . . . . . . . . . Perturbation stream function ψ̂ . . . . . . . . . . . Perturbation stream function eigenfunction ψK . . . . . . . . . Kinematic portion of perturbation stream function ψ̂K . . . . . . . . . Kinematic portion of perturbation stream function eigenfunction ψB . . . . . . . . . Baroclinic portion of perturbation stream function ψ̂B . . . . . . . . . Baroclinic portion of perturbation stream function eigenfunction ω . . . . . . . . . . . Wave frequency List of Abbreviations: DNS . . . . . . . . Direct numerical simulation KH . . . . . . . . . Kelvin-Helmholtz H . . . . . . . . . . . Holmboe T . . . . . . . . . . . Taylor TG . . . . . . . . . Taylor-Goldstein xAcknowledgements I would like to thank both of my supervisors Greg Lawrence and Neil Balmforth for giving me a great deal of support and advice. The interactions that I’ve had with them has changed my attitude and approach to research and teaching for the better. I also want to thank Ted Tedford for suggesting that we work together on a couple of confer- ence papers, which turned out to be a sizable portion of this thesis, and an enjoyable experience. Half of that collaboration would not have been possible if Rich Pawlowicz had not invited us to take part in the Fraser River project. Roger Pieters has also been a great source of help in all sorts of areas during my stay at UBC. The numerical model that has been used to generate a signif- icant portion of the results in this thesis was kindly provided by Kraig Winters and Bill Smyth. In addition, Bill Smyth has provided the Taylor-Goldstein code, as well as contributed valuable discussions in the early stages of this thesis. I’ve had the opportunity to take part in a number of informal seminar series which have given me excellent feedback on my work. I’d therefore like to thank the Complex Fluids group, EFM group, and the Physical Oceanography group for their feedback. And of course thank you to my family and friends. xi Co-authorship statement The four chapters comprising the main body of this thesis have each been part of a collaborative effort. Each is intended, or has already been submitted, as a journal publication. The contributions of each author involved are outlined below. Chapter 2 has been submitted for publication in the Journal of Fluid Mechanics with myself as the primary author as well as E.W. Tedford, M. Rahmani and G.A. Lawrence as coauthors. The research program was initiated by E.W. Tedford and myself under the guidance of G.A. Lawrence. The data analysis was performed primarily by myself and E.W. Tedford with additional help from M. Rahmani. I am responsible for writing the initial manuscript with considerable edits and revi- sions from all of the coauthors, E.W. Tedford in particular. The authors of chapter 3 are myself and G.A. Lawrence. It is currently in preparation for publication. I am responsible for initiating the research and implementing the simulations and analysis. I am also responsible for writing the manuscript with editing and revisions by G.A. Lawrence. Chapter 4 has been accepted for publication, subject to revision, in the Journal of Geophysical Research, with E.W. Tedford, myself, R. Pawlowicz, R. Pieters and G.A. Lawrence as the authors. R. Pawlowicz, E.W. Tedford and myself are responsible for initiating the research. I helped with the collection and early analysis of the data. The initial manuscript was prepared by E.W. Tedford and myself with substantial revisions and refinements from all of the coauthors. The research of Chapter 5 was initiated by myself, under the guidance of G.A. Lawrence. The analysis was developed by myself in conjunction with N.J. Balmforth. I am responsible for the writing of the manuscript with editing and revisions by N.J. Balmforth and G.A. Lawrence. 1Chapter 1 Introduction and overview 1.1 Introduction In recent years, increasing attention has been paid to the processes shaping our weather and cli- mate. This has led to the widespread use of global models to predict oceanic and atmospheric circulation. However, the utility of global models relies on our ability to represent small ‘subgrid- scale’ mixing processes in terms of large-scale flow properties. Therefore, an understanding of the mechanisms in which small-scale mixing is generated is necessary in order to accurately predict large-scale transport processes. Fundamental to large-scale transport processes is the generation of instability by shear. The instability process is responsible for extracting energy from large-scale motions and transferring it to smaller scales where it will ultimately end up being converting into viscous dissipation (i.e. frictional heating) and mixing of the buoyancy field. Such a small-scale process is not merely a consequence of the large-scale dynamics, but also a cause. For example, small scale mixing is known to be responsible for the upwelling of deep water that drives the thermohaline circulation of the worlds oceans. Despite the crucial role played by small-scale mixing in density stratified environments, our present understanding of the process is incomplete (Ivey et al., 2008). The complexity of the problem is nicely illustrated by considering the evolution of the stratified shear layer. This ide- alized flow represents a simplified model of a shear-driven mixing event as it may occur in the oceans, lakes, or atmosphere. Previous studies of stratified shear layers have shown that the linear analysis originally performed by Taylor (1931) and Goldstein (1931) is sufficient to accurately capture the onset of instability (Thorpe, 1973; Pouliquen et al., 1994) as well as providing a first order description of the developing flow structures (Thorpe, 1973; Lawrence et al., 1991; Tedford et al., 2009). For these reasons, results of the linear theory are often used as a basis for developing conceptual models of mixing. Therefore, a background of the linear analysis of stratified shear flows is now briefly described. Description of linear stability analysis As a first step in building the linear model of a sheared stratified current, it is necessary to describe the equilibrium (background state) of the flow. This state must be a solution of the equations of motion ux + wz = 0 (1.1) Chapter 1. Introduction and overview 2 ut + uux + wuz = −px ρ (1.2) wt + uwx + wwz = −pz ρ − g (1.3) ρt + uρx + wρz = 0 (1.4) where (u,w) are the velocities in the horizontal and vertical directions (x, z), t is time, p represents the pressure, ρ the fluid density, g is the gravitational acceleration, and the subscripts represent partial differentiation. The first equation (1.1) expresses the law of conservation of mass in the case of an incompressible fluid. Both (1.2) and (1.3) express the conservation of horizontal and vertical momentum in the absence of viscous forces, respectively. The final equation (1.4) is a statement of the conservation of a scalar quantity in an incompressible flow that exhibits no transport of mass by molecular diffusion. However, it has been written in terms of ρ assuming that the equation of state relating them is a linear function. We now consider a parallel shear flow with u = U(z) and w = 0 that is also density stratified in the vertical so that ρ = ρ̄(z). Substituting this into (1.1 - 1.4) above demonstrates that, as long as there are no imposed pressure gradients other than hydrostatic equilibrium, any form of the horizontal velocity U(z), and density stratification ρ̄(z), are acceptable solutions. In other words, we are free to choose the U(z) and ρ̄(z) profiles as desired. The choice of U(z) and ρ̄(z) profiles is made to represent the typical background flow of interest, and is usually done with a particular application in mind. Regardless of U(z) and ρ̄(z), we now move beyond a description of the background flow to consider the evolution of perturbations from this equilibrium state. In keeping with our linear approximation, we assume that the perturbations are small superpositions on the much larger background flow, so that we can write u = U + ũ, w = w̃, p = p̄+ p̃, ρ = ρ̄+ ρ̃, (1.5) where dp̄/dz = −ρ̄g describes a hydrostatic background pressure, and the tilde quantities repre- sent the perturbation quantities. A single equation describing the evolution of the perturbations can be obtained after performing the following steps: 1. Substitute (1.5) into (1.1 - 1.4) and neglect the product of perturbation quantities (as justified by the linearity assumption). 2. Define a perturbation streamfunction ψ, such that ũ = ψz and w̃ = −ψx, and the continuity equation (1.1) is satisfied. 3. Reduce the set of equations to a single partial differential equation for ψ. 4. Take the normal mode form of the solution ψ(x, z, t) = ψ̂(z)eik(x−ct). Chapter 1. Introduction and overview 3 These steps lead to the Taylor-Goldstein (TG) equation ψ̂′′ + [ N2 (U − c)2 − U ′′ U − c − k 2 ] ψ̂ = 0, (1.6) where the primes indicate ordinary differentiation with respect to z. The TG equation describes an eigenvalue problem for the eigenvalue c and the eigenfunction ψ̂(z). While ψ̂ describes the vertical structure of the mode, it is c that determines the stability of the flow: when Im(c) ≡ ci 6= 0 then the flow is said to be unstable since perturbations grow exponentially in time at the rate kci. The perturbations may also propagate with respect to the flow with a phase speed given by Re(c) ≡ cr. The stratified shear layer The TG equation is sufficiently general to describe the linear stability of any stratified shear flow from which we can define U(z) and ρ̄(z). However, we shall focus on a U(z) that has two free- stream velocities that are connected by a smooth (shear) layer in which the majority of the shear is concentrated. Similarly, ρ̄(z) is taken as a statically stable stratification that is composed of homogeneous upper and lower layers that are separated by a diffuse density interface. The flow is also taken to be vertically unbounded so that it is free from any influence of the boundaries. It is common to represent this simple stratified shear layer by hyperbolic tangent profiles (e.g. Holmboe, 1962; Hazel, 1972; Alexakis, 2005) of the form U(z) = ∆U 2 tanh (2(z − d) h ) and ρ̄(z) = ρ0 − ∆ρ2 tanh (2z δ ) (1.7) where ∆U represents the total difference in velocity between the upper and lower streams, h gives a measure of the thickness of the shear layer, ρ0 is a reference density, and ∆ρ denotes the total change in density across the interface of thickness δ. The profiles of (1.7) also allow for a vertical offset d, between the shear layer and the density interface. From these scales we are able to identify three important dimensionless parameters defined as J ≡ g∆ρh ρ0(∆U)2 , R ≡ h δ , and a ≡ 2d h . (1.8) The bulk Richardson number J is a measure of the relative strength of the density stratification to the shear, whereas R and a represent the relative thickness of the shear layer to that of the density interface, and the vertical asymmetry between the profiles, respectively. These three dimensionless parameters (J,R, a) have important consequences for the stability of the flow, and for the subsequent development of nonlinear structures and the turbulent mixing of the density field. These consequences are now briefly described. Chapter 1. Introduction and overview 4 (d) Density field h Δρ ΔU JR (a) U(z) (b) ρ(z) (c) Ri(z) Figure 1.1: Profiles of an unstable stratified shear layer with R = 1 and a = 0. The U(z) and ρ̄(z) profiles are shown in (a,b) as well as Ri(z) in (c). For J < 1/4 the flow develops Kelvin- Helmholtz instabilities. A representative density field of the instability is shown in (d), once it has reached a large amplitude nonlinear form of development. R ≈ 1: Kelvin-Helmholtz instability In the first instance we shall take profiles with both R = 1 and a = 0 (although only R ≈ 1 is required), as shown diagrammatically in figure 1.1(a,b). In this very idealized circumstance, where the shear layer and the density interface have the same thickness and are vertically coincident, the ultimate outcome of the shear layer is determined entirely by J : when J < 1/4 the flow is subject to the Kelvin-Helmholtz (KH) instability, which leads to the formation of vorticies that eventually break down and drive turbulent mixing (figure 1.1d). When J > 1/4 the density stratification (∆ρ/h) dominates the shear (∆U/h) and the flow is stable. The onset of instability in the simple profiles of figure 1.1(a,b) once J < 1/4 is a particular result of a more general criteria that may be expressed in terms of the gradient Richardson number Ri ≡ N 2 (U ′)2 , (1.9) where N2(z) = (−g/ρ0)ρ̄′ is the squared buoyancy frequency. It was shown by Miles (1961) and Howard (1961) that if Ri > 1/4 for all z, then the flow must be stable to small amplitude (linear) disturbances. The stability of the profiles in figure 1.1(a,b) once J > 1/4 can be seen with reference to the Ri(z) plot in figure 1.1(c). Hazel (1972) has shown that as long as R < √ 2 and a = 0, the Ri(z) profile exhibits a single minimum at the centre of the layers, with a minimum value of Rimin = JR. Hence, it is necessary for J < 1/4 for instability to be possible when R = 1. Direct solution of the TG equation has found that J < 1/4 is in fact a sufficient condition for instability in this case. The very general and practical result, thatRi > 1/4 everywhere in the profile ensures stability, has come to be known as the Miles-Howard criteria. It is often used as the basis of various mixing Chapter 1. Introduction and overview 5 (d) Density field δ h Δρ ΔU JR (a) U(z) (b) ρ(z) (c) Ri(z) Figure 1.2: Profiles of an unstable stratified shear layer with R = 3 and a = 0. The U(z) and ρ̄(z) profiles are shown in (a,b) as well as Ri(z) in (c). In this case the flow is susceptible to the Holmboe mode of instability, a schematic of which is shown in the density field of (d). Arrows on the upper and lower wave crests in (d) indicate the direction of wave propagation. parameterizations which often involve a formulation of mixing that is dependent exclusively on Ri, which must be computed as a local bulk value between grid points. It is also common to utilize a critical Richardson number (sometimes taken as 1/4), above which the turbulence is suppressed by buoyancy forces (e.g. Ivey et al., 2008; Thorpe, 2005, §7.3). R > 2: The Holmboe instability The breakdown of simple mixing parameterizations such as those described above becomes ap- parent when considering the slightly more complicated, yet still highly idealized, set of symmetric (a = 0) background profiles shown in figure 1.2(a,b). In this case, the thickness of the density interface separating the homogeneous upper and lower layers is thinner than the shear layer, and has been chosen such that R = 3. Profiles that possess a relatively sharp density gradient re- gion, become susceptible to the Holmboe (1962) mode of instability, in addition to the KH. (In the case of ‘tanh’ profiles Alexakis (2005) has found that R > 2 is necessary for the presence of the Holmboe instability.) The Holmboe instability is present when J is sufficiently large (i.e. for relatively strong stratification), and consists of cusp-like propagating internal waves, as shown in figure 1.2(d). In contrast to the R = 1 profiles shown in figure 1.1(a,b), the Holmboe instability is not stabilized by increases in J , and may persist in flows where J is large (there being no upper limit on J when the flow is inviscid). This is in full agreement with the Miles-Howard criteria, as can be seen from the plots of Ri(z) in figure 1.2(c). Whereas the KH profile of Ri in figure 1.1(c) reaches a minimum in the centre of the shear layer where Rimin = JR, on the other hand, the Holmboe profile leads to a maximum at the shear layer centre with Rimax = JR, and Ri vanishing on either side whenever R > 2. This ability of the Holmboe instability to develop at large J is a result of significant shear being present above and below the density interface, and makes it a potentially relevant process in many geophysical flows. Chapter 1. Introduction and overview 6 (d) Density field h Δρ ΔU (a) U(z) (b) ρ(z) (c) Ri(z) d Figure 1.3: Profiles of an asymmetric stratified shear layer with R = 3 and a = 0.5. The U(z) and ρ̄(z) profiles are shown in (a,b) as well as Ri(z) in (c). In this case the flow is susceptible to propagating asymmetric ‘one-sided’ instabilities, a schematic of which is shown in the density field of (d). Arrows on the upper and lower waves in (d) indicate the direction of wave propagation. Asymmetric flows: ‘One-sided’ instability An additional means of generating instability at relatively large J is through a vertical displace- ment of the shear layer relative to the density interface position (figure 1.3a,b). This asymmetry, measured by the parameter a, produces lowRi on the strongly sheared side of the density interface and greater Ri on the reverse side (e.g. figure 1.3c). Instability is then found to be concentrated in the strongly sheared region of lowRi (figure 1.3d). These ‘one-sided’ flows, where instability and mixing are confined to a region either above or below the density interface, have been observed in laboratory experiments for some time (e.g. Thorpe, 1968; Koop and Browand, 1979). However, it was not confirmed that ‘one-sided’ behaviour was caused by asymmetry in the profiles until the stability analysis and laboratory experiments of Lawrence et al. (1991). Since then, Haigh (1995) and Haigh and Lawrence (1999) have performed a detailed linear stability analysis of asymmetric flows for different ranges of J , R, and a. A central question that remains is what relation these instabilities have to the better known KH and Holmboe instabilities of the symmetric (a = 0) stratified shear layer. This question will be returned to in chapter 5 of this thesis. Implications for mixing In the symmetric (a = 0) case, the presence of two different instability types when R > 2 – the KH at low J , and the Holmboe at high J – has important consequences on the mixing of mass and momentum that are not fully understood. Smyth et al. (2007) show that the eddy viscosity and diffusivity can drop by more than an order of magnitude as the flow transitions from the KH to the Holmboe instability. On the other hand, Smyth and Winters (2003) and Carpenter et al. (2007) show that under certain conditions the mixing in KH and Holmboe instabilities can be comparable. Not only the amount of mixing, but also the vertical distribution of mixing is found Chapter 1. Introduction and overview 7 to be dependent on the type of instability that results (Smyth and Winters, 2003; Carpenter et al., 2007). Mixing in KH instabilities is focused within the region of strong stratification centred on the density interface, and tends to promote the formation of an intermediate mixed layer (Caulfield and Peltier, 2000). In contrast, the Holmboe instability concentrates mixing on either side of the density interface, tending to broaden the interface (Smyth and Winters, 2003; Carpenter et al., 2007). In our analysis thus far, no mention has been made of viscous or diffusive effects. These are of course crucial to the mixing process, and are represented by the dimensionless Reynold’s and Prandtl numbers Re ≡ ∆Uh ν and Pr ≡ ν κ , (1.10) respectively. Here ν represents the kinematic viscosity and κ the diffusivity of the stratifying agent. In geophysical mixing events the stratification is usually caused by gradients in tempera- ture or salinity resulting in Prandtl numbers of Pr ≈ 9 and Pr ≈ 700, respectively. The approxi- mate range of Re that is of interest in geophysical scale applications is large, and can roughly be bounded by 103 . Re . 108. It is therefore important to note that the mixing studies discussed above utilize DNS, and are limited to low values of Re and Pr by computational constraints. Nonetheless these preliminary investigations of mixing suggest that significant changes in mix- ing behaviour may occur with only small changes in the flow parameters (i.e. J , R, or a), and illustrates that the characteristics of a shear-driven mixing event are linked to the complicated and subtle instability mechanisms of the stratified shear layer. In other words, these results suggest that the mixing of mass and momentum can depend sensitively on details of the background density and velocity structure. This is perhaps not surprising given that the stratified shear layer gives rise to a number of instability types, with stratification acting alternatively as a stabilizing and desta- bilizing influence on the shear flow (Howard and Maslowe, 1973). This dependence of mixing characteristics on the type of instability underscores the need to achieve a better understanding of the basic nonlinear behaviour of the instabilities. These results may then be used to develop more sophisticated parameterizations of mixing that account for the various instability modes that may occur. 1.2 Overview In comparison to the large body of literature devoted to the KH instability (see Peltier and Caulfield, 2003; Thorpe, 2005, §3, for recent reviews), the Holmboe instability has received relatively less attention. Previous studies have generally focused on the applicability of linear theory to observa- tions (Lawrence et al., 1991; Pouliquen et al., 1994; Hogg and Ivey, 2003; Tedford et al., 2009), the modeling of secondary circulations (Smyth and Peltier, 1991; Smyth, 2006), or preliminary investigations of the turbulence and mixing characteristics (Koop and Browand, 1979; Smyth and Chapter 1. Introduction and overview 8 Winters, 2003; Smyth et al., 2007; Carpenter et al., 2007). In particular, the evolution of the wave fields that result once the Holmboe instability reaches a finite amplitude have not previously been described, and are therefore the subject of the following chapter. This chapter also confronts the issues of comparing laboratory results to those obtained by direct numerical simulations (DNS). Since the earliest laboratory observations of the Holmboe instability (Thorpe, 1968; Koop and Browand, 1979) it has been noted that mixing occurs largely by entrainment of interfacial fluid into the upper and lower layers (see also Smyth and Winters, 2003; Carpenter et al., 2007). This has been observed to occur through a number of different processes, however, perhaps the most conspicuous of these is by the ejection of interfacial fluid from the wave crests. Although this ‘ejection process’ is generally considered to be important for turbulence and mixing (Koop and Browand, 1979; Smyth and Peltier, 1991; Smyth and Winters, 2003), as well as the evolution of the wave field (Chapter 2), no mechanism has been described to explain this process. Chapter 3, is devoted to providing such a description using the high-resolution observations from DNS. Chapter 4 is devoted to observations of shear instabilities in the highly stratified Fraser River estuary. Previous work of Geyer and Smith (1987) is extended by performing a direct compari- son of observations of shear instability in the estuary with the predictions of linear theory. Close agreement of the linear predictions with the observations demonstrate that the theory is applicable even in an evolving turbulent environment. Asymmetry between the regions of high density gra- dients and strong shear is found to be the dominant source of instability. This motivates the study of asymmetric instabilities discussed in chapter 5, as they are of relevance to the dynamics of the Fraser River estuary, and are found to be important in other highly stratified estuaries (e.g. Seim and Gregg, 1994; Yoshida et al., 1998; Sharples et al., 2003). Unlike certain viscous wall-bounded instabilities of pipe and channel flows (see Trefethen et al., 1993), the essentially inviscid instability of a stratified shear layer is well described, to first order, by linear stability theory (Thorpe, 1973; Lawrence et al., 1991; Pouliquen et al., 1994; Ted- ford et al., 2009). This provides a useful point of departure for studying the nonlinear aspects of the instabilities. This approach is adopted in chapter 5, where linear theory is used to infer charac- teristics of unstable modes that do not fit nicely into the KH/Holmboe classification. A blurring of the boundary between KH and Holmboe instabilities appears to be a general feature of stratified shear layers that possess some type of asymmetry. This asymmetry may be produced in a number of ways such as in non-Boussinesq flows (Umurhan and Heifetz, 2007), in spatially developing in- stabilities (Pawlak and Armi, 1998; Ortiz et al., 2002), or if the density interface is vertically offset from the shear layer centre (Lawrence et al., 1991; Haigh and Lawrence, 1999; Carpenter et al., 2007), to name a few. Choosing an asymmetric stratified shear layer with displaced profiles, the relationship between unstable asymmetric modes and the more familiar KH and Holmboe modes is examined. It is found that the asymmetric instabilities exhibit a gradual transition from modes of KH- to Holmboe-type as the stratification is increased (i.e. for increasing J). The analysis Chapter 1. Introduction and overview 9 is also applied to profiles measured in the Fraser River estuary to demonstrate the approach in a geophysically relevant scenario. 10 Bibliography Alexakis, A. (2005). On Holmboe’s instability for smooth shear and density profiles. Phys. Fluids, 17:084103. Carpenter, J., Lawrence, G., and Smyth, W. (2007). Evolution and mixing of asymmetric Holm- boe instabilities. J. Fluid Mech., 582:103–132. Caulfield, C. and Peltier, W. (2000). The anatomy of the mixing transition in homogenous and stratified free shear layers. J. Fluid Mech., 413:1–47. Geyer, W. and Smith, J. (1987). Shear instability in a highly stratified estuary. J. Phys. Oceanogr., 17:1668–1679. Goldstein, S. (1931). On the stability of superposed streams of fluids of different densities. Proc. R. Soc. Lond. A, 132:524–548. Haigh, S. (1995). Non-symmetric Holmboe Waves. PhD thesis, University of British Columbia. Haigh, S. and Lawrence, G. (1999). Symmetric and nonsymmetric Holmboe instabilities in an inviscid flow. Phys. Fluids, 11(6):1459–1468. Hazel, P. (1972). Numerical studies of the stability of inviscid stratified shear flows. J. Fluid Mech., 51:39–61. Hogg, A. M. and Ivey, G. (2003). The Kelvin-Helmholtz to Holmboe instability transition in stratified exchange flows. J. Fluid Mech., 477:339–362. Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geofys. Publ., 24:67–112. Howard, L. (1961). Note on a paper of John W. Miles. J. Fluid Mech., 10:509–512. Howard, L. and Maslowe, S. (1973). Stability of stratified shear flows. Boundary-Layer Met., 4:511–523. Ivey, G., Winters, K., and Koseff, J. (2008). Density stratification, turbulence, but how much mixing? Ann. Rev. Fluid Mech., 40:169–184. Koop, C. and Browand, F. (1979). Instability and turbulence in a stratified fluid with shear. J. Fluid Mech., 93:135–159. Bibliography 11 Lawrence, G., Browand, F., and Redekopp, L. (1991). The stability of a sheared density interface. Phys. Fluids, 3(10):2360–2370. Miles, J. (1961). On the stability of heterogeneous shear flows. J. Fluid Mech., 10:496–508. Ortiz, S., Chomaz, J., and Loiseleux, T. (2002). Spatial Holmboe instability. Phys. Fluids, 14:2585–2597. Pawlak, G. and Armi, L. (1998). Vortex dynamics in a spatially accelerating shear layer. J. Fluid Mech., 376:1–35. Peltier, W. and Caulfield, C. (2003). Mixing efficiency in stratified shear flows. Ann. Rev. Fluid Mech., 35:135–167. Pouliquen, O., Chomaz, J., and Huerre, P. (1994). Propagating Holmboe waves at the interface between two immiscible fluids. J. Fluid Mech., 266:277–302. Seim, H. and Gregg, M. (1994). Detailed observations of naturally occurring shear instability. J. Geophys. Res., 99:10,049–10,073. Sharples, J., Coates, M., and Sherwood, J. (2003). Quantifying turbulent mixing and oxygen fluxes in a Mediterranean-type, microtidal estuary. Ocean Dyn., 53:126–136. Smyth, W. (2006). Secondary circulations in Holmboe waves. Phys. Fluids, 18:064104. Smyth, W., Carpenter, J., and Lawrence, G. (2007). Mixing in symmetric Holmboe waves. J. Phys. Oceanogr., 37:1566–1583. Smyth, W. and Peltier, W. (1991). Instability and transition in finite-amplitude Kelvin-Helmholtz and Holmboe waves. J. Fluid Mech., 228:387–415. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Taylor, G. (1931). Effect of variation in density on the stability of superposed streams of fluid. Proc. R. Soc. Lond. A, 132:499–523. Tedford, E., Pieters, R., and Lawrence, G. (2009). Holmboe instabilities in a laboratory exchange flow. J. Fluid Mech., page submitted. Thorpe, S. (1968). A method of producing a shear flow in a stratified fluid. J. Fluid Mech., 32:693–704. Thorpe, S. (1973). Experiments on instability and turbulence in a stratified shear flow. J. Fluid Mech., 61:731–751. Bibliography 12 Thorpe, S. (2005). The Turbulent Ocean. Cambridge University Press, first edition. Trefethen, L., Trefethen, A., Reddy, S., and Driscoll, T. (1993). Hydrodynamic stability without eigenvalues. Science, 261:578–584. Umurhan, O. and Heifetz, E. (2007). Holmboe modes revisited. Phys. Fluids, 19:064102. Yoshida, S., Ohtani, M., Nishida, S., and Linden, P. (1998). Mixing processes in a highly strat- ified river. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 389–400. American Geophysical Union. 13 Chapter 2 Holmboe wave fields in simulation and experiment 1 2.1 Introduction Geophysical flows often exhibit stratified shear layers in which the region of density variation is thinner than the thickness of the shear layer (e.g. Wesson and Gregg, 1994; Yoshida et al., 1998; Sharples et al., 2003). In these circumstances, when the stratification is sufficiently strong (mea- sured by an appropriate Richardson number), Holmboe’s (1962) instability develops. At finite amplitude the instability is characterized by cusp-like internal waves (referred to herein as Holm- boe waves) that propagate with respect to the mean flow. Accurate modelling of these instabilities is important for the correct parameterization of momentum and mass transfers occurring in flows of this nature. Previous studies on the nonlinear behaviour of Holmboe waves have adopted one of two meth- ods: either an experimental approach in which the instability is studied under specified laboratory settings (Pouliquen et al., 1994; Zhu and Lawrence, 2001; Hogg and Ivey, 2003), or a numerical approach that allows for a detailed description of the flow in an idealized stratified mixing layer (Smyth et al., 1988; Smyth and Winters, 2003; Smyth, 2006; Smyth et al., 2007). It is difficult to make a meaningful comparison of laboratory and numerical results for a number of reasons. In the case of laboratory experiments, Holmboe waves often arise as a local feature of a larger-scale flow, such as an exchange flow between two basins of different density (Pawlak and Armi, 1996; Zhu and Lawrence, 2001; Hogg and Ivey, 2003), or an arrested salt wedge flow (Sargent and Jirka, 1987; Yonemitsu et al., 1996). In many of these experiments the mean flow varies appreciably over length scales that are comparable to the wavelength of the waves. For this reason, it can be difficult to isolate the dynamics of the waves from that of the mean flow. Investigations utilizing numerical simulations comprises a great majority of the literature on the nonlinear dynamics of Holmboe waves. The first verification of two oppositely propagating cusp-like waves of equal amplitude, predicted by the Holmboe (1962) theory, was made through the numerical simulations of Smyth et al. (1988). Since then, increases in computational resources have led to fully three-dimensional direct numerical simulations (DNS) of Holmboe waves that resolve the smallest scales of variability. These simulations have been used to understand turbu- lence and mixing characteristics (Smyth and Winters, 2003; Smyth et al., 2007; Carpenter et al., 2007), as well as the growth of secondary circulations and the transition to turbulence (Smyth, 1A version of this chapter has been submitted for publication. J.R. Carpenter, E.W. Tedford, M. Rahmani and G.A. Lawrence (2009) Holmboe wave fields in simulation and experiment. Chapter 2. Holmboe wave fields in simulation and experiment 14 2006). However, partly due to computational constraints, only a single wavelength of the primary instability has been reported in the literature. Furthermore, no attempt has been made to compare the results of numerical simulations with laboratory experiments. In this chapter, we undertake a combined numerical and experimental study of Holmboe waves. The experiments, originally described by Tedford et al. (2009), consist of an exchange flow through a relatively long channel with a rectangular cross-section. The experimental design allows for a detailed study of the Holmboe wave field within a steady mean flow that exhibits gradual spatial variation relative to the wave properties. The DNS of the present study were de- signed to correspond as closely as possible to the conditions present in the experiments to effect a meaningful comparison between the two methods. To our knowledge, this is the first study to compare experimental and numerical results, as well as the first to perform DNS for multiple wavelengths of the instability. We focus on comparing basic descriptors of the wave fields such as phase speed, wavenumber, and wave amplitude, in order to gain a fuller understanding of the processes affecting the nonlinear behaviour of the waves. The chapter is organized as follows. Section 2.2 gives a background on the stability of strat- ified shear flows. This is followed by a description of the numerical simulations, and laboratory experiments in §2.3. We then discuss comparisons between the simulations and experiments in terms of the basic wave structure (§2.4), phase speed (§2.5), wave spectral evolution (§2.6), and wave amplitude and growth (§2.7). Conclusions are stated in the final section. 2.2 Linear stability of stratified shear layers In both experiment and simulation, the mean flow exhibits the characteristics of a classic stratified shear layer. The velocity profile undergoes a total change of ∆U , over a length scale h, that is closely centred with respect to the density interface. Similarly, the density profile changes by ∆ρ between the two layers, over a scale of δ. This suggests using an idealized model of the horizontal velocity and density profiles that is given by U(z) = ∆U 2 tanh [2(z − z0) h ] and ρ̄(z) = ρ0 − ∆ρ2 tanh [2(z − z0) δ ] , (2.1) respectively. The density profile ρ̄(z) is measured relative to a reference density ρ0, with z the vertical coordinate. A necessary condition for the growth of Holmboe’s instability is that the thickness ratio R ≡ h/δ > 2 (Alexakis, 2005). In addition to R, we may define three more important dimensionless parameters Re ≡ ∆Uh ν , J ≡ g ′h (∆U)2 , and Pr ≡ ν κ , where g′ = ∆ρg/ρ0 is the reduced gravitational acceleration, ν is the kinematic viscosity, and Chapter 2. Holmboe wave fields in simulation and experiment 15 0 0.5 1 1.5 0 0.5 1 1.5 c r ( cm s -1 ) (a) 0 0.5 1 1.5 0 0.4 0.8 1.2 (b) 0 0.5 1 1.5 0 0.05 0.1 0.15 (c) σ ( ra d s -1 ) k c i (r ad s -1 ) k (rad cm-1)k (rad cm-1)k (rad cm-1) Figure 2.1: Plots of (a) phase speed cr, (b) frequency σ, and (c) growth rate kci, of the profiles in (2.1). Conditions in the experiment are shown as thick lines, and the three-dimensional simulation at Pr = 25 and R = 5 as thin lines. No noticeable difference between the simulation and experiment can be seen in (b). The locations of the wavenumber of maximum growth are marked with a vertical dotted line in (c). The dashed line in (b) indicates σ ∝ k. κ the diffusivity of the stratifying agent. These are the Reynolds, bulk Richardson, and Prandtl numbers, respectively. Linear stability analysis of the profiles in (2.1) has been performed in numerous studies (e.g. Hazel, 1972; Smyth et al., 1988; Haigh, 1995). For the flows considered here, the effects of viscosity and mass diffusion have been included. The resulting equation is a sixth order eigenvalue problem originally described by Koppel (1964). Like the better known Taylor–Goldstein equation, Koppel’s equation gives predictions of the complex phase speed c = cr + ici, and vertical mode shape, as a function of wavenumber k. Results of the stability analysis are shown in figure 2.1, which includes the dispersion relation in terms of phase speed cr(k), and frequency σ(k), as well as the temporal growth rate, kci. This is done for the idealized profiles (2.1) using Re = 630, J = 0.30, Pr = 700, and R = 8, matching the dimensionless parameters in the laboratory exchange flow. As discussed in the next section, computational constraints limited our three- dimensional DNS to a Pr = 25 and R = 5, resulting in slightly different results (figure 2.1). Although no appreciable changes are seen in the predicted phase speed cr and frequency σ, there are differences in maximum growth rate and the location of the wavenumber of maximum growth kmax. 2.3 Methods 2.3.1 Description of the numerical simulations Numerical simulations were performed using the DNS code described by Winters et al. (2004), which has been modified to include greater resolution of the density scalar field by Smyth et al. Chapter 2. Holmboe wave fields in simulation and experiment 16 -2 0 2 0 2 4 6 8 10 U (cm s-1) z ( cm ) -1 0 1 ρ−ρ0 (kg m-3) (a) (b) Experiment t = 0 s t = 100 s t = 195 s Simulation: Figure 2.2: Evolution of the background profiles in both simulation and experiment. Plots (a) and (b) show temporal changes to U(z) and ρ̄(z) − ρ0 profiles in the simulation with experimental profiles taken from the channel centre (x = 0) in thick lines overlain for comparison. (2005). The simulations were designed to reproduce conditions present in the laboratory exper- iment as closely as possible, while still conforming to the general methodology used in recent investigations of nonlinear Holmboe waves (Smyth and Winters, 2003; Smyth, 2006; Smyth et al., 2007; Carpenter et al., 2007). The boundary conditions are periodic on the streamwise (x) and transverse (y) boundaries, and free-slip on the vertical (z) boundaries. Simulations are initialized with profiles in the form of (2.1) that closely match what is observed in the experiment. Figure 2.2 shows a sequence of representative U and ρ̄−ρ0 profiles at three different times during a simulation, as well as profiles from the experiment for comparison. The periodic boundary conditions of the simulations cause the flow to ‘run down’ over time, i.e. there is a continual loss of kinetic energy from the shear layer due to viscous dissipation and mixing. This results in an increase of h and δ over time, as can be seen in figure 2.2. To indicate conditions at the initial time step (t = 0 s) of the simulations we will use a zero subscript (e.g. h0). In order to initiate growth of the primary Holmboe instability, the flow is perturbed with a random velocity field at the first time step. The noise is distributed evenly in the x, y directions, but given greater amplitude near the centre of the shear layer and density interface, in the same manner as Smyth and Winters (2003). The amplitude of the random perturbation was chosen large enough such that the instability grows to finite amplitude with minimal diffusion of the background profiles, yet is still small enough to satisfy the conditions for numerical stability given our choice of time step. While an ideal comparison between simulation and experiment would involve matching all Chapter 2. Holmboe wave fields in simulation and experiment 17 Parameters Linear Theory Results Pr R Lx Ly kmax cr arms (cm) (cm) (rad cm−1) (cm s−1) (cm) Experiment 700 8 200 10 0.91 0.79 0.31 (0.51-0.84) I Simulation (3D) 25 5 128 5 0.79 0.84 0.62 II Simulation (2D) 700 8 64 0 0.91 0.79 0.48 III Simulation (3D) 25 5 64 10 0.79 0.84 0.62 IV Simulation (2D) 25 5 64 0 0.79 0.84 0.74 Table 2.1: Values of various important parameters for the simulations and experiment. The pa- rameters listed in the simulations are evaluated using the initial conditions. In all cases we have J0 = 0.3, Re0 = 630, and Lz = 10.8 cm. Also included are kmax and cr from the results of the linear stability analysis, and the root mean square saturated amplitude observations. The values in parentheses under the experimental arms column indicate the predicted range of wave amplitudes once wave stretching is accounted for (see §2.7). four of the relevant dimensionless parameters, we are constrained by the high computational de- mands of DNS. Of particular difficulty is the fine grid resolution required for high Pr flows. For this reason we have chosen a Pr = 25, opposed to Pr = 700 for the laboratory salt stratification. Large values of R also place a high demand on the computational resources, and we have there- fore chosen R0 = 5, opposed to the R = 8 observed in the experiments. The effects of Pr and R have been tested by performing a two-dimensional simulation at Pr = 700 and R0 = 8. The remaining two parameters, Re0 = 630 and J0 = 0.30, have been matched to the experimental values. Computational constraints also limit the size of the simulation domain. In all cases the vertical depth Lz , has been matched to the 10.8 cm of the experiments. The simulation width Ly = 5 cm, has been reduced to half of that in the experiment (Ly = 10 cm), but was not found to effect the results presented. This reduction in the width of the computational domain enabled a larger length Lx = 128 cm, allowing for approximately 16 wavelengths of the most amplified mode, and compares well with the Lx = 200 cm in the experiment. A summary of the parameters in the experiment and the simulations is shown in table 2.1. In addition to the three- and two-dimensional simulations already mentioned (labeled I and II in table 2.1, respectively), two supplementary simulations (III and IV) were also performed to test the effects of Ly, R and Pr. Unless explicitly stated, we will refer to simulation I simply as ‘the simulation’, hereafter. 2.3.2 Description of the laboratory experiment The laboratory experiment was performed in the exchange flow facility described by Tedford et al. (2009). A complete discussion of the experimental procedures and apparatus can be found in that study, however, we now provide a summary of the pertinent features. Chapter 2. Holmboe wave fields in simulation and experiment 18 0 5 10 z ( cm ) (a) -100 -50 0 50 100 -1 -0.5 0 0.5 1 x (cm) U ( cm s -1 ) (b) Experiment Simulation Linear fit ρ1 ρ2>ρ1 U 1 U 2 Figure 2.3: Spatial changes in U(z) and layer depths that occur along the laboratory channel are shown in (a), along with the corresponding distribution of Ū(x) in (b). A linear fit to Ū(x) in the central portion of the channel is shown as the dashed line, and the mean velocity in the simulation domain is given by the thin solid line. The apparatus consists of two reservoirs connected by a rectangular channel 200 cm in length, and 10 cm in width. The reservoirs are initially filled with fresh and saline water (∆ρ = 1.41 kg m−3) such that the depth in the channel is 10.8 cm. A bi-directional exchange flow is initiated by the removal of a gate from the centre of the channel. After an initial transient period in which gravity currents propagate to each reservoir, and mixed interfacial fluid is advected from the chan- nel, the flow enters a period of steady exchange where the density interface is found to display an abundance of Holmboe wave activity. In contrast to the run-down conditions in the DNS, the storage of unmixed water in the reservoirs maintains a steady exchange flow for approximately 600 s. Our comparison is restricted to instabilities observed during the period of steady exchange. The exchange flow exhibits internal hydraulic controls at the entrance to each of the reser- voirs, effectively isolating the channel from disturbances in the reservoirs, and enforcing radiation boundary conditions at the channel ends. Friction between the layers leads to a gradually sloping density interface (figure 2.3) that produces an x-dependent mean velocity, Ū = (U1 +U2)/2, with the upper (U1) and lower (U2) layer velocities defined as the maximum and minimum free-stream velocities (Tedford et al., 2009). The gradual variation of Ū(x) along the laboratory channel is a result of the acceleration in each of the layers due to the sloping interface. This variation is shown in figure 2.3(b), and is found to be a near-linear function of x for the central portion of the channel. In contrast, Ū is identically zero throughout the domain in the simulation, due to the periodic boundary conditions. This difference in mean flow is found to have important effects on Chapter 2. Holmboe wave fields in simulation and experiment 19 z ( cm ) 0 5 10 z ( cm ) 0 5 10 z ( cm ) 0 5 10 z ( cm ) x (cm) -20 -10 0 10 20 0 5 10 (a) (b) (c) (d) Figure 2.4: Representative plots of the density field for the experiment (a) along with simulation II (b), and simulation I (c,d). The plot in (d) is taken at a later time when two ejections are underway, indicated by arrows. The x-axis has been shifted by Lx/2 to the left in (d) to better display the ejection process. the nonlinear development of the Holmboe wave field. 2.4 Wave structure In the first instance, it is beneficial to perform a simple visual comparison of the density struc- ture of the waves. This is shown in figure 2.4, where a representative photograph of the labora- tory waves is displayed above plots of the density field from the two-dimensional simulation II (figure 2.4b) and three-dimensional simulation I (figure 2.4c,d). The density structure in figures 2.4(a,b) is very similar, as each has an identical set of dimensionless parameters, differing only in the initial and boundary conditions. In all panels of figure 2.4 it can be seen that although many of the waves display the typical form of the Holmboe instability, and consist of cusps projecting into Chapter 2. Holmboe wave fields in simulation and experiment 20 the upper and lower layers, others appear sinusoidal in form. The waves focused above the density interface (upwards pointing cusps) are moving from left to right, in the same direction as the flow in the upper layer, while the waves in the lower layer (downwards cusps) move at an equal but opposite speed with respect to the mean velocity. An important feature of nonlinear Holmboe waves is the occasional ejection of stratified fluid from the wave crests into the upper and lower layers. Two such ejections are shown in figure 2.4(d) where indicated, and can be characterized by thin wisps of fluid being drawn from the wave crest and advected by the mean flow. These wisps often settle back to the interface level, contributing to the accumulation of mixed fluid there. This accumulation is observed to a much greater extent in figure 2.4(c,d), and should be expected due to the larger value of R0, as well as the higher diffusion that comes with the lower Pr used in this simulation. Although ejections are observed in both the laboratory experiment and high Pr and R simulation (II), there is a greater frequency of occurrence in the lower Pr = 25 and R = 5 simulation (I). Holmboe’s instability has the uncommon property that, under certain conditions, the growth of the primary instability may take place as a three-dimensional wave. Such a wave would travel obliquely to the orientation of the shear, and produce significant departures from a two- dimensional wave. One of the conditions for this three-dimensional growth is that Re be suffi- ciently low (Smyth and Peltier, 1990). As the laboratory experiments are carried out at low Re, and show some variation in the spanwise direction, it must be questioned whether the growth of the primary Holmboe instability is three-dimensional. This is easily tested by the simulation results, which show a clear two-dimensional growth (see §2.7 as well), even to an initially ran- dom perturbation as described above. We can therefore confirm that the primary instability is two-dimensional for the conditions examined in the present study. 2.5 Phase speed Many of the basic features in the wave field are revealed by an x− t characteristics diagram of the density interface elevation, shown in figure 2.5 for both the simulation and experiment. Although the interface consists of contributions from both upper and lower Holmboe wave modes (each travelling in opposite directions), we have filtered the characteristics using a two-dimensional Fourier transform to reveal only the upper, rightward propagating wave modes. Certain differences between the simulation characteristics (figure 2.5a) and the experimental characteristics (figure 2.5b) are immediately apparent. The experimental characteristics exhibit a greater degree of irregularity. Since each plot represents a two dimensional slice from a three- dimensional field, this may be a result of greater variability in the transverse direction in the case of the experiments. Since the waves in the simulation develop from an initial random perturbation at t = 0, there is also a temporal growth of the average wave amplitude in figure 2.5(a) that is not present in the experimental characteristics. Chapter 2. Holmboe wave fields in simulation and experiment 21 x (cm) t (s ) (a) Simulation 0 20 40 60 80 100 120 140 160 180 x (cm) (b) Experiment -20 0 20 Figure 2.5: Rightward propagating wave characteristics for the simulation (a) and experiment (b). Shading represents the elevation of the density interface with red indicating a high (crest) and blue indicating a low (trough), and has a different scale in each of (a,b). Solid black lines indicate the characteristic slope given by the linear prediction of phase speed cr. In the case of the laboratory experiment, the cr has a slight curvature since the changes in Ū along the channel have been included. The dark circles indicate locations and times of ejections. Despite these apparent differences in the characteristics, the phase speeds (inferred from the slope of the characteristics) are in good agreement. The observed phase speeds in both the simu- lation and experiment are found to be slightly greater than the predictions of linear theory (solid lines), which has been noted in previous studies (Haigh, 1995; Hogg and Ivey, 2003; Tedford et al., 2009). However, the observations also suggest an increase in phase speed with wave am- plitude. This is a quintessential feature of nonlinear wave behaviour (e.g. Stokes waves). Note that a ‘pulsing’ of the wave amplitude and phase speed is present in both sets of characteristics in figure 2.5. This is a well known feature of Holmboe waves due to the interaction between the two oppositely propagating modes (Smyth et al., 1988; Zhu and Lawrence, 2001; Hogg and Ivey, 2003). Sudden decreases in wave amplitude can be seen in both sets of characteristics at a number of times and locations. It is often the case (though not always) that these sudden amplitude changes are a result of the ejection process. Instances where ejections occur have been identified in fig- ure 2.5, and are denoted by circles. It is generally observed that the ejection process preferentially acts on the largest amplitude waves, similar to the breaking of surface waves. Chapter 2. Holmboe wave fields in simulation and experiment 22 x (cm) t (s ) (a) Simulation -50 0 50 0 50 100 150 x (cm) (b) Experiment -50 0 50 Figure 2.6: Rightward propagating wave characteristics for the simulation (a) and experiment (b). White indicates a wave crest while grey indicates a wave trough. In each panel a number of wave crests are indicated by solid and dashed lines. In (a), the dashed lines correspond to waves that are ‘lost’ over the duration of the simulation, whereas in (b), the dashed lines correspond to waves that have formed within the channel. Only the central portion of the laboratory channel corresponding to the simulation domain has been shown. Circles and squares indicate locations and times of ejections and pairing events, respectively. 2.6 Spectral evolution This section concerns the distribution and evolution of wave energy with k. It will be shown that there are two different processes acting separately in the simulation and experiment that are responsible for a shifting of wave energy to lower k (i.e. longer waves). 2.6.1 Frequency shifting In order to gain an understanding of the wave spectrum, it is first useful to carefully examine the characteristic diagrams. Figure 2.6 shows rightward propagating characteristics from both simulation and experiment that highlight the location of wave crests (in white) and troughs (in grey). Characteristics from the simulation (figure 2.6a) are discussed first, and are shown for the entire computational domain. Beginning with the initial random perturbation at t = 0, energy is extracted from the mean Chapter 2. Holmboe wave fields in simulation and experiment 23 x (cm) z ( cm ) -40 -35 -30 -25 0 5 10 x (cm) -30 -25 -20 x (cm) -20 -15 -10 (a) (b) (c) -0.5 0 0.5 1 1.5 Figure 2.7: Vorticity field from the simulation illustrating the vortex pairing process. Colour bar indicates the magnitude of the spanwise (y) vorticity field nondimensionalized by ∆U/h with positive values corresponding to clockwise rotation. Times shown are (a) t = 80 s, (b) t = 87 s, and (c) t = 97 s. Black contours denote isopycnals with a spacing of ∆ρ/8, and arrows indicate the position of vorticies in the upper layer that undergo pairing. Note that two vorticies are pairing in the lower layer in (c) as well. flow by the instability and fed into the wave field at, or very close to, the wavenumber of maximum growth, kmax. This results in approximately 16 waves in the computational domain (given by Lxkmax/2pi) for early times. We see however, that as the simulation proceeds wave crests are continually being ‘lost’ over time. This feature is highlighted by the solid and dashed lines that are used to trace the wave crests in figure 2.6(a). The dashed lines indicate wave crests that are ‘lost’, while the solid lines correspond to crests that persist. This process of losing waves results in an observed frequency, ω, that is continually shifted downwards. Because previous numerical studies of Holmboe waves simulated only a single wavelength, this process has not been described before. This ‘frequency downshifting’ or ‘wave coarsening’ has, however, been noted previously in other nonlinear wave systems (e.g. Huang et al., 1999; Balmforth and Mandre, 2004). In general, for all of the simulations performed, the ejection process typically results in a loss of waves and a downshift in frequency (figure 2.6a). This observation mirrors similar findings in the frequency downshifting of nonlinear surface gravity waves, where the occurrence of wave breaking is related to lost waves (Huang et al., 1996; Tulin and Waseda, 1999). Close examination of the vorticity field also suggests that the Holmboe waves undergo a vortex pairing process. Each interfacial wave is led by a vortex formed from the ‘rolling up’ of the initial shear layer vorticity. This leading vortex propagates with the waves in the upper and lower layers, and can be seen in the vorticity field plotted in figure 2.7(a). The pairing process is seen to act on the two vorticies (indicated by arrows) in the upper layer of figure 2.7; they draw closer together, rotating about each other (figure 2.7b) and eventually merge into a single vortex (figure 2.7c). Although the pairing of adjacent vorticies is a well known feature of homogenous and weakly stratified shear layers (Browand and Winant, 1973), it has received little attention in Holmboe waves, the only observation being that of Lawrence et al. (1991). This is an additional means to Chapter 2. Holmboe wave fields in simulation and experiment 24 effect a shift of wave frequency, and is denoted by square symbols in figure 2.6(a). In contrast, figure 2.6(b) shows that the experimental characteristics display a distinctly dif- ferent behaviour. In this case, new wave crests are continually being formed as the waves traverse the channel. Again, this process is highlighted by the tracing of crests by solid and dashed lines. Now, the dashed lines represent new wave crests that have been formed within the channel. This process results in an increasing ω with x in the experiments. Tedford et al. (2009) explain the formation of new waves as follows. As waves propagate through the channel they are accelerated by the increasing mean velocity Ū(x). This leads to a ‘stretching’ of the waves that decreases k from near kmax, where the waves initially formed, to lower values (i.e. longer wavelengths). Once a sufficiently low k is achieved, the Holmboe instability mechanism acts between the wave crests to form additional waves. This feeds energy back into the wave field near kmax, resulting in an average k that is constant across the channel, and an increasing ω. The two processes, wave coarsening in the simulation, and wave stretching in the experiments, can be described quantitatively using wave spectra. 2.6.2 Wave energy spectra Differences between the processes responsible for modifying k in both the simulation and experi- ment can be seen in figure 2.8. It demonstrates how wave energy (indicated by the dark bands) is redistributed in k over time. The spectra of the simulation (figure 2.8a), which has been normalized by the variance in order to remove the time-dependent growth of the waves, shows a discrete transfer of wave energy to lower k. The simulation spectra is required to evolve in discrete steps due to the periodic boundary conditions (i.e. in wavenumber increments of ∆k = 2pi/Lx). As a point of comparison, the kmax prediction from linear stability theory is plotted in red. The predicted kmax has been discretized according to the boundary conditions, and decreases in time due to the diffusion of the background profiles, i.e. the increase in the shear layer thickness h(t). The spectral evolution plot (figure 2.8a) compliments the characteristics diagram of figure 2.6(a), showing an initial input of energy at kmax(t = 0), and a subsequent shifting of that energy to lower k. It is interesting to note that the kmax(t) curve shows the same general trend as the concentration of wave energy (shown by the dark ‘blocks’ in figure 2.8a). Although the details are unclear, we speculate that the shift in wave energy to lower k is the result of nonlinear processes such as the ejections and vortex pairing. It is apparent from the wave spectra in figure 2.8(b) that the process responsible for the redis- tribution of wave energy in the experiments is a continuous one. Energy at any given time is found to be focused in a number of ‘bands’. These bands originate near kmax, and move towards lower k in time. In addition, they all appear to have a similar trajectory in kt-space. Tedford et al. (2009) Chapter 2. Holmboe wave fields in simulation and experiment 25 0 t (s ) (a) Simulation 0 0.5 1 1.5 0 100 200 k (rad cm-1) (b) Experiment 0 0.5 1 1.5 0 100 200 300 400 500 600 stretching prediction k max prediction k (rad cm-1) 10 -1 10 cm 2 cm 2 Figure 2.8: Spectral evolution of the rightward propagating waves from simulation (a) and ex- periment (b). Dark colours denote a high in energy which is proportional to the mean square amplitude of the interface displacement. The wave energy has been normalized by the variance in (a) to remove the time dependent wave growth. Linear stability theory is used to predict kmax (red lines), which changes in time for the simulations. The predicted stretching of wave energy in the experiment by Ū(x) to lower k is shown as the yellow dashed lines in (b). The beginning of the steady period of exchange is referenced to t = 0 s in (b). hypothesize that these bands are a result of the stretching of wave energy to lower k by Ū(x). We now formulate a simple model in order to quantify this hypothesis. Wave stretching prediction The changes in k that result from wave stretching by Ū(x) can be described by an application of gradually varying wave theory. This theory assumes that the density interface elevation η(x, t), may be expressed in terms of a gradually varying amplitude a(x, t), and a rapidly varying sinu- soidal component viz. η(x, t) = Re{a(x, t)eiθ(x,t)}. (2.2) The local wavenumber and frequency are defined in terms of the phase function θ(x, t) by k ≡ ∂θ/∂x and ω ≡ −∂θ/∂t, respectively. We assume, for the moment, that θ(x, t) is a smooth Chapter 2. Holmboe wave fields in simulation and experiment 26 function. This implies that waves are conserved, giving ∂k ∂t + ∂ω ∂x = 0. (2.3) Recognizing that ω, which is the frequency that a stationary observer would measure, includes both an intrinsic portion σ(k), and an advective portion kŪ , leads to ω = σ(k) + kŪ(x). (2.4) Substituting into (2.3) gives Dk Dt = −Sk, (2.5) where the material derivative, defined as D Dt ≡ ∂ ∂t + (cg + Ū) ∂ ∂x , denotes changes in time while moving at the speed cg + Ū , and cg ≡ dσ/dk is the intrinsic group speed. This is the speed that wave energy, i.e. the dark bands in figure 2.8(b), is expected to propagate through the channel. We have also defined S ≡ dŪ/dx, which is found to be very nearly constant in the central portion of the laboratory channel (see figure 2.3b). Choosing a Lagrangian frame of reference, that moves at the speed cg + Ū through the channel, allows for a simple integration of (2.5) to give k(t) = k∗e−S(t−t∗), (2.6) where k∗ = k(t∗) is some initial value of k that wave energy begins the stretching process at. A direct comparison is now possible between the prediction of (2.6) and the bands of energy in the observed spectral evolution. The prediction is shown by the yellow dashed lines in figure 2.8(b), and is found to be in agreement with the observations. This validates the hypothesis that the spectral shift towards lower k is a result of wave stretching. The excellent agreement between the predictions and observations also reveals that our assumption of wave conservation is justi- fied. This is not in contradiction with the formation of new waves described in §2.6.1 since wave conservation is applied only after energy is fed into the wave field by the instability mechanism. 2.7 Wave growth and amplitude The final basic parameter that we intend to compare is the wave amplitude, a. This feature of the wave field is determined when the linear growth reaches some level where it saturates. It is a nonlinear property of the waves, and may involve three-dimensional effects as well as interaction Chapter 2. Holmboe wave fields in simulation and experiment 27 with the mean flow. We begin by discussing the various phases of wave growth. 2.7.1 Wave growth In the simulation, the instability mechanism causes the growth of waves from an initial random perturbation into a large-amplitude nonlinear wave form. This growth process is best illustrated by considering the kinetic energy of the waves, K. Following Caulfield and Peltier (2000), we partitionK into a two-dimensional kinetic energyK2d associated with the primary Holmboe wave, and a three-dimensional component K3d, that provides a measure of the departures from a strictly two-dimensional wave. By this partitioning we have K = K2d +K3d, (2.7) where K2d = 〈u2d · u2d/2K0〉xz and K3d = 〈u3d · u3d/2K0〉xyz, (2.8) and we have used u1d(z, t) = 〈u〉xy, (2.9) u2d(x, z, t) = 〈u− u1d〉y, (2.10) u3d(x, y, z, t) = u− u1d − u2d, (2.11) with 〈·〉i representing an average in the direction i, and K0 the total kinetic energy at t = 0. The K2d and K3d components are plotted on a log-scale in figure 2.9 for the simulation. The plot indicates that after a start up period where the energy of the initial perturbation rapidly de- cays, a stage of exponential growth is achieved in K2d. This stage of exponential growth can be compared to the prediction of linear theory, and is found to be less than the prediction (kci = 0.10 s−1 versus the observed kci = 0.06 s−1). One possible explanation of this discrepancy in kci is due to the fairly rapid spreading rate of the shear layer, measured by r = h−1dh/dt, during the stage of linear growth up to t ≈ 65 s. This spreading is due to the viscous diffusion of vorticity in the shear layer, and may be compared with the growth rate of the instability by the ratio r/kci. Linear theory implicitly assumes that r/kci  1 indicating that there is no time dependence of the background profile. While r/kci ≈ 0.1 is relatively small during the linear growth period, it may be sufficient to explain the differences in observed growth rates. The wave growth is entirely two-dimensional until the waves have reached a finite amplitude (t ≈ 65 s), at which point the growth of three-dimensional secondary structures results (see Smyth 2006 for a discussion of this process in Holmboe waves). However, the waves remain primarily two-dimensional, with K3d at least an order of magnitude smaller than K2d. There is not a well defined transition to turbulence, as is found in other studies of stratified shear layers (e.g. Caulfield Chapter 2. Holmboe wave fields in simulation and experiment 28 0 50 100 150 200 10-6 10-5 10-4 10-3 10-2 10-1 K in et ic E n er g y t (s) K 2d K 3d linear growth Figure 2.9: Growth of K2d and K3d for the simulation. The thick line gives the linear growth rate prediction of the growth of K2d, which is a function of time due to the changing background profiles. and Peltier, 2000; Smyth et al., 2001), likely due to the low Re. Once the saturated amplitude is reached, there is a slow decline of K over the remainder of the simulation. In the laboratory experiments we have focused only on the period of steady exchange, and therefore do not observe a time-dependent growth of the wave field on average. However, as discussed previously, the instability is constantly acting to produce new waves along the channel. It is difficult to measure the growth rate of these waves, but they appear to reach a saturated amplitude rapidly, suggesting that they are strongly forced by disturbances within the channel (Tedford et al., 2009). 2.7.2 Comparison of saturated amplitudes Although the transient growth of the instability is difficult to quantify in the experiments, it is possible to measure the mean amplitude of the waves. This is done by using the root mean square amplitude of the interface elevation η(x, t), given by arms(x) = √ 1 T ∫ T η2 dt, (2.12) where T denotes the duration of the steady period of exchange. When averaged over a number of experiments arms is found to display little dependence on x. A similar arms can be defined for the simulations, however, the temporal average is replaced by a spatial average in x. The Chapter 2. Holmboe wave fields in simulation and experiment 29 growth of arms in time in the simulations shows a similar behaviour to K2d; an exponential initial growth, followed by a saturation, and subsequent decay. The saturated (maximum) amplitude reached during each of the simulations is shown in table 2.1, along with the mean amplitude in the experiments. The first feature to note is that the waves of the two-dimensional simulation (II) at R = 8 and Pr = 700 (matching the conditions in the experiment) have a lower amplitude than of all the other simulations, especially the two-dimensional simulation (IV) at R = 5, Pr = 25. This indicates that there is a possible dependence of the saturated amplitude on R,Pr. Most impor- tantly, the amplitude measured in the experiments is significantly smaller than any of the saturated amplitudes reached in the simulations. The small amplitudes observed in the experiments can be explained by, once again, appealing to the effects of wave stretching. Wave stretching effects on amplitude To understand the effects of wave stretching on amplitude in the experiments, we apply principles that have been established for waves on slowly varying currents (e.g. Peregrine, 1976). In doing so, we assume that the Holmboe waves may be represented by a simple train of linear internal waves that satisfy the dispersion relation in figure 2.1(c). We are then able to track the changes in wave amplitude that occur as a result of the spatially varying mean velocity Ū(x), i.e. the wave stretching. In this simplified model it is the conservation of wave action density that is relevant. This is given as E/σ, where E is the wave energy density, and recall that σ(k) is the intrinsic wave frequency. Substitution into the conservation law, and following a similar procedure to §2.6.2 leads to a similar result D Dt (E σ ) = −S (E σ ) , (2.13) which describes changes in action density due to the stretching by Ū . In arriving at (2.13) we have neglected a term that is proportional to d2σ/dk2, which is small in the range of k that we are interested in (see figure 2.1b). Taking S to be constant once again, allows for simple integration of (2.13) to give (E σ ) = (E σ ) ∗ e−S(t−t∗). For linear internal waves E ∝ a2, so that we have an estimate of the reduction in wave amplitude due to stretching of a a∗ = √ σ σ∗ e−S(t−t∗)/2. (2.14) If we now take the intrinsic frequency σ ∝ k, as suggested by the linear dispersion relation in figure 2.1(b), it is possible to write the right hand side of (2.14) as e−S(t−t∗), where we have used the spectral prediction in (2.6). By inspection of figure 2.8(b), we can estimate a time interval, ∆t, that wave energy spends in the channel (i.e. the average time interval that the dark bands Chapter 2. Holmboe wave fields in simulation and experiment 30 appear for) to be between 100 and 200 s. The amplitude reduction is therefore in the range 0.37 < e−S∆t < 0.61. Given this reduction, and assuming that no other processes are taking place that may affect the wave amplitudes, we would expect amplitudes in the range of that shown in table 2.1, given in parentheses. This adjusted amplitude is comparable to the observations in the simulations, and demonstrates that – in the absence of other processes – the stretching of wave action is significant in reducing the experimental wave amplitudes. 2.8 Conclusions We have compared simulations of Holmboe wave fields with the results of laboratory experiments. A meaningful comparison was possible since both methods exhibit only gradual variations in the mean flow. In the laboratory experiment, the mean flow is spatially varying, whereas the numerical simulations display a temporal variation. Focusing on basic descriptors of the waves, such as phase speed, wavenumber, and amplitude, we have identified a number of processes affecting the nonlinear behaviour of Holmboe wave fields. Similarities between results of the two methods include the basic structure of the waves, and the phase speeds. The observations show slightly greater phase speeds when compared with the predictions of linear theory, in agreement with previous studies (Haigh, 1995; Hogg and Ivey, 2003; Tedford et al., 2009). Further departures from the linear predictions are attributed to a nonlinear dependence of the phase speed on amplitude. The greatest differences between simulation and experiment are found in the spectral evolution and wave amplitudes. In simulations, a transfer of wave energy to lower k was found to result from wave coarsening, which caused waves to be ‘lost’ through discrete merging events. These events were typically found to result from both vortex pairing as well as the ejection process. The latter of which is suggested to be similar to wave breaking in surface waves, since it appears to act preferentially to reduce the amplitude of the largest (steepest) waves. A detailed investigation of both ejections and vortex pairing in Holmboe waves is currently underway. The shift of wave energy to lower k that was observed in the experiments can be attributed to the ‘stretching’ of the wave field by the spatially accelerating mean flow. This suggestion of Tedford et al. (2009) has been confirmed by a simple application of gradually varying wave theory, which is able to accurately predict the time dependence of the spectral shift. The wave stretching process is also expected to have a significant effect in reducing the wave amplitudes observed in the experiments. This conclusion appears sufficient to explain discrepan- cies between wave amplitudes in experiment and simulation, and is based on a simple application of the conservation of wave action. In this application we have assumed that no other processes are actively influencing the wave amplitude, however, a dependence of wave amplitude on R and Pr has been noted. Chapter 2. Holmboe wave fields in simulation and experiment 31 Our comparisons of simulation and experiments reveal ways in which the nonlinear evolution of a Holmboe wave field is dependent on the mean flow, and hence, on the boundary conditions. This must be considered when studying the nonlinear behaviour of the Holmboe instability. 32 Bibliography Alexakis, A. (2005). On Holmboe’s instability for smooth shear and density profiles. Phys. Fluids, 17:084103. Balmforth, N. and Mandre, S. (2004). Dynamics of roll waves. J. Fluid Mech., 514:1–33. Browand, F. and Winant, C. (1973). Laboratory observations of shear-layer instability in a strat- ified fluid. Boundary-Layer Met., 5:67–77. Carpenter, J., Lawrence, G., and Smyth, W. (2007). Evolution and mixing of asymmetric Holm- boe instabilities. J. Fluid Mech., 582:103–132. Caulfield, C. and Peltier, W. (2000). The anatomy of the mixing transition in homogenous and stratified free shear layers. J. Fluid Mech., 413:1–47. Haigh, S. (1995). Non-symmetric Holmboe Waves. PhD thesis, University of British Columbia. Hazel, P. (1972). Numerical studies of the stability of inviscid stratified shear flows. J. Fluid Mech., 51:39–61. Hogg, A. M. and Ivey, G. (2003). The Kelvin-Helmholtz to Holmboe instability transition in stratified exchange flows. J. Fluid Mech., 477:339–362. Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geofys. Publ., 24:67–112. Huang, N., Long, S., and Shen, Z. (1996). The mechanism for frequency downshift in nonlinear water wave evolution. Adv. Appl. Mech., 32:59–117. Huang, N., Shen, Z., and Long, S. (1999). A new view of nonlinear water waves: the Hilbert spectrum. Ann. Rev. Fluid Mech., 31:417–457. Koppel, D. (1964). On the stability of a thermally stratified fluid under the action of gravity. J. Meth. Phys., 5:963–982. Lawrence, G., Browand, F., and Redekopp, L. (1991). The stability of a sheared density interface. Phys. Fluids, 3(10):2360–2370. Pawlak, G. and Armi, L. (1996). Stability and mixing of a two-layer exchange flow. Dyn. Atmos. Oceans, 24:139–151. Peregrine, D. (1976). Interaction of water waves and currents. Adv. Appl. Mech., 16:9–117. Bibliography 33 Pouliquen, O., Chomaz, J., and Huerre, P. (1994). Propagating Holmboe waves at the interface between two immiscible fluids. J. Fluid Mech., 266:277–302. Sargent, F. and Jirka, G. (1987). Experiments on saline wedge. J. Hydraul. Engng ASCE, 113(10):1307–1324. Sharples, J., Coates, M., and Sherwood, J. (2003). Quantifying turbulent mixing and oxygen fluxes in a Mediterranean-type, microtidal estuary. Ocean Dyn., 53:126–136. Smyth, W. (2006). Secondary circulations in Holmboe waves. Phys. Fluids, 18:064104. Smyth, W., Carpenter, J., and Lawrence, G. (2007). Mixing in symmetric Holmboe waves. J. Phys. Oceanogr., 37:1566–1583. Smyth, W., Klaassen, G., and Peltier, W. (1988). Finite amplitude Holmboe waves. Geophys. Astrophys. Fluid Dyn., 43:181–222. Smyth, W., Moum, J., and Caldwell, D. (2001). The efficiency of mixing in turbulent patches: In- ferences from direct simulations and microstructure observations. J. Phys. Oceanogr., 31:1969– 1992. Smyth, W., Nash, J., and Moum, J. (2005). Differential diffusion in breaking Kelvin-Helmholtz billows. J. Phys. Oceanogr., 35:1004–1022. Smyth, W. and Peltier, W. (1990). Three-dimensional primary instabilities of a stratified, dissi- pative, parallel flow. Geophys. Astrophys. Fluid Dyn., 52:249–261. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Tedford, E., Pieters, R., and Lawrence, G. (2009). Holmboe instabilities in a laboratory exchange flow. J. Fluid Mech., page submitted. Tulin, M. and Waseda, T. (1999). Laboratory observations of wave group evolution, including breaking effects. J. Fluid Mech., 92:197–232. Wesson, M. and Gregg, M. (1994). Mixing at the Camarinal Sill in the strait of Gibraltar. J. Geophys. Res., 99:9847–9878. Winters, K., MacKinnon, J., and Mills, B. (2004). A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Ocean. Tech., 21:69–94. Yonemitsu, N., Swaters, G., Rajaratnam, N., and Lawrence, G. (1996). Shear instabilities in arrested salt-wedge flows. Dyn. Atmos. Oceans, 24:173–182. Bibliography 34 Yoshida, S., Ohtani, M., Nishida, S., and Linden, P. (1998). Mixing processes in a highly strat- ified river. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 389–400. American Geophysical Union. Zhu, D. and Lawrence, G. (2001). Holmboe’s instability in exchange flows. J. Fluid Mech., 429:391–409. 35 Chapter 3 Note on the ejection process in Holmboe waves 1 3.1 Introduction The Holmboe instability is one of a number of possible instability types that may evolve from the stratified shear layer (Holmboe, 1962; Howard and Maslowe, 1973). It results from the in- teraction between disturbances in the shear layer vorticity and the density stratification (Baines and Mitsudera, 1994; Caulfield, 1994), and requires a relatively thin density interface compared to the shear layer thickness. It is able to grow in strong stratifications (i.e. large bulk Richardson numbers), making it a potentially important process in geophysical flows. At finite amplitude, the instability develops into a series of propagating internal waves that often appear cusp-like in form. Depending on conditions, the Holmboe waves may form on one, or both sides of the density interface. Since the earliest laboratory studies of Holmboe waves, investigators have observed a pro- cess whereby interfacial fluid is drawn away from the density interface and entrained into the homogeneous layers on either side (see e.g. Thorpe, 1968; Browand and Winant, 1973; Koop and Browand, 1979). The resulting structures have been referred to by many different names such as wisps, plumes, puffs and ejections, and it is likely that they result from different physical pro- cesses. However, we shall focus on what we observe to be the most conspicuous of these events, in which interfacial fluid is ejected from the wave crests and often travels a considerable dis- tance from the density interface. A snap-shot of this ‘ejection process’ is shown in the laboratory photograph of figure 3.1. Although the ejection process is generally considered to be important for the generation of three-dimensional motions (Smyth and Peltier, 1991), mixing of the density field (Koop and Browand, 1979; Smyth and Winters, 2003), and the basic evolution of the Holm- boe wave field (chapter 2), relatively little is known of the mechanism by which it occurs. It is therefore, the goal of this chapter to provide a mechanistic description of the ejection process. In doing so, we will also address a number of questions raised in the literature regarding the ejections (see Smyth and Peltier, 1991; Hogg and Ivey, 2003, in particular). These include: (i) whether or not ejections are dependent on three-dimensional effects; (ii) if the interaction between the up- per and lower waves may trigger an ejection; and (iii) whether the ejections can be a source of three-dimensional motion. The present work also compliments that of Smyth (2006) on secondary circulations in Holmboe waves. 1A version of this chapter is in preparation for publication. J.R. Carpenter and G.A. Lawrence (2009), Note on the ejection process in Holmboe waves. Chapter 3. Ejection process in Holmboe waves 36 Figure 3.1: Laboratory photograph of a Holmboe wave ejection. The colours are representative of the fluid density and are measured by the light intensity of fluorescing dye. See Tedford et al. (2009) for details of the experimental apparatus. 3.2 Methods Our examination is based on direct numerical simulations (DNS) of the stratified shear layer. The code has been described previously by Winters et al. (2004), but has been modified as in Smyth et al. (2005). The DNS solves the three-dimensional incompressible Boussinesq equa- tions of motion in a rectangular domain that has periodic boundary conditions in the streamwise (x), and spanwise (y) directions, and a free-slip condition on the vertical (z) boundaries. The flow is initialized with random noise in the velocity field in the manner described by Smyth and Winters (2003). This is superimposed on ‘tanh’ background profiles of horizontal velocity U(z) = (∆U/2) tanh[2(z − d)/h], and density ρ̄(z) = ρ0 − (∆ρ/2) tanh(2Rz/h), where ∆U is the velocity difference across the shear layer, h is the shear thickness, ∆ρ is the density dif- ference across the interface, ρ0 is a reference density, and R is the ratio of the shear thickness to the thickness of the density interface. For reasons discussed below, we allow the density inter- face to be vertically offset from the centre of the shear layer by a distance d, and measure this degree of asymmetry by a = 2h/d. Denoting the kinematic viscosity of the fluid by ν, and the diffusivity of the stratifying agent by κ, we identify three more important dimensionless parame- ters J = ∆ρgh/ρ0(∆U)2, Re = ∆Uh/ν, and Pr = ν/κ, as bulk Richardson, Reynold’s, and Prandtl numbers, respectively. We have chosen conditions in the laboratory experiments of Tedford et al. (2009) as a guide, giving Re = 630 and J = 0.30. Since the large values of R = 8 and Pr ≈ 700 in Tedford et al. (2009) place a high demand on computational resources, we have taken R = 5 and Pr = 25; sufficiently large to ensure that Holmboe’s instability develops, yet low enough to adequately re- solve the density interface within our limited computational resources. Two different simulations Chapter 3. Ejection process in Holmboe waves 37 x (a) (b) (c) leading vortexbaroclinic negative vorticity positive shear layer vorticity density interface Figure 3.2: Illustration of the initial formation of asymmetric Holmboe waves. Time is increasing in each panel from (a-c), and each is in a frame of reference moving with the wave speed. are performed, one with a = 0 that results in ‘symmetric’ Holmboe waves both above and below the density interface, and the other an ‘asymmetric’ case with a = 0.5 where waves only form in the upper layer. Ejections occurred most often when multiple waves were present in the domain. We therefore chose a domain length of Lx = 64h for the symmetric simulation, and Lx = 16h in the asymmetric case, which allow for the initial growth of 16 waves and 4 waves, respectively. To accommodate the large number of waves in the a = 0 simulation we choose a width Ly = 2.5h, whereas Ly = 5h in the asymmetric case. In both simulations the depth is Lz = 5.4h as in Tedford et al. (2009). 3.3 Ejection mechanism From the initial random noise of the velocity perturbation, the Holmboe instability acts together with the periodic streamwise boundary conditions to selectively amplify the most unstable wave- length from linear theory. The instability growth consists of a ‘rolling-up’ of the shear layer vorticity on one or both sides of the density interface, depending on a. This process is illustrated in figure 3.2 for a = 0.5, where the waves are formed only in the upper layer. The positive (clockwise) shear layer vorticity indicated by the light grey shading can be seen to concentrate into spanwise oriented leading vorticies that travel with the waves on the density interface (in this case from left to right). The baroclinic generation of vorticity that occurs in regions of horizontal density gradients results in the accumulation of negative (counter-clockwise) vorticity in the wave crests, and positive vorticity in the troughs, similar to a linear internal wave. However, the negative vorticity becomes concentrated in a narrow region about the wave crests, as seen in figure 3.2(c). The time sequence of density and spanwise vorticity fields shown in figure 3.3 demonstrates the ejection process in asymmetric Holmboe waves. A thin dense wisp of fluid is drawn from the Chapter 3. Ejection process in Holmboe waves 38 10-1-2 z/ h 2 7 2 5 2 7 z/ h 5 10 2 5 5 10 z/ h 7 12 2 5 7 12 z/ h x/h10 15 2 5 x/h10 15 (a) (e) (b) (f) (c) (g) (d) (h) Figure 3.3: Sequence of density (a-d) and vorticity (e-h) fields taken at times of t∆U/h = {88, 97, 105, 116} demonstrating the ejection process. The vorticity scale, made dimensionless by h/∆U , is given below with the thin dark contours in (e-h) representing isopycnals spaced every ∆ρ/8. Chapter 3. Ejection process in Holmboe waves 39 x/h t Δ U /h (a) 0 20 40 60 80 100 120 140 160 180 200 x/h (b) 0 20 40 60 Figure 3.4: Characteristics of the right- and left-propagating waves for the symmetric simulation in (a) and (b), respectively. Colour shading represents the elevation of the density interface, but has been reversed in (a) and (b) so that red indicates a wave crest (either upwards deflections in (a), and downwards deflections in (b)) and blue a wave trough. The approximate locations and times of ejections are marked by circles, which occur within the region indicated by the dashed rectangle. wave crest and transported into the homogeneous upper layer where it becomes advected with the background flow. Inspection of the vorticity field in figure 3.3 shows that the vertical transport of ejected fluid against buoyancy forces is accomplished by the formation of a vortex couple, i.e. two vorticies of opposite sign. The vortex couple is composed of the negative baroclinically generated vorticity of the wave crest and the positive vorticity in the leading vortex that originated from the rolling-up of the shear layer. The two oppositely signed vorticies induce an upwards velocity in each other that results in an upwards translation of the couple. Simple models of buoyant vortex couples have been proposed by Turner (1960) and Saffman (1972), as well as others. In these models, a dense vortex couple is able to rise against the downward buoyant force by consuming the linear impulse of the vorticity distribution. Although the application of these models in the present context is complicated by many factors such as complex vorticity and density fields, as well as the presence of a mean shear flow, the basic mechanism is the same. Our observation of asymmetric ejections, that occur in the absence of waves in the lower layer, indicates that the interaction between upper and lower wave modes is not a necessary ingredient in the process. However, it is found that the presence of a second wave mode is effective in triggering ejections. This is demonstrated in figure 3.4, with an x–t characteristic diagram of the density interface elevation in the a = 0 simulation. A two-dimensional Fourier filter has been applied to separate the wave characteristics into only the right-propagating waves in figure 3.4(a) Chapter 3. Ejection process in Holmboe waves 40 and the left-propagating waves in 3.4(b). Four ejections are observed in both the upper and lower waves at similar times and locations. They are triggered as two large-amplitude wave packets in the upper and lower waves meet in the region identified by the dashed box. Interaction of the upper and lower wave modes is known to produce oscillations in phase speed and wave amplitude (Smyth et al., 1988; Zhu and Lawrence, 2001; Hogg and Ivey, 2003), which can also be seen in the ‘wiggling’ and ‘pulsing’ of the characteristics in figure 3.4. These periodic changes in amplitude will also periodically increase the strength of the negative vorticity in the wave crests, thus leading to a stronger vortex couple and a greater chance of ejection. A feature of the ejections, as seen in the wave characteristics, is the sudden reduction in wave amplitude that results. This can also be seen in the reduced amplitude of the wave in the final stage of the ejection process shown in figure 3.3(d,h). Since the ejections are found to occur in only the largest amplitude (steepest) waves, it was proposed in chapter 2 that the process be considered a type of Holmboe wave breaking. This is in general agreement with the mechanism discussed above since a steep wave crest is able to generate stronger concentrations of negative baroclinic vorticity required in the formation of the vortex couple. Previous investigations have noted the ejection process to occur erratically, raising specula- tion that the cause may be linked to three-dimensional effects (Hogg and Ivey, 2003). Since the Holmboe instability is found to be a two-dimensional instability of the basic stratified shear layer over most of parameter space (Smyth and Peltier, 1990; Haigh, 1995), and the mechanism just described is also two-dimensional, there is no need to appeal to the third dimension in order to explain the occurrence of ejections. Indeed, the ejections have been observed in strictly two- dimensional simulations (not shown). However, three-dimensional motions are observed to form on the ejected fluid even when the wave crests are highly two-dimensional, suggesting that they may be a source of three-dimensionality. This three-dimensionality has also been observed in laboratory experiments (Thorpe, 1968; Maxworthy and Browand, 1975). It is possible that the er- ratic occurrence of the ejections is related to both interactions between the upper and lower modes as shown above, as well as from energy transfers occurring from within the wave field (e.g. a modulational instability). 3.4 Conclusions We have a described a mechanism to explain the ejection process that is commonly observed in Holmboe waves. The basic mechanism is two-dimensional, relying on the formation of a vortex couple to transport interfacial fluid against buoyancy. Although it is not necessary for two Holm- boe wave modes to be present, ejections can be triggered by an interaction between the upper and lower wave modes. The ejections also appear to be a source of three-dimensional motions, in agreement with observations made in laboratory experiments. 41 Bibliography Baines, P. and Mitsudera, H. (1994). On the mechanism of shear flow instabilities. J. Fluid Mech., 276:327–342. Browand, F. and Winant, C. (1973). Laboratory observations of shear-layer instability in a strat- ified fluid. Boundary-Layer Met., 5:67–77. Caulfield, C. (1994). Multiple linear instability of layered stratified shear flow. J. Fluid Mech., 258:255–285. Haigh, S. (1995). Non-symmetric Holmboe Waves. PhD thesis, University of British Columbia. Hogg, A. M. and Ivey, G. (2003). The Kelvin-Helmholtz to Holmboe instability transition in stratified exchange flows. J. Fluid Mech., 477:339–362. Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geofys. Publ., 24:67–112. Howard, L. and Maslowe, S. (1973). Stability of stratified shear flows. Boundary-Layer Met., 4:511–523. Koop, C. and Browand, F. (1979). Instability and turbulence in a stratified fluid with shear. J. Fluid Mech., 93:135–159. Maxworthy, T. and Browand, F. (1975). Experiments in rotating and stratified flows: oceano- graphic application. Ann. Rev. Fluid Mech., 7:273–305. Saffman, P. G. (1972). The motion of a vortex pair in a stratified atomosphere. Stud. Appl. Maths, 51(2):107–119. Smyth, W. (2006). Secondary circulations in Holmboe waves. Phys. Fluids, 18:064104. Smyth, W., Klaassen, G., and Peltier, W. (1988). Finite amplitude Holmboe waves. Geophys. Astrophys. Fluid Dyn., 43:181–222. Smyth, W., Nash, J., and Moum, J. (2005). Differential diffusion in breaking Kelvin-Helmholtz billows. J. Phys. Oceanogr., 35:1004–1022. Smyth, W. and Peltier, W. (1990). Three-dimensional primary instabilities of a stratified, dissi- pative, parallel flow. Geophys. Astrophys. Fluid Dyn., 52:249–261. Bibliography 42 Smyth, W. and Peltier, W. (1991). Instability and transition in finite-amplitude Kelvin-Helmholtz and Holmboe waves. J. Fluid Mech., 228:387–415. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Tedford, E., Pieters, R., and Lawrence, G. (2009). Holmboe instabilities in a laboratory exchange flow. J. Fluid Mech., page submitted. Thorpe, S. (1968). A method of producing a shear flow in a stratified fluid. J. Fluid Mech., 32:693–704. Turner, J. S. (1960). A comparison between buoyant vortex rings and vortex pairs. J. Fluid Mech., 7:419–432. Winters, K., MacKinnon, J., and Mills, B. (2004). A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Ocean. Tech., 21:69–94. Zhu, D. and Lawrence, G. (2001). Holmboe’s instability in exchange flows. J. Fluid Mech., 429:391–409. 43 Chapter 4 Observation and analysis of shear instability in the Fraser River estuary 1 4.1 Introduction Shear instabilities occur in highly stratified estuaries and can influence the large scale dynamics by redistributing mass and momentum. Specifically, shear instabilities have been found to influence salinity intrusion in the Fraser River estuary (Geyer and Smith, 1987; Geyer and Farmer, 1989; MacDonald and Horner-Devine, 2008). We describe recent observations in this estuary and exam- ine the shear and stratification that lead to instability. The influence of long time scale processes such as freshwater discharge and the tidal cycle are also discussed. Rather than relying on a bulk or gradient Richardson number to assess stability we use numer- ical solutions of the Taylor-Goldstein (TG) equation based on observed profiles of velocity and density. This approach has been used with some success in the ocean (e.g. Moum et al., 2003) but, with the exception of the simplified application by Yoshida et al. (1998), has not been applied in estuaries. Solving the TG equation provides the growth rate, wavelength, phase speed and mode shape of the instabilities. We compare these predicted wave properties with instabilities observed using an echosounder. Geyer and Farmer (1989) found that instabilities in the Fraser River estuary were most ap- parent during ebb tide when strong shear occurred over the length of the salinity intrusion. They outlined a progression of three phases of increasingly unstable flow that occurs over the course of the ebb. In the first phase, strain sharpens the density interface; shear is stronger than during flood, but insufficient to cause shear instability. In the second phase, the lower layer reverses direction causing shear between the fresh and saline layers to increase. Shear instability and tur- bulent mixing are concentrated at the pycnocline rather than in the bottom boundary layer. By the third phase of the ebb, shear instability has completely mixed the two layers leaving homoge- neous water throughout the depth. During flood there is some mixing, however it is concentrated at the front located at the landward tip of the salinity intrusion. Similarly, MacDonald and Horner- Devine (2008), studying mixing at high fresh water discharge (7000 m3s−1), found that two to three times more mixing occurred during ebb tide than during flood. The present analysis is fo- cused on the ebb tide at high and low freshwater discharge, although some results during flood tide are also presented. 1A version of this chapter has been accepted for publication. E.W. Tedford, J.R. Carpenter, R. Pawlowicz, R. Pieters and G.A. Lawrence (2009) Observation and analysis of shear instability in the Fraser River estuary, J. Geophys. Res. Chapter 4. Shear instability in the Fraser River estuary 44 Sand Heads (0 km) Massey Tunnel (18 km) T1 T2T3 T4 T5T6 N6 km0 2 4 Figure 4.1: Map of the lower 27 km of the Fraser River. The locations of the six transects are marked T1-T6. The mouth of the river (Sand Heads) is located at 49◦ 6′ N and 123◦ 18′ W. Discharge Tide x ∆U ∆ρ h J (m3 s−1) (km) (m s−1) (kg m−3) (m) 1 6400 Ebb 8.6 1.6 14.3 5.2 0.29 2 6500 Ebb 11 1.65 20 3.5 0.25 3 5700 Flood 2.2 1.5 23.1 3.5 0.35 4 850 Ebb 24.5 1.5 12.9 12 1.3 5 850 Ebb 19 1.5 12.9 12 1.3 6 850 Ebb 10.5 2.5 7.3 12 0.3 Table 4.1: Details of transects shown in figures 4.1 and 4.2. The location indicates the distance upstream from the mouth (Sand Heads). This chapter is organized as follows. The setting and field methods are described in section 4.2. The general structure of the salinity intrusion is described in section 4.3. In section 4.4 we present the background theory needed to perform stability analysis in the Fraser River estuary. In section 4.5 predictions from the stability analysis are compared with observations. In section 4.6 the source of relatively small scale overturning is briefly discussed. In section 4.7 the results of the stability analysis are discussed followed by conclusions in section 4.8. 4.2 Site Description and Data Collection Data were collected in the main arm of the Fraser River estuary, British Columbia, Canada (fig- ure 4.1). The estuary is 10 to 20 m deep with a channel width of 600 to 900 m. Cruises were conducted on June 12, 14 and 21, 2006 and March 10, 2008. Here we present one transect from each of the June 2006 cruises and three transects from the March 2008 cruise (see Table 4.1). The freshwater discharge during the June 2006 transects was typical of the freshet at approxi- mately 6000 m3s−1. During the March 2008 transects, freshwater discharge was near the annual minimum at 850 m3s−1. In June 2006, transects were made during both ebb and flood tide. In March 2008, transects cover most of a single ebb tide (figure 4.2). The tides in the Strait of Geor- Chapter 4. Shear instability in the Fraser River estuary 45 0 2 4 6 8 10 12 -2 0 2 H ei g h t (m ) T1 June 12, 2006. Q fresh = 6400 m3s-1 2 4 6 8 10 12 14 -2 0 2 H ei g h t (m ) T2 8 10 12 14 16 18 -2 0 2 H ei g h t (m ) T3 6 8 10 12 14 16 18 -2 0 2 H ei g h t (m ) Local Time (hours, PDT) T4 T5 T6 June 14, 2006. Q fresh = 6500 m3s-1 June 21, 2006. Q fresh = 5700 m3s-1 March 10, 2008. Q fresh = 850 m3s-1 Figure 4.2: Observed tides at Point Atkinson (heavy line) and New Westminster (thin line) for the four days of field observations. The Point Atkinson data is representative of the tides in the Strait of Georgia beyond the influence of the Fraser River. New Westminster is located 37 km upstream of the mouth of the river at Sand Heads (see figure 4.1). The records are both referenced to mean sea level at Point Atkinson. The duration of the six transects are marked T1-T6. gia have M2 and K1 components of similar amplitude (approximately 1 m) resulting in strong diurnal variations. The tidal range varies from approximately 2 m during neap tides to approxi- mately 4.5 m during spring tides. During both the 2006 and 2008 observations the tidal range was approximately 3 m. The distance salinity intrudes landward of Sand Heads, i.e. the total length of the salinity intrusion, varies considerably with tidal conditions and freshwater discharge. Ward (1976), found the maximum length of the intrusion occurred just after high tide and varied from 8 km at high discharge (9000 m3s−1) to 31 km at low discharge (850 m3s−1). Geyer and Farmer (1989) found that, at average discharge (3000 m3s−1), the maximum length of the intrusion matched the hori- zontal excursion of the tides (10 to 20 km) and, similar to Ward (1976), occurred just after high tide. Kostachuk and Atwood (1990) found that the minimum length of the salinity intrusion typi- Chapter 4. Shear instability in the Fraser River estuary 46 cally occurred approximately one hour after low tide. The longest intrusion they observed at low tide was approximately 20 km. They predicted that complete flushing of salt from the estuary would occur on most days during the freshet (freshwater discharge > 5000 m3s−1). Field methods Data along the six transects were collected by drifting seaward with the surface flow while logging velocity and echosounder data and yoyoing a CTD (conductivity, temperature and depth) profiler. The velocity measurements were made with a 1200 kHz RDI Acoustic Doppler Current Profiler (ADCP) sampling at 0.4 Hz with a vertical resolution of 250 mm. The velocities were averaged over 60 seconds to remove high frequency variability. The echo soundings were made with a 200 kHz Biosonics sounder sampling at 5 Hz with a vertical resolution of 18 mm. Profile data was collected with a Seabird 19 sampling at 2 Hz. Selected echosounder, ADCP and CTD data are shown in figure 4.3. As indicated by the superimposed density profiles, strong gradients in density are generally associated with a strong echo from the sounder. The CTD was profiled on a load bearing data cable that provided constant monitoring of conductivity, temperature and depth. These data allowed us to quickly identify the front of the salinity intrusion and avoid direct contact of the instrument with the bottom. To increase the vertical resolution of the profiles, the CTD was mounted horizontally with a fin to direct the sensors into the flow. In this configuration, the instrument was allowed to descend rapidly and then was raised slowly (0.2 - 0.4 m s−1) relying on horizontal velocity of the water relative to the CTD to flush the sensors. The upcast, which had higher vertical resolution, was in reasonable agreement with the echo intensity from the sounder. On the few occasions that the higher resolution upcast did not coincide with the appearance of instabilities in the echosounder, we used the downcast. The total number of CTD casts we were able to perform varied from transect to transect depending on field conditions (surface velocity, shear, ship traffic, woody debris). 4.3 General Description of the Salinity Intrusion We observed important differences in the structure of the salinity intrusion between high and low freshwater discharge. At high discharge, our observations were similar to those described by Geyer and Farmer (1989), where the salinity intrusion had a two-layer structure resembling a classic salt-wedge. At low discharge, however, the salinity intrusion exhibited greater complexity. 4.3.1 High Discharge During flood tide, mixing was concentrated near the steep front at the landward tip of the salt- wedge (2.7 to 3.03 km in figure 4.3c). During ebb tide, the steep front was replaced by a gently Chapter 4. Shear instability in the Fraser River estuary 47 (a) Transect 1 8.6 8.7 8.8 8.9 9 9.1 9.2 9.3 9.4 9.5 9.6 2 4 6 8 10 12 D ep th ( m ) (b) Transect 2 1 m s-1 25 kg m 10.8 10.9 11 11.1 11.2 11.3 11.4 11.5 11.6 11.7 11.8 6 8 10 12 14 16 Distance (km) (c) Transect 3 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 2 4 6 8 10 12 -3 Figure 4.3: Echo soundings observed during high discharge on: (a) transect 1, ebb tide; (b) tran- sect 2, ebb tide; (c) transect 3, flood tide. The shading scales with the log of the echo intensity with black corresponding to the strongest echos. Selected velocity profiles (red) from the ADCP and density profiles (blue dash) from the CTD are superimposed (not all are shown). The black line indicates the location of the boat in the middle of the cast, as well as the zero reference for the velocity and σT . The velocity profile was calculated as a 1 minute average centred on the time of the CTD cast. The undulations in the bed of the river (thick black line at the bottom of the echosoundings) are a result of sandwaves. Chapter 4. Shear instability in the Fraser River estuary 48 sloping pycnocline (figure 4.3b, landward of 11.6 km) and there was no apparent concentration of mixing at the landward tip of the salt-wedge (not shown). We will focus on the wave-like disturbances that occur on the pycnocline especially during ebb tide. The largest of these were observed during transect 1 (figure 4.3a, 8.7 to 8.9 km between depths of 3 and 10 m). These disturbances occurred within the upper layer as it passed over the nearly stationary water below a depth of 10 m. Smaller amplitude disturbances were observed during transect 2 (figure 4.3b, 11.05 km). In our application of the TG equation we will show that disturbances like these are a result of shear instability. Not all of the disturbances on the pycnocline are a result of shear instability. For example, for most of the velocity and density profiles collected during transect 3 (figure 4.3c) the TG equation does not predict instability. The disturbances seen from 2.5 to 3.0 km are caused by the large sand waves on the bottom (the thick black line in the echo sounding). The crests of the sand waves were typically 30 m apart and 1 to 2 m high, and were found over most of the river surveyed during high discharge (2.5 km to 15 km). During flood tide, flow over these sand waves caused particularly regular disturbances on the pycnocline. 4.3.2 Low Discharge At low discharge, at the beginning of the ebb, the front of the salinity intrusion was located be- tween 28 and 30 km from Sand Heads. Unlike the observations at high discharge a well defined front was not visible in the echosounder, and CTD profiles were needed to identify its location. Seaward of the front (figure 4.4a), the echosounder and the CTD profiles show a multilayered structure with more complexity than was observed at high discharge. At this early stage of the ebb, the CTD profiles generally show partially mixed layers separated by several weak density interfaces. Later in the ebb, during transect 5 (figure 4.4b), near-bottom velocities turn seaward and the velocity shear between the top and the bottom increases. At maximum ebb (transect 6, figure 4.4c), the shear increases further, reaching a maximum of approximately 2.5 m s−1 over a depth of 12 m. Mixed water occurs at both the surface and the bottom resulting in an overall decrease in the vertical density gradient. By the time transect 6 is complete the ebb flow is decelerating. The salinity intrusion continues to propagate seaward until low tide but, given its length and velocity it does not have sufficient time to be completely flushed from the estuary. During the next flood the mixed water remaining in the estuary allows a complex density structure to develop similar to that seen early in the observed ebb. This differs from the behaviour at high freshwater discharge when nearly all of the seawater is flushed completely from the estuary at least once a day. Chapter 4. Shear instability in the Fraser River estuary 49 (a) Transect 4 21.5 22 22.5 23 23.5 24 24.5 25 0 5 10 15 D ep th ( m ) (b) Transect 5 17.5 18 18.5 19 19.5 20 20.5 0 5 10 15 Distance (km) (c) Transect 6 7.5 8 8.5 9 9.5 10 10.5 11 0 5 10 15 2 m s-1 20 kg m -3 Figure 4.4: Echo soundings during low discharge observed during: (a) transect 4, early ebb; (b) transect 5, mid ebb; and (c) transect 6, late ebb. The shading scales with the log of the echo intensity with black corresponding to the strongest echos. Note that the scale of the shading is the same in all three panels. Velocities (red) from the ADCP and densities (blue dashed) from the CTD are superimposed. The black line indicates the location of the boat in the middle of the cast, as well as the zero reference for the velocity and σT . The velocity profile was calculated as a 1 minute average centred on the time of the CTD cast. Chapter 4. Shear instability in the Fraser River estuary 50 4.4 Stability of Stratified Shear Flows 4.4.1 Taylor-Goldstein Equation Following Taylor (1931) and Goldstein (1931) we assess the stability of the flow by consider- ing the evolution of perturbations on the background profiles of density and horizontal velocity, denoted here by ρ̄(z) and U(z), respectively. If the perturbations to the background state are sufficiently small they are well approximated by the linear equations of motion. It then suffices to consider sinusoidal perturbations, represented by the normal mode form eik(x−ct), where x is the horizontal position and t is time. Here k = 2pi/λ is the horizontal wave number with λ the wavelength, c = cr + ici is the complex phase speed. If we further assume that the flow is incompressible, Boussinesq, inviscid, and non-diffusive, we arrive at the Taylor-Goldstein (TG) equation ψ̂′′ + [ N2 (U − c)2 − U ′′ U − c − k 2 ] ψ̂ = 0, (4.1) where primes represent ordinary differentiation with respect to height (z), the stream function is given by ψ(x, z, t) = ψ̂(z)eik(x−ct) and N2(z) = −gρ̄′/ρ0 represents the Boussinesq form of the squared buoyancy frequency with a reference density, ρ0. We will find solutions to the TG equation that consist of eigenfunction-eigenvalue sets {ψ̂(z), c}, for each value of k. Each set {ψ̂(z), c} is referred to as a mode, and the solution may consist of the sum of many such modes for a single k. The background flow, represented by U(z) and ρ̄(z), is then said to be unstable if any modes exist that have ci 6= 0. In this case the small perturbations grow exponentially at a rate given by kci. In general, unstable modes are found over a range of k, and it is the mode with the largest growth rate that is likely to be observed. Although they are based on linear analysis, TG predictions of the wave properties, k and c, typically match those of finite amplitude instabilities observed in the laboratory (Thorpe, 1973; Lawrence et al., 1991; Tedford et al., 2009). 4.4.2 Miles-Howard criterion A useful criterion to assess the stability of a given flow without solving the TG equation was derived by Miles (1961) and Howard (1961). They found that if the gradient Richardson number, Ri(z) = N2/(U ′)2, exceeds 1/4 everywhere in the profile, then the TG equation has no unstable modes, i.e. ci must be zero for all modes. In other words, Ri ≥ 1/4 everywhere is a sufficient condition for stability, referred to as the Miles-Howard criterion. Note that if Ri < 1/4 at some location, instability is possible, but not guaranteed. Despite the inconclusive nature of the Miles-Howard criterion for determining instability, it is often employed as a sufficient condition for instability in density stratified flows, and has been found to have reasonable agreement with observations (Thorpe, 2005, p. 201-204). Looking Chapter 4. Shear instability in the Fraser River estuary 51 specifically at the Fraser River estuary, Geyer and Smith (1987) were able to compute statistics of Ri and show that decreases in Ri were accompanied by mixing in the estuary. 4.4.3 Mixing Layer Solution Since the TG equation is an eigenvalue problem with variable coefficients, analytical solutions can only be obtained for the simplest profiles, and recourse is usually made to numerical methods (e.g. Hazel, 1972). However, the available analytical solutions are often a useful point of departure. We look at one such solution that closely approximates conditions found in the estuary during high discharge. This solution is based on the simple mixing layer model of Holmboe (described in Miles, 1963). In a general form of the model, the velocity and density profiles are represented by the hyper- bolic tangent functions U(z) = ∆U 2 tanh [2(z − d) h ] and ρ̄(z) = −∆ρ 2 tanh (2z δ ) + ρ0. (4.2) where h is the shear layer thickness, δ is the thickness of the density interface, and ∆U , ∆ρ corre- spond to the total change in velocity and density across the layers, respectively. The parameter d allows for a vertical offset in the positions of the shear layer and density interface. In the simplest case the shear layer and density stratification have equal thickness, giving R = h/δ = 1, and they coincide in their vertical positions so that the asymmetry a = 2d/h = 0. In this case, Ri(z) is at its minimum at the centre of the mixing layer (z = 0), and is equal to the bulk Richardson number J = g∆ρh/ρ0(∆U)2. When the bulk Richardson number (i.e. the minimum Ri) drops below 1/4, flows with R = 1 and a = 0 become unstable. The resulting instabilities are of the Kelvin-Helmholtz (KH) type, in which the shear layer rolls up to form an array of billows that are stationary with respect to the mean flow, and which display large overturns in density (Thorpe, 1973). It is not generally the case that J ≥ 1/4 results in stability. For example, if the density interface is relatively sharp (R > 2) an additional mode of instability, the Holmboe mode, is excited (Alexakis, 2005). In this case, the range of J over which instability occurs extends above 1/4. That is, Ri < 1/4 somewhere in z at the same time as J > 1/4. While it is generally true that flows with higher J are subject to less mixing by shear instabilities, by itself, J does not indicate whether or not a flow is unstable. For simplicity, the analytical solution of Holmboe’s mixing layer model assumes the flow is unbounded in the vertical. In our analysis we include boundaries at the top and bottom where ψ̂ must satisfy the boundary condition ψ̂ = 0. The presence of these boundaries tends to extend the range of unstable wavenumber to longer wavelengths (Hazel, 1972). However, in the cases considered here, at the wavenumber of maximum growth, the boundaries have little or no impact Chapter 4. Shear instability in the Fraser River estuary 52 on k and c. 4.4.4 Solution of the TG equation for observed profiles We use the numerical method described in Moum et al. (2003) to generate solutions to the TG equation based on measured velocity and density profiles. Whenever possible we use velocity and density profiles collected at the upstream edge of apparent instabilities in the echosoundings. The velocity profile, a 60 second average, is an average over one or more instabilities (the instabilities have periods< 60 seconds). This averaging reduces the influence of individual instabilities on the velocity profile, which in the TG equation, is taken to represent the background velocity profile. The velocity profile is then smoothed in the vertical using a low pass filter (removing vertical wavelengths < 2 m). The density profile is smoothed by fitting a linear function, and one or more tanh functions (one for each density interface). By using smooth profiles we are effectively ignoring instability associated with small scale variations in the profiles. Because the point of observation moves in time, i.e. the boat is drifting seaward, predicted wavelengths from the TG equation cannot be compared directly to the wavelength of instabilities as they appear in the echosoundings. The wavelength predicted with the TG solution must be shifted to account for the speed of the instabilities with respect to the speed of the boat: λ∗ = ∣∣∣ vb cr − vb ∣∣∣λ. (4.3) Here vb is the velocity of the boat and cr and λ are the phase speed and wavelength predicted with the TG equation. The predicted apparent wavelength, λ∗, is directly comparable to observations made from the moving boat. Seim and Gregg (1994) used a similar approach for estimating the wavelength of observed features. As well as giving a wavelength, phase speed, and growth rate for each unstable mode, the TG solutions also give an eigenfunction that describes the vertical structure of the growing mode. The vertical displacement eigenfunction η̂(z) = −ψ̂/(U − c) is particularly useful. At the location in z where |η̂| is a maximum we expect to see evidence of instabilities in the echosoundings. 4.5 Results In this section we use J , Ri(z) and solutions of the TG equation to assess the stability of six sets of velocity and density profiles (one from each of the six transects). Each set of profiles was chosen to coincide with evidence of instability in the echosoundings. Ebb during high discharge: Transect 1 The selected velocity and density profiles from transect 1 are shown in figure 4.5. The corre- Chapter 4. Shear instability in the Fraser River estuary 53 -2 -1 0 0 2 4 6 8 10 12 14 D ep th ( m ) U (m s ) (a) 0 10 20 ρ - ρ (kg m )-3 (b) Distance (m) (c) ↓ 0 50 100 150 200 250 0 -1 Figure 4.5: Velocity (a) and density (b) profiles observed during transect 1 (June 12, 2006, 8h05 PDT, 8.9 km upstream of Sand Heads). The smooth profiles used in the stability analysis are shown as thick black lines and the observed data are plotted as points. The gray shading indicates regions in which Ri < 1/4. The black horizontal line indicates the location of maximum displacement (|η̂|) for the most unstable mode predicted with the TG equation. The thin lines in (b) show the displacement functions for each of the unstable modes. The functions are scaled in proportion to the growth rate. A close up of the echosounding logged near the location of the profiles is shown in (c), and includes a scale indicating the apparent wavelength predicted by the TG equation. The arrow at the top of image indicates the approximate location of the density and velocity measurements. In this case, the velocity is averaged over a distance of approximately 130 m. sponding value of J for these profiles is 0.29 (see table 4.1). The stability analysis yields two modes of instability. The fastest growing mode is unstable for wavelengths greater than 11 m and has a peak growth rate of 0.025 s−1 (doubling time of 28 s) occurring at a wavelength of 21 m. The phase speed of the instability at this wavelength is -1.02 m s−1, where the negative indicates a seaward direction. Given this phase speed and the seaward drift of the boat (-2.2 m s−1), an apparent wavelength of 39 m is calculated. Chapter 4. Shear instability in the Fraser River estuary 54 Echosoundings collected at the same time, figure 4.5(c), show clear evidence of instabilities. Our interpretation of the echosoundings follows that described by Browning (1971) and Browning et al. (1973). The prediction is found to be similar to, although shorter than, the approximately 50 m wavelength of the observed instabilities. The maximum displacement of the predicted instabil- ities is located at a depth of 7.6 m (indicated by the horizontal line), closely matching the depth of the observed instabilities. Both the observed and predicted instability occur within the region of shear above the maximum gradient in ρ̄ (at a depth of 9 m). As indicated by the gray shading, this region of high shear and low gradient in ρ̄ corresponds to Ri < 1/4. For the set of profiles shown in figure 4.5 the TG equation predicts a second, weaker, unstable mode located at a depth of 2.5 m. This mode is associated with the inflection point (U ′′ = 0) in the velocity profile at this depth. Because there is very little density stratification and hence weak echo intensity at this depth we are unable to confirm or deny the presence of this mode in the echosoundings. Ebb during high discharge: Transect 2 In transect 2 a single hyperbolic tangent gives a good fit to the measured density profile (figure 4.6b). Due to difficulties in profiling, the density profile at this location was missing data below 12 m. Data from the previous cast, taken 60 m upstream, was used below 12 m. This cast is expected to be sampling water of similar density below this depth. In this case the stability analysis of the profiles results in a single mode of instability. The mode is unstable for wavelengths from 10 m to 35 m with a peak growth rate of 0.02 s−1 (doubling time of 35 s) occurring at a wavelength of 17 m. The phase speed of the instability at this wavelength is -0.51 m s−1. Given the drift velocity of -1.9 m s−1, an apparent wavelength of 24 m is calculated. This prediction is found to be similar to, although longer than, the approximately 18 m wavelength of the small instabilities appearing in the echosounding (figure 4.6c). The maximum displacement of the predicted instabilities is located at a depth of 10.6 m, closely matching the depth of the observed instabilities. Flood during high discharge: Transect 3 Despite the occurrence of Ri < 1/4 the stability analysis of the profiles in figure 4.7(a) and 4.7(b) does not find any unstable modes. Echosoundings collected during the flood generally show features on the pycnocline that were well correlated with sand waves (figure 4.7c). These correlated features are likely controlled by the hydraulics of the flow over the sand waves. There was very little evidence of instabilities independent of these sand waves. There appear to be some wave-like features on the pycnocline that are shorter (≈ 10 m) than the sandwaves, however, these are not well resolved by the echosounder (e.g. depth of 9 m at x = 60 m). Properly Chapter 4. Shear instability in the Fraser River estuary 55 -1.5 -1 -0.5 0 0 2 4 6 8 10 12 14 16 0 10 20 ↓ 0 20 40 60 80 100 D ep th ( m ) U (m s ) (a) ρ - ρ (kg m )-3 (b) Distance (m) (c) 0 -1 Figure 4.6: Velocity (a) and density (b) profiles observed during transect 2 (June 14, 2006, 8h21 PDT, 11.1 km upstream of Sand Heads). See figure 4.5 for details. In this case, the velocity is averaged over approximately 110 m. assessing the flow over these sandwaves would require at least two or three sets of density and velocity profiles per sandwave, more than we were able to obtain. Low freshwater discharge Early ebb during low discharge: Transect 4 At low discharge, during the ebb tide, shear and density stratification are spread over the entire depth (see figure 4.4). The bulk shear layer thickness, h, is therefore greater than at high discharge, where shear and stratification were concentrated at a single, relatively thin interface. The increase in the vertical extent of the shear results in a greater bulk Richardson number despite a decrease in the overall strength of the density stratification, ∆ρ (see Table 4.1). The density and velocity profiles collected early in the ebb (transect 4, figure 4.8) consist of a number of layers. The stability analysis yields two modes of instability. The most unstable mode Chapter 4. Shear instability in the Fraser River estuary 56 -0.5 0 0.5 0 2 4 6 8 10 12 14 0 10 20 ↓ 0 50 100 150 200 D ep th ( m ) U (m s ) (a) ρ - ρ (kg m )-3 (b) Distance (m) (c) 0 -1 Figure 4.7: Velocity (a) and density (b) profiles observed during transect 3 (June 21, 2006, 12h38 PDT, 2.66 km upstream of Sand Heads). See figure 4.5 for details. In this case, the velocity is averaged over approximately 30 m. has a peak growth rate of 0.023 s−1 occurring at a wavelength of 10.3 m with a phase speed of -0.86 m s−1. Given this phase speed and the seaward drift of the boat (1.6 m s−1), an apparent wavelength of 22 m is calculated. This is very similar to the wavelength of the largest instability in figure 4.8(c). This mode has a maximum displacement at a depth of 2.5 m, closely matching the location of the observed instabilities. Mid ebb during low discharge: Transect 5 The instabilities in figure 4.9(c) were observed one hour later and approximately 3 km downstream from Transect 4. The ρ̄ profile (figure 4.9b) again displays a number of layers consisting of high- gradient steps. However, the layers are not evident in the measured velocity profile (figure 4.9a), as was the case in figure 4.8, and the overall shape of the velocity profile is more linear. The CTD cast is one of the few collected during the study where the instrument passed through Chapter 4. Shear instability in the Fraser River estuary 57 -1 -0.5 0 0 2 4 6 8 10 12 14 0 5 10 0 10 20 30 40 50 60 70 D ep th ( m ) U (m s ) (a) ρ - ρ (kg m )-3 (b) Distance (m) (c) 0 -1 Figure 4.8: Velocity (a) and density (b) profiles observed during transect 4 (March 10, 2008, 11h20 PDT, 22.4 km upstream of Sand Heads). See figure 4.5 for details. In this case, the velocity is averaged over approximately 90 m. an overturn in the pycnocline (depth of approximately 3.8 m). Consistent with the small ampli- tude of the instabilities in the echosounder, the overturn in the density profile has only water of intermediate density, i.e. no surface or bottom water is observed in the overturn. The TG equation predicts an unstable mode with a peak growth rate (0.03 s−1) at a wave- length of 14 m with a phase speed of -1.2 m s−1. The apparent wavelength is predicted to be 32 m, whereas the features in the echosounder range in horizontal length from approximately 10 to 40 m, with the largest being near the TG prediction (≈ 30 m). The predicted maximum in the displacement eigenfunction occurs at a depth of 4.2 m closely matching the depth of the instabilities. Chapter 4. Shear instability in the Fraser River estuary 58 -1.5 -1 -0.5 0 2 4 6 8 10 12 0 5 10 15 ↓ 0 50 100 150 200 D ep th ( m ) U (m s ) (a) ρ - ρ (kg m )-3 (b) Distance (m) (c) 0 -1 Figure 4.9: Velocity (a) and density (b) profiles observed during transect 5 (March 10, 2008, 12h21 PDT, 19.6 km upstream of Sand Heads). See figure 4.5 for details. In this case, the velocity is averaged over approximately 140 m. Late ebb during low discharge: Transect 6 In the later stages of the ebb, during transect 6 (figure 4.10), the shear has increased such that J is reduced to approximately 0.3. Unlike most of the other profiles collected during low or high discharge the density profile has no homogeneous layers, and shows small scale (i.e. on the scale of the instrument resolution) overturning throughout the depth. In these profiles Ri is below critical throughout most of the depth aside from at the density interface. The most unstable mode predicted with the TG equation is located at a depth of 5.6 m and has a maximum growth rate of 0.019 s−1 at an apparent wavelength of 65m. This is close to, but longer than, the largest features in the echosounder (approximately 50 m). Chapter 4. Shear instability in the Fraser River estuary 59 -2 -1 0 0 2 4 6 8 10 12 5 10 15 ↓ 0 50 100 150 200 250 D ep th ( m ) U (m s ) (a) ρ - ρ (kg m )-3 (b) Distance (m) (c) 0 -1 Figure 4.10: Velocity (a) and density (b) profiles observed during transect 6 (March 10, 2008, 14h34 PDT, 7.6 km upstream of Sand Heads). See figure 4.5 for details. In this case, the velocity is averaged over approximately 130 m. 4.6 Small scale overturns and bottom stress In figure 4.10 there are no features in the echosoundings that are associated with the small scale overturns in ρ̄ below a depth of 7 m, and although our solutions to the TG equation suggest unstable modes, these are both located well above a depth of 7 m. To further examine the source of these overturns we compare selected density profiles from each of the low discharge transects (figure 4.11). In the density profile from transect 4, small scale overturns are rare or completely absent (figure 4.11, T4). Approximately two hours later, during transect 5, just one profile exhibits these small scale overturns (figure 4.11, T5). This cast was performed at the shallow constriction in the river associated with the Massey Tunnel (figure 4.4b, 18 km). In this case the small scale overturns in the profile occur only below the pycnocline suggesting that the stratification within the pycnocline is confining the overturns to the lower layer. By maximum ebb, small scale overturns Chapter 4. Shear instability in the Fraser River estuary 60 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 D ep th ( m ) ρ − ρ (kg m ) T4 T5 T6 −3− 0 Figure 4.11: Selected density profiles from transects performed at low freshwater discharge. The profiles were collected at t=10h53, 12h36 and 14h25, at x=26.2, 17.9 and 8.8 km (transects 4, 5 and 6 respectively). occur throughout the depth (figure 4.11, T6). The presence of these small scale overturns is apparent, although not immediately obvious, in the echosoundings in figure 4.4. Note that the scale of the shading is the same in all three panels of figure 4.4 and that there is a gradual increase (darkening) in background echo intensity from early to late ebb (transects 4 to 6). This increase in echo intensity is attributed to the small scale overturning observed in the density profiles. Early in the ebb the dark shading associated with high echo intensity is concentrated at the density interfaces (transect 4). Otherwise, at this time, echo intensity is low (light shading) corresponding to an absence of small scale overturns in the density profiles (e.g. figure 4.11, T4). At this stage of the ebb, near-bottom velocities are close to zero and bottom stress is expected to be negligible. In transect 5 (figure 4.4b) there is an increase in echo intensity as the flow passes over the Massey Tunnel (18 km). At this location and during this stage of the ebb, near bottom velocity increases to approximately 0.2 m s−1 at 1 m above the bed. In this case the small scale overturns in the profile occur only below the pycnocline (figure 4.11, T5) suggesting that the stratification within the pycnocline is confining bottom generated turbulence to the lower layer. Near maximum ebb, during transect 6, near bottom velocities reach 0.5 m s−1 at 1 m above the bed. By this stage, high echo intensity and small scale overturns occur throughout the depth (figure 4.11, T6) suggesting that bottom generated turbulence has reached the surface Chapter 4. Shear instability in the Fraser River estuary 61 despite the presence of stratification. 4.7 Discussion Although combining echosoundings, velocity, and density measurements to study shear flows is not in itself novel, even for studies in the Fraser estuary (e.g. Geyer and Smith, 1987), efforts in the present study were focussed on simultaneously measuring the details of the flow and the shear instabilities. Our strategy of drifting slowly with the upper-level flow allowed acoustic imaging to capture shear instabilities similar to those observed in laboratory and numerical simulations (e.g. Tedford et al., 2009). Density and velocity measurements also allowed us to analyze these features using a method more typically applied to laboratory experiments, namely direct application of the TG equation. This analysis has refined our understanding of instability and mixing in the Fraser River estuary. One-sided instability In all five of the cases that the TG equation predicted the occurrence of unstable modes, the bulk Richardson number, J , was greater than 1/4. This result suggests the mixing layer model and associated J (see section 4.4.3) are not adequate for describing the stability of the measured profiles. In all of these unstable cases, both the region of Ri(z) < 1/4 and the depth of the maximum in the displacement eigenfunction (|η̂(z)|) were vertically offset from the maximum gradient in density (ρ̄′). This offset between the depth of the predicted region of instability and the density interface is due to asymmetry between the density and velocity profiles. This suggests that a minimum of three bulk parameters (J , R, a) are required if the stability is to be represented by the simplified profiles of equation 4.2. Laboratory models and direct numerical simulations (DNS) of asymmetry result in one-sided instabilities that resemble the features in the echosoundings in figures 4.5(c), 4.6(c) and 4.8(c) (e.g. Lawrence et al., 1991; Yonemitsu, 1991; Carpenter et al., 2007). Similar observations were made in the Strait of Gibraltar by Farmer and Armi (1998) and in a strongly stratified estuary by Yoshida et al. (1998). In these cases the instabilities were attributed to one-sided modes. One- sided modes are part of a general class of instability that includes the Holmboe mode. In contrast to the classic KH mode, the Holmboe mode is a result of the unstable interaction of gradients in density and gradients in shear (N2 and U ′′ in equation 4.1) and can occur at relatively high values of J (Holmboe, 1962). There are a number of potential sources of asymmetry in the Fraser estuary. The most obvious is the difference in the bottom and surface boundary condition. If the stress acting on these bound- aries is not equal and opposite, i.e. if it is unbalanced, then there is the potential for asymmetry. During low freshwater discharge the presence of multiple layers of varying thickness adds further Chapter 4. Shear instability in the Fraser River estuary 62 Kelvin−Helmholtzin- el holtz One-sided Instability Density interface Figure 4.12: Schematic of a Kelvin-Helmholtz and a one-sided instability. The gray shading indicates mixed fluid and the solid black line indicates the position of the central isopycnal (i.e. density interface). irregularity and potential asymmetry to the profiles. Although some laboratory models of strat- ified flows successfully generate symmetric conditions (e.g. Thorpe, 1973; Tedford et al., 2009) many others result in asymmetry (e.g. Lawrence et al., 1991; Yonemitsu et al., 1996; Pawlak and Armi, 1998; Zhu and Lawrence, 2001). In most of these cases asymmetry in the flow results from the geometry of the channel, such as a sill causing localized acceleration of the lower layer. In the arrested salt wedge experiments of Yonemitsu et al. (1996) asymmetry was associated with secondary circulation. Given the common occurrence of asymmetry in the laboratory it is not surprising to find asymmetry in nature. Mixing Linear stability analysis does not provide quantitative predictions of mixing. When one-sided instabilities are modeled using DNS at the values of J observed here the complete overturning of the density interface normally associated with KH billows does not occur. Figure 4.12 shows a schematic of a one-sided instability and a Kelvin-Helholtz instability. Although one-sided in- stabilities are offset from the region of maximum density gradient they have been found to be responsible for considerable mixing (Smyth and Winters, 2003; Carpenter et al., 2007). Unlike the mixed fluid that results from the KH instability, the mixed fluid that results from one-sided in- stabilities is not concentrated at the density interface, but, is instead drawn away from the density Chapter 4. Shear instability in the Fraser River estuary 63 interface (Carpenter et al., 2007). MacDonald and Horner-Devine (2008) quantified mixing in the Fraser estuary at high fresh- water discharge over approximately two tidal cycles. Using a control volume approach and over- turning analysis they estimated mean buoyancy flux, B, during the ebb to be 2.2 ×10−5 m2s−3. The associated mean turbulent eddy diffusivity, K = B/N2, was estimated to be 9×10−4 m2 s−1 (MacDonald, 2003). Smyth et al. (2007) proposed parametrizing the mixing caused by Holmboe instabilities as K = 0.8 × 10−4 h∆U . For transect 1 of the present study, which most closely matches the conditions of MacDonald and Horner-Devine (2008) (see Table 4.1), this results in K = 6.7 × 10−4 m2 s−1 (0.8 × 10−4 × 5.2 m×1.6 m s−1). The parametrization of Smyth et al. (2007) represents the effect of a uniform distribution of instabilities and has not been validated at high Reynolds number. Mindful of these inherent limitations of DNS and the complexity of the field conditions, the similarity between the observed (K = 9×10−4 m2 s−1) and predicted mixing (K = 6.7× 10−4 m2 s−1) is promising. A more rigorous analysis would include a description of the spatial and temporal distribution of instabilities. Unfortunately, our sampling was inadequate to comprehensively describe this distribution particularly at high fresh water discharge. During our survey of the estuary the mixing was apparently caused by shear instabilities acting within the interior of the flow and, to a lesser extent, by turbulence associated with the bottom boundary. Although we have addressed these two types of mixing separately they both originate as a shear instability. Unlike instability predicted with the TG equation the instability associated with the bottom boundary layer relies on viscous effects and the presence of the solid boundary. In some cases, for example during late ebb at low freshwater discharge (transect 6, figure 4.10), the two mechanisms (TG-type instabilities and viscous shear instabilities) are acting together to generate mixing. At high freshwater discharge, during the ebb, MacDonald and Horner-Devine (2008) found that mixing at the pycnocline causes a collapse of the salt wedge which leads to complete flushing of seawater from the estuary. The well defined salt wedge is then regenerated during the subse- quent flood. Although we also see a well defined salt wedge at high discharge our observations at low freshwater discharge suggest that during the ebb mixing caused by both shear instabilities at the pycnocline and bottom generated turbulence is not able to homogenize the water column. We therefore expect the presence of mixed water in the estuary at the beginning of the subsequent flood. The presence of this mixed water will prevent the formation of a well defined salt wedge and the estuary will remain in a partially mixed state. Wave Height Unlike KH instabilities, the deflection of the density interface (wave amplitude) caused by one- sided instabilities is usually smaller than the amplitude of the billows (see figure 4.12). It is therefore difficult to assess the amplitude of the instabilities using echosoundings (e.g. figure 4.5). Chapter 4. Shear instability in the Fraser River estuary 64 Nevertheless, taking the vertical distance between the trough and the cusp, the observed instabili- ties vary in height from approximately 0.5 m to 2 m. The maximum height to wavelength aspect ratio of the observed instabilities varies between approximately 0.025 (0.5/20, figure 4.6c) and 0.1 (2/20, figure 4.5c). In the tilting tube experiments of Thorpe (1973) the maximum aspect ratio of KH instabilities varied between 0.05 and 0.6. Given the low values of J (< 1/4) in Thorpe’s experiments this difference in aspect ratio is not surprising. Unfortunately, other than the case of the KH instability (symmetric density and velocity profiles and J < 1/4) the height of shear instabilities in stratified flows is not well documented. Use of echosoundings to identify instability Our analysis focused on periods when instabilities were evident in the echosoundings. There were instances where predictions from the TG equation suggested instabilities would occur, but none were visible in the echosounder. In some cases (e.g. the secondary mode in figure 4.5), the lack of apparent instabilities in the echosoundings can be explained by the absence of the strong variations in salinity and temperature (i.e. density stratification) that are responsible for most of the back scatter of sound to the instrument (see Seim, 1999; Lavery et al., 2003, for a thorough description of acoustic scattering in similar environments). The quality of the visualization of the instabilities also depends on the speed of the boat relative to the speed of the instabilities. For profiles collected at 2.2 km, during transect 3 (figure 4.3c), the TG equation predicted instability close to the depth of the pycnocline (results not shown). In this region the boat speed and predicted instability speed were almost the same (-0.28 m s−1 versus -0.24 m s−1). Considering equation 4.3, the resulting apparent wavelength would be 250 m. The corresponding apparent period of approximately 15 minutes (250 m / -0.28 m s−1) would likely distort the appearance of an instability beyond recognition. This highlights an important challenge in identifying instabilities in echosoundings: if the point of observation is moving at a speed similar to the instabilities, the appearance of the instabilities becomes greatly distorted. On the other hand, if the observer is moving at a much different velocity than the instabilities, i.e. the apparent wavelength and period are relatively short, the sampling rate of the echosounder may not be sufficient to resolve the instabilities. In addition, our ability to detect shear instabilities depends on the timing of the echosounding relative to the stage of development of the instability. In DNS of symmetric and asymmetric instabilities there are several stages of development beginning with rapid growth and finishing with a breakdown into three-dimensional turbulence. Only during the stage where the instabilities have large two dimensional structures, e.g. billows, will they be easily recognizable with the echosounder. For example, during transect 2 (figure 4.3b) instabilities were only recognizable over a distance of approximately 100 m (11 to 11.1 km). However it is possible that instabilities are at a less recognizable stage of development throughout most of this transect. Chapter 4. Shear instability in the Fraser River estuary 65 Because of these challenges, the echosounder is able to confirm only the presence and not the absence of shear instability. We therefore limited our application of the TG equation to cases where instabilities were apparent. 4.8 Conclusions After performing a detailed stability analysis on six sets of velocity and density profiles using the Taylor-Goldstein equation and comparing with the echosoundings we conclude the following. 1. All of the instabilities observed in the echosoundings coincided with the most unstable mode in the TG analysis. This confirms the applicability of the TG equation in predicting instability, even in cases as complex as the Fraser River estuary. 2. The location of each of the observed instabilities occurs in a region of depth where Ri < 1/4. However, there are also cases that haveRi < 1/4 in which no unstable modes were ob- served. This result is in full agreement with the Miles-Howard criterion, but also highlights the inconclusive nature of this criterion. 3. Although the observed instabilities all act on a well defined density interface, they appear to be concentrated on only one side of the interface. The maximum vertical displacement occurs either above or below the density interface in a region of z where Ri < 1/4. None of the observations show Ri < 1/4 across the thickness of a density interface. This is in contrast to the archetypal KH instability described by the simple mixing layer model, in which Ri < 1/4 where ρ̄′ (N2) is greatest. The observed instabilities might therefore be better described by the so-called ‘one-sided’ modes of Lawrence et al. (1991); Carpenter et al. (2007). 4. During the majority of the survey the observed mixing was due to shear instabilities at the pycnocline. In other stratified estuaries with moderate to strong tidal forcing, such as the Columbia and Hudson rivers, turbulence generated at the bottom is considered the dominant source of mixing (Nash et al., 2008; Peters and Bokhorst, 2000). In the present study we only observe mixing due to bottom generated turbulence during late ebb at low freshwater discharge. 66 Bibliography Alexakis, A. (2005). On Holmboe’s instability for smooth shear and density profiles. Phys. Fluids, 17:084103. Browning, K. (1971). Structure of the atmosphere in the vicinity of large-amplitude Kelvin- Helmholtz billows. Quart. J. Roy. Met. Soc., 97:283–298. Browning, K., Bryant, G., Starr, J., and Axford, D. (1973). Air motion within Kelvin-Helmholtz billows determined from simultaneous Dopler radart and aircraft measurements. Quart. J. Roy. Met. Soc., 99:608–618. Carpenter, J., Lawrence, G., and Smyth, W. (2007). Evolution and mixing of asymmetric Holm- boe instabilities. J. Fluid Mech., 582:103–132. Farmer, D. and Armi, L. (1998). The flow of Atlantic water through the Strait of Gibraltar. Prog. Oceanogr., 21:1–98. Geyer, W. and Farmer, D. (1989). Tide-induced variation of the dynamics of a salt wedge estuary. J. Phys. Oceanogr., 19:1060–1672. Geyer, W. and Smith, J. (1987). Shear instability in a highly stratified estuary. J. Phys. Oceanogr., 17:1668–1679. Goldstein, S. (1931). On the stability of superposed streams of fluids of different densities. Proc. R. Soc. Lond. A, 132:524–548. Hazel, P. (1972). Numerical studies of the stability of inviscid stratified shear flows. J. Fluid Mech., 51:39–61. Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geofys. Publ., 24:67–112. Howard, L. (1961). Note on a paper of John W. Miles. J. Fluid Mech., 10:509–512. Kostachuk, R. and Atwood, L. (1990). River discharge and tidal controls on salt-wedge posi- tion and implications for channel shoaling: Fraser River British Columbia. Can. J. Civil Eng., 17:452–459. Lavery, A., Schmitt, R., and Stanton, T. (2003). High-frequency acoustic scattering from tur- bulent oceanic microstructure: The importance of density fluctuations. J. Acoust. Soc. Am., 114:2685–2697. Bibliography 67 Lawrence, G., Browand, F., and Redekopp, L. (1991). The stability of a sheared density interface. Phys. Fluids, 3(10):2360–2370. MacDonald, D. (2003). Mixing processes and hydraulic control in a highly stratified estuary. PhD thesis, Mass. Inst. of Technol. MacDonald, D. and Horner-Devine, A. (2008). Temporal and spatial variability of vertical salt flux in a highly stratified estuary. J. Geophys. Res., 113. Miles, J. (1961). On the stability of heterogeneous shear flows. J. Fluid Mech., 10:496–508. Miles, J. (1963). On the stability of heterogeneous shear flows. part 2. J. Fluid Mech., 16:209– 227. Moum, J., Farmer, D., Smyth, W., Armi, L., and Vagle, S. (2003). Structure and generation of turbulence at interfaces strained by internal solitary waves propagating shoreward over the continental shelf. J. Phys. Oceanogr., 33:2093–2112. Nash, J., Kilcher, L., and Moum, J. (2008). Turbulent mixing in the columbia river estuary: Structure and consequences for plume composition. J. Geophys. Res., page submitted. Pawlak, G. and Armi, L. (1998). Vortex dynamics in a spatially accelerating shear layer. J. Fluid Mech., 376:1–35. Peters, H. and Bokhorst, R. (2000). Microstructure observations of turbulent mixing in a partially mixed estuary. Part 1: Dissipation. J. Phys. Oceanogr., 30:1232–1244. Seim, H. (1999). Acoustic backscatter from salinity microstructure. J. Atmos. Ocean. Technol., 16:1491–1498. Seim, H. and Gregg, M. (1994). Detailed observations of naturally occurring shear instability. J. Geophys. Res., 99:10,049–10,073. Smyth, W., Carpenter, J., and Lawrence, G. (2007). Mixing in symmetric Holmboe waves. J. Phys. Oceanogr., 37:1566–1583. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Taylor, G. (1931). Effect of variation in density on the stability of superposed streams of fluid. Proc. R. Soc. Lond. A, 132:499–523. Tedford, E., Pieters, R., and Lawrence, G. (2009). Holmboe instabilities in a laboratory exchange flow. J. Fluid Mech., page submitted. Bibliography 68 Thorpe, S. (1973). Experiments on instability and turbulence in a stratified shear flow. J. Fluid Mech., 61:731–751. Thorpe, S. (2005). The Turbulent Ocean. Cambridge University Press, first edition. Ward, P. (1976). Seasonal salinity changes in the Fraser River Estuary. Can. J. Civil Eng., 3:342–348. Yonemitsu, N. (1991). The stability and interfacial wave phenomena of a salt wedge flow. PhD thesis, University of Alberta. Yonemitsu, N., Swaters, G., Rajaratnam, N., and Lawrence, G. (1996). Shear instabilities in arrested salt-wedge flows. Dyn. Atmos. Oceans, 24:173–182. Yoshida, S., Ohtani, M., Nishida, S., and Linden, P. (1998). Mixing processes in a highly strat- ified river. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 389–400. American Geophysical Union. Zhu, D. and Lawrence, G. (2001). Holmboe’s instability in exchange flows. J. Fluid Mech., 429:391–409. 69 Chapter 5 Unstable modes in asymmetric stratified shear layers 1 5.1 Introduction When considering the possibility of turbulence production and mixing in a sheared density strati- fied environment, it is important to determine whether or not a particular flow configuration repre- sents a stable solution of the equations of motion. This problem has a long history that started with the work of Helmholtz (1868) and Kelvin (1871) on the stability of homogeneous and stratified vortex sheets in the 19th century. Since then, numerous authors have examined cases of increasing complexity in an attempt to understand the basic instability mechanisms that are present in flows with greater physical relevance. Significant progress was made by Rayleigh (1880) who examined a piecewise-linear repre- sentation of the homogeneous shear layer, represented by U(z), that consists of a finite shear thickness h (see figure 5.1). A stability analysis on this idealized shear layer produced results that are in qualitative agreement with subsequent studies of smooth profiles, like the hyperbolic tangent shear layer sketched in figure 5.1 (Michalke, 1964). This suggests that the piecewise shear layer is sufficient to capture the basic instability mechanism that is also present in more realistic smooth profiles. Motivated by geophysical flows, Holmboe (1962) extended Rayleigh’s analysis to include a stable density stratification. The piecewise-linear shear layer was retained and a lay- ered, piecewise-constant density profile ρ̄(z) was added with a step change in density of ∆ρ at the shear layer centre (figure 5.1, with d = 0). As in the homogeneous (Rayleigh) case, these ideal- ized profiles give qualitatively similar results to the smooth profiles shown in figure 5.1 (Hazel, 1972). In representing Holmboe’s (1962) piecewise profiles with smooth functions, it is necessary to introduce another length scale, given by the density interface thickness δ. Agreement between the smooth and piecewise profiles is achieved when δ is sufficiently small (i.e. R ≡ h/δ & 3 in general) with the piecewise profiles representing the limit of vanishing δ. Holmboe’s (1962) stability analysis shows the presence of two distinct types of unstable modes. Which mode of instability develops is found to be dependent on the wavenumber k, made dimensionless by α ≡ kh/2, and the relative strength of the stratification, measured by the dimensionless bulk Richardson number J ≡ g′h/(∆U)2, where g′ = ∆ρg/ρ0 is the reduced gravity, and ρ0 a reference density. The resulting stability diagram is shown in figure 5.2(a) for smooth ‘tanh’ profiles with R = 5. At low J , the stratification is relatively unimportant, and 1A version of this chapter has been submitted for publication. J.R. Carpenter, N.J. Balmforth and G.A. Lawrence (2009) Identifying unstable modes in stratified shear layers. Chapter 5. Unstable modes in asymmetric stratified shear layers 70 ΔU/2 ΔU/2 h/2 h/2 z U(z) z ρ0+Δρ/2ρ0 d δ ρ(z) Figure 5.1: Profiles of the stratified shear layers to be considered. The thick lines indicate the smooth profiles and the thin lines represent the piecewise profiles. the resulting mode of instability is essentially a stratified analogue of the Rayleigh instability. It is common in the literature to refer to this mode as Kelvin-Helmholtz (KH), despite the closer association with Rayleigh’s shear layer, and we keep with this convention throughout. When J is sufficiently large, the qualitative behaviour of the instability changes. Unstable Holmboe (H) modes develop due to a destabilizing influence of the stratification (Holmboe, 1962; Howard and Maslowe, 1973), reaching a peak growth rate at finite J . This distinct change in the stability properties that occurs across the KH-H transition can also be seen in the nonlinear development of the instabilities (Smyth et al., 1988; Smyth and Peltier, 1991; Hogg and Ivey, 2003). Figure 5.2(b-d) shows plots of the density field for KH and H in- stabilities, taken from a series of direct numerical simulations2 (DNS), once they have reached a large amplitude nonlinear stage of development. The KH instability (figure 5.2b) exhibits the well-known billows of overturning fluid caused by a rolling up of the shear layer vorticity (as in the homogeneous shear layer). These billows become susceptible to secondary instabilities and subsequently breakdown to drive turbulent mixing of the density field (Thorpe, 1973; Caulfield and Peltier, 2000). The nonlinear form of H instabilities (figure 5.2c,d) consist of cusp-like prop- agating waves that protrude into the upper and lower layers. These generally do not involve a complete overturning of the density interface. It is perhaps not surprising that this different finite amplitude behaviour between the KH and H modes has recently been found to have a pronounced effect on the mixing of mass and momentum in the shear layer; in some cases changing the ef- 2The DNS results in this chapter were performed using the code of Winters et al. (2004). A two-dimensional grid was utilized as well as periodic boundary conditions in x and free-slip on the top and bottom. Since the model is viscous and diffusive, moderate values of Re = ∆Uh/ν = 1200 and Pr = ν/κ = 25, with ν the kinematic viscosity and κ the diffusivity, were used. These do not result in significant changes in stability properties compared to the inviscid nondiffusive theory. Chapter 5. Unstable modes in asymmetric stratified shear layers 71 α J 0.0 6 0. 04 0. 02 0.16 0.1 2 0. 06 0.0 4 0.25 0.5 0.75 1 0 0.1 0.2 0.3 (a) b c d (b) (c) (d) 0.2 0.3 + + 0.4+ 0.5+ Figure 5.2: Stability diagram (a) of the Holmboe model of the stratified shear layer for smooth (tanh) profiles with R = 5. Dark contours are of growth rate and grey contours of phase speed with the dark grey shading representing regions of stability. The thick lines correspond to the stability boundaries and the transition between stationary (below) and propagating (above) modes. Representative density fields from DNS are shown in (b-d). The approximate location of the instabilities on the stability diagram is indicated with letters. The different widths of (b-d) are due to changes in the wavenumber of maximum growth rate. The vertical domain height is taken to be 10h for the stability diagram and in the DNS, which is sufficiently large to approximate unbounded domains. fective diffusivity by an order of magnitude across the KH-H transition (Smyth et al., 2007). Not only the amount of mixing, but also the character of the mixing have been found to depend on the resulting mode type (Carpenter et al., 2007). It is therefore important to predict which mode of instability is to occur when quantifying mixing and momentum transfers. An important feature of Holmboe’s profiles is that they exhibit a symmetry about the shear layer centre (once the Boussinesq approximation has been made). However, observations of strat- ified shear instabilities in the field often display some asymmetry between the shear layer centre and the vertical location of the density interface (Armi and Farmer, 1988; Wesson and Gregg, 1994; Yoshida et al., 1998, chapter 4). It is also common to observe this asymmetry in laboratory experiments (Koop and Browand, 1979; Lawrence et al., 1991; Pawlak and Armi, 1998). The im- plications that this asymmetry has on the stability of the flow was first studied by Lawrence et al. (1991) using the piecewise model of Holmboe (1962) with a density interface located a distance Chapter 5. Unstable modes in asymmetric stratified shear layers 72 J 0 .1 6 0 .1 2 0 .0 8 0 .0 4 0 .1 2 0 .0 8 0 .0 4 0 0.1 0.2 0.3 α 0.5 1 1.5 (a) b c d (b) (c) (d) 0. 3 0. 2 0.1 Figure 5.3: Stability diagram of the asymmetric stratified shear layer for smooth profiles with R = 5 and a = 0.5. See caption of figure 5.2 for further details. d below the shear layer centre, as shown in figure 5.1. They found that two distinct branches of instability were present, each consisting of propagating modes. No distinct transition from KH to H modes is apparent. One of the two modes always consists of larger growth rates than the other (referred to as the dominant mode), and has often been found to be the only mode observed at large amplitudes (Lawrence et al., 1991; Haigh, 1995). Stability results for the dominant mode of the asymmetric stratified shear layer using smooth profiles with R = 5, and the asymmetry parameter a ≡ 2d/h = 0.5, are shown in figure 5.3. No distinct transition between KH and H modes is apparent from the results of linear stability theory when an asymmetry is present (a 6= 0). The asymmetric laboratory observations of Lawrence et al. (1991, 1998) show a continuous change in behaviour from overturning billows to cusp-like waves. This can also be seen in the simulation results shown in figure 5.3(b-d), where instabilities resemble KH at low J , and become more like H instabilities at larger J . The lack of any distinct transition in the stability properties when the flow is asymmetric demonstrates the difficulty in distinguishing between KH- and H-like instabilities in this case, and raises the question of how to appropriately define what is meant by KH and H modes. When the flow is perfectly symmetric the marginal curve separating stationary (cr = 0) modes from propagating (cr 6= 0) modes also coincides well with distinct changes in growth rate. These Chapter 5. Unstable modes in asymmetric stratified shear layers 73 changes in stability characteristics have been found to match reasonably well with changes in the nonlinear behaviour of the instabilities, with stationary billows occurring in the KH (cr = 0) region and cusp-like propagating waves in the H (cr 6= 0) region. It should be noted however, that nonlinear behaviour characteristic of both types of instabilities has been observed near the transition between the two mode types (Smyth and Peltier, 1991; Hogg and Ivey, 2003). This led Smyth and Peltier (1991) to hypothesize that nonlinear effects can modify the location of the transition. The fact that a transition between KH and H behaviour, based on the nonlinear dynamics, is observed in asymmetric instabilities, suggests that we must reexamine the distinction based on phase speed between the two modes. We therefore develop a diagnostic to interpret the unstable modes of asymmetric stratified shear layers. The purpose of this chapter is to utilize linear theory to predict the occurrence of KH- and H-type modes in asymmetric flows where the distinction between these two modes is blurred. The motivation for this is based on three previous findings: 1. KH and H instabilities result from two different linear growth mechanisms (to be discussed fully in the next section) (Holmboe, 1962; Baines and Mitsudera, 1994; Caulfield, 1994); 2. it is common to find asymmetry in geophysically relevant flows (Armi and Farmer, 1988; Wesson and Gregg, 1994; Yoshida et al., 1998, chapter 4); and 3. the type of instability that develops can have a significant influence on turbulent mixing and vertical transports (Smyth and Winters, 2003; Smyth et al., 2007; Carpenter et al., 2007). We use results from a linear stability analysis of both piecewise and smooth profiles of the asym- metric stratified shear layer to predict the occurrence of KH and H modes. These results are based on the ‘wave interaction’ interpretation of shear instability that is reviewed in the following sec- tion. This is followed by the formulation of a diagnostic in §5.3 that is used to distinguish between the contributions that KH and H modes make to the stability of the flow. Results of applying the diagnostic to the symmetric (a = 0) and asymmetric (a 6= 0) cases shown above are outlined and discussed in §5.4. We then examine a set of profiles measured from a highly stratified estuary in order to illustrate the applicability of the formulation in a geophysical context. Conclusions are stated in the final section. 5.2 Wave interaction interpretation of instability The wave interaction interpretation attributes instability in stratified and homogeneous shear flows to a mutual interaction between otherwise freely propagating stable waves in the profiles. The majority of work that has utilized the wave interaction interpretation has been carried out on the idealized piecewise -linear representation of the stratified shear layer (Holmboe, 1962; Cairns, 1979; Caulfield, 1994; Redekopp, 2001), since this is the easiest possible geometry to understand Chapter 5. Unstable modes in asymmetric stratified shear layers 74 and apply the theory. Piecewise profiles of U and ρ̄ are particularly simple because they have delta function behaviour of vorticity gradients U ′′, and density gradients, represented by N2 = −gρ̄′/ρ0, where primes denote differentiation with respect to z. At these locations, referred to as interfaces, wave motion may occur. The phase speed of waves on vorticity and density interfaces in a stationary frame of reference is given by cv = ∆q 2k and cd = ± ( g′ 2k )1/2 , (5.1) respectively, where ∆q = U ′(z+v ) − U ′(z−v ) denotes the jump in vorticity across the vorticity interface. Note that the vorticity interface supports a single unidirectional mode of propagation, whereas the density interface supports two oppositely propagating modes. The wave interaction interpretation requires that two interfaces must be present, each sup- porting an oppositely propagating wave mode, in order for instability to be possible. The two interfacial waves are then able to interact such that (i) they are stationary relative to one another, and (ii) in a ‘phase-locked’ position such that they may cause mutual growth. It is only possi- ble for (i) to occur between two oppositely propagating wave modes when there is shear in the background profile. Condition (i) suggests that the region of instability in the αJ-plane should be close to the locus of points where the two freely propagating interfacial wave modes have equal phase speeds. Although this is not strictly true, as each wave will interact and adjust the others phase speed, in the limit of large α this interaction vanishes and the approximation becomes ac- curate. This large-α approximation has proven useful in identifying different instability modes in previous studies (Caulfield, 1994; Baines and Mitsudera, 1994; Redekopp, 2001). This method of identifying the wave interactions that lead to instability provides an interpre- tation of the KH and H modes that are observed in the piecewise symmetric stratified shear layer profiles of Holmboe (1962). Figure 5.4(a) shows the resulting stability diagram with the large-α approximation shown as a dashed line. The curve is obtained by equating the upper (lower) vor- ticity wave speed with the rightward (leftward) propagating internal wave speed on the density interface. The close agreement between the large-α approximation and the region of instabil- ity indicates that the H modes are caused by an interaction between vorticity and internal wave modes. As J vanishes, we approach Rayleigh’s homogeneous shear layer, where instability must result from the interaction of the two vorticity modes. This was clearly illustrated by Baines and Mitsudera (1994) in considering the same profiles as Holmboe (1962) except with the lower vor- ticity interface removed, shown in figure 5.4(b). By comparing the stability diagrams in figure 5.4, we see that the H region of instability remains largely unchanged, however, the KH region is completely eliminated. Since only two interfacial waves are able to interact – the upper vorticity wave and the rightward propagating internal wave – only one type of unstable mode (the H mode) is present. The KH instability is not present since it is caused by the interaction of the upper and lower vorticity modes. Chapter 5. Unstable modes in asymmetric stratified shear layers 75 α 0 1 2 3 4 z U(z) ρ1 ρ2 > ρ1 α J 0 1 0 0.3J 0 1 2 α 0 1 2 3 4 α J 0 1 0 0.3 z U(z) ρ1 ρ2 > ρ1 (a) (b) Figure 5.4: Stability diagram after Baines and Mitsudera (1994) showing (a) the piecewise sym- metric stratified shear layer of Holmboe (1962), and (b) the profiles of Baines and Mitsudera (1994). Contours are of growth rate with the grey shading representing regions of stability. Thick dashed line represents the large-α approximation obtained by equating the speeds of an internal wave and a vorticity wave. Profiles are shown as insets in the upper left. The dotted lines denote the boundary of the close-up regions shown as insets in the lower right. In the above description, we have only concentrated on piecewise profiles. However, the same mechanisms are still believed to apply to smooth profiles. In this case, rather than having delta function behaviour of U ′′ and N2 at the interfaces, those functions take on smooth distributions that attain extrema in an ‘interfacial region’. The KH instability is now a result of the inflection point in the U -profile that separates two regions of oppositely signed vorticity gradients. Like- wise, the H instability is the interaction of the region of strong vorticity gradients (U ′′) with a strong density gradient region (N2). It does not require the presence of an inflection point. Note that this is not a violation of Rayleigh’s (1880) inflection point theorem since it applies only to homogeneous flows. Similarities between smooth profiles and piecewise profiles in terms of wave interactions has been discussed previously by Baines and Mitsudera (1994), and will also be seen in the results to follow. Stratified shear layers consisting of two density interfaces may also be susceptible to a third instability type that was first described by Taylor (1931). These unstable Taylor (T) modes re- sult from the interaction of two oppositely propagating waves on the density interfaces that may become phase-locked due to a background shear (Caulfield, 1994). Similar to the H modes, the T modes do not require the presence of an inflection point (in fact U ′′ may be identically zero throughout the domain). In general, the instability of a stratified shear layer may be described in terms of these three interaction types (i.e. KH, H, and T). In the following section we formulate a diagnostic in order to quantify the strength of the three types of wave interactions. Chapter 5. Unstable modes in asymmetric stratified shear layers 76 5.3 Formulation of a diagnostic In this section we utilize condition (ii) from §5.2, that the interacting waves must cause mutual growth in each other, to formulate a diagnostic used to interpret unstable modes of the stratified shear layer. The formulation is general, and may be applied to any profiles in which distinct interfaces can be identified. This allows for a classification of the unstable modes in terms of KH-, H- and T-types, following Caulfield (1994). 5.3.1 Taylor–Goldstein equation We will be concerned with the small amplitude motions of an incompressible inviscid Boussinesq fluid, with perturbations taken about the basic profiles that are small enough to be well approxi- mated by the linearized equations of motion. Following the framework of Holmboe (1962), we partition the total perturbation vorticity of the flow, q, into a kinematic portion qK , and a baroclinic portion, qB . The kinematic vorticity is created by the vertical displacement of vorticity gradients in the U -profile, and is given in the linear approximation by qK = −U ′′η, (5.2) where η denotes the vertical displacement field. In Boussinesq fluids, baroclinic vorticity is pro- duced by the horizontal tilting of the constant density surfaces of the ρ̄-profile. This baroclinic production of vorticity may be written, within the linear approximation, as DqB Dt = N2 ∂η ∂x , (5.3) where the material derivative here and throughout the remainder of the paper has been linearized, and is given by D Dt ≡ ∂ ∂t + U ∂ ∂x . The total perturbation vorticity can now be written as the sum of the kinematic and baroclinic portions viz. q = ∇2ψ = qK + qB. (5.4) Here we have used a stream function representationψ of the perturbation velocity field u = (u,w), such that u = ∂ψ/∂z andw = −∂ψ/∂x. This ensures that the incompressible continuity equation for the perturbation velocity field is satisfied. Changes in the vertical displacement field may be related to the vertical velocity through the kinematic condition Dη Dt = −∂ψ ∂x , (5.5) Chapter 5. Unstable modes in asymmetric stratified shear layers 77 and allows the problem, given by (5.2) through (5.4), to be expressed in terms of a single equation for the stream function. Perturbations can now be taken to be of the normal mode form, i.e. ψ(x, z, t) = ψ̂(z)eik(x−ct), (5.6) for the stream function, where k is the horizontal wave number and c = cr + ici is the complex phase speed. Substituting this form results in the well known Taylor–Goldstein (TG) equation ψ̂′′ + [ N2 (U − c)2 − U ′′ U − c − k 2 ] ψ̂ = 0. (5.7) This equation, together with the condition that ψ̂ vanishes on the boundaries (which may be taken at z = ±∞), describes an eigen-problem for the eigenvalue c, and the eigenfunction ψ̂(z). The flow given by the basic profiles is said to be unstable if the eigenvalue has ci > 0, which indicates that the perturbations grow at an exponential rate σ = kci. 5.3.2 Partitioning into kinematic and baroclinic fields Once the normal modes are determined from solving the TG equation, it is possible to use this information to examine the roles that the kinematic and baroclinic fields play in the growth of the resulting instabilities. This may be accomplished by first noting that the ψ due to the kinematic vorticity alone, may be determined directly from η. By defining this stream function field as ψK(x, z, t), (5.4) implies that ψ = ψK + ψB, and similarly ψ̂ = ψ̂K + ψ̂B, (5.8) so that (5.2) may be written as ψ̂′′K − k2ψ̂K = −U ′′η̂. (5.9) A similar form may be found for the baroclinic field from (5.3), and expressed as ψ̂′′B − k2ψ̂B = N2 U − c η̂. (5.10) Since η̂ and c are both known from the solution to the TG equation, (5.9) and (5.10) can be expressed in the general form L[ψ̂] = f(z), (5.11) where the linear operator L = d2/dz2 − k2, and f(z) may be regarded as some known forcing function given by the right hand sides of (5.9, 5.10). We may now solve for ψ̂(z) by inverting the Chapter 5. Unstable modes in asymmetric stratified shear layers 78 operator L by the relation ψ̂ = ∫ D G(s, z)f(s)ds, (5.12) where D is the domain, and G(s, z) is the appropriate Green’s function for L, which depends on the type of domain and boundary conditions. For an unbounded domainG(s, z) = −e−k|z−s|/2k, otherwise G(s, z) = −1 k sinh[k(zu − zl)] · { sinh[k(zu − z)] sinh[k(s− zl)] , z > s sinh[k(zu − s)] sinh[k(z − zl)] , z < s for domains with upper and lower boundaries at zu and zl, respectively. From (5.12) it is possible to partition the ψ̂ into kinematic and baroclinic effects. The contri- bution of each field to the growth rate and phase speed of the normal mode disturbance can now be explicitly solved for by rearranging the kinematic condition (5.5) and substituting the normal mode form (5.6), to give cr = U + Re ( ψ̂K + ψ̂B η̂ ) and σ = kIm ( ψ̂K + ψ̂B η̂ ) . (5.13) The relation for σ in (5.13) will be used throughout the remainder of the chapter to assess the contributions of the kinematic and baroclinic fields in the growth of unstable modes. 5.3.3 Piecewise profiles We now consider the idealized piecewise profiles in which the vorticity and density gradients (U ′′ and N2) exhibit delta function behaviour at a number of discrete interface locations. This allows us to write U ′′(z) = n∑ j=1 ∆qj δ(z − zj) and N2(z) = m∑ `=1 g′` δ(z − z`) (5.14) where we are considering general piecewise profiles consisting of n vorticity interfaces with jumps ∆qj , and m density interfaces with jumps g′`, at the vertical locations zj and z`, respectively. When the delta-function forms in (5.14) are substituted into (5.12) the integrals for ψ̂K and ψ̂B reduce to sums, where each term represents the contribution of a particular interface. If we now choose an interface of interest, the pth interface say, and apply (5.13), we are able to break the total growth rate of the normal mode σ, into the individual contributions of each interfacial wave. This allows us to write σ = n∑ j=1 σpKj + m∑ `=1 σpB`, (5.15) where each term of the sums will be referred to as a partial growth rate, and are given explicitly Chapter 5. Unstable modes in asymmetric stratified shear layers 79 by σpKj = −k Im { ∆qj G(zj , zp) η̂(zj) η̂(zp) } and σpB` = k Im { g′` G(z`, zp) U(z`)− c η̂(z`) η̂(zp) } . (5.16) Note that for piecewise profiles a vorticity interface cannot cause growth in itself, i.e. σpKp = 0, since ∆qp is a real number. 5.3.4 Smooth profiles The same partial growth rate diagnostic is now developed for continuous distributions of U ′′ and N2. In doing so we presume that the domain can be split into a number of ‘interfacial regions,’ where either U ′′ or N2 reaches some extrema. A ψ̂K,B can then be defined for each interface, using (5.12), as ψ̂Kj(z) = − ∫ Dj G(s, z)U ′′(s)η̂(s) ds (5.17) and ψ̂B`(z) = ∫ D` G(s, z) N2(s) U(s)− c η̂(s) ds, (5.18) where the Dj,` denotes the domain of the jth kinematic vorticity region and the `th baroclinic vorticity region. Using (5.13) we can write σ = n∑ j=1 piKj(z) + m∑ `=1 piB`(z), (5.19) with piKj(z) = k Im(ψ̂Kj/η̂) and piB`(z) = k Im(ψ̂B`/η̂). Finally, to apply this condition to the pth interfacial region, we multiply both sides of (5.19) by a suitable weight function F (z), integrate over Dp, and rearrange to give a direct analogy to (5.15) for smooth profiles, σ = n∑ j=1 〈FpiKj〉p 〈F 〉p + m∑ `=1 〈FpiB`〉p 〈F 〉p (5.20) = n∑ j=1 σpKj + m∑ `=1 σpB` (5.21) where 〈〉p indicates integration overDp. A natural choice for the weight function is either F = U ′′ if p corresponds to a vorticity interface, or F = N2 if p is a density interface, and we will keep with this convention throughout. Chapter 5. Unstable modes in asymmetric stratified shear layers 80 5.3.5 Classification of modes We have described above, a method whereby the growth rate of each interface σ, can be expressed as a sum of the contributions from all of the other interfaces present. This allows us to quantify whether or not a particular growth mechanism is present if two interfaces are causing mutual growth in one another. If the interfaces are both vorticity interfaces interacting across an inflection point in the U profile, then the growth mechanism is classified as KH-type, for example. The other two possible mechanisms (following Caulfield, 1994) can be due to the interaction of a vorticity- internal wave (i.e. of the H-type), or an internal wave-internal wave (i.e. of the T-type). However, for reasons of brevity we shall limit our application of the partial growth rates in quantifying the strength of only KH and H mechanisms in the following sections. 5.4 Results We now apply the partial growth rate formulation for piecewise and smooth profiles to the stratified shear layers in figure 5.1. The symmetric case is examined first, followed by the asymmetric case where we let a = 0.5. This particular value for a was chosen to match the asymmetry that Lawrence et al. (1991) observe in their laboratory experiments. In the case of smooth profiles, we will take U(z) = ∆U 2 tanh (2z h ) and ρ̄(z) = ρ0 − ∆ρ2 tanh (2Rz h + a ) , (5.22) and examine a single value of the interfacial thickness ratio R = 5, which is large enough to permit a significant region of unstable propagating modes when a = 0. The profiles consist of two vorticity interfaces and a single density interface. Therefore, the T mode that results from the interaction of two density interfaces can immediately be disregarded. The partial growth rate diagnostics in (5.15) or (5.20) reduces to a sum of three terms for each interface. However, we will limit our attention to only the rightward propagating H mode, the KH mode, as well as the dominant asymmetric mode. In doing so, it will suffice to apply the partial growth rate diagnostic only to the upper vorticity interface, since each unstable mode will consist of mutual growth between the upper vorticity interface, and one of the lower vorticity or density interfaces. The diagnostic equation then becomes σ = σKH + σH + σself , (5.23) where the terms on the right hand side represent the partial growth rate due to the interaction between the upper vorticity interface and (from left to right) the lower vorticity, density, and upper vorticity interfaces. These terms are identified with KH, H, and self-interaction components, respectively. Chapter 5. Unstable modes in asymmetric stratified shear layers 81 J (a) σ 0 .1 2 0 .0 8 0. 04 0. 160. 12 0. 08 0 0.1 0.2 0.3 0.12 0.2 0.08 0.16 -0.04 -0.08 0 0 .0 8-0.08 α 0.25 0.5 0.75 α 0.25 0.5 0.75 α 0.25 0.5 0.750 00 (b) σKH (c) σH Figure 5.5: Partial growth rates for the symmetric case using piecewise profiles. Contours of σ are shown in (a), which is given as the sum of σKH in (b), and σH in (c). The thick solid line denotes the stability boundary and transition from stationary (cr = 0) to propagating (cr 6= 0) modes. Dark grey shading indicates a region of stable modes, whereas light grey shading denotes regions where the partial growth takes negative values. 5.4.1 Symmetric profiles Beginning with the stability properties of the piecewise symmetric stratified shear layer, we plot the partial growth rates from (5.23) in figure 5.5. Recall that a vorticity interface cannot interact with itself, therefore, for all the piecewise results σself = 0 and only σKH and σH are required in the sum, i.e. the normal mode growth rate can be expressed as a KH and H component according to (5.23). In the region of stationary (cr = 0) instability, growth is due entirely to σKH (figure 5.5b), while the density interface acts as a stabilizing influence, indicated by the negative values of σH (figure 5.5c). In addition, the asymptotic result of an unstable H mode consisting predominantly of a σH component is recovered for large α (and J), as expected. However, as J increases to- wards the transition to propagating modes, σKH increases to a maximum value on the transition. Once we cross over to the propagating region of the diagram σKH gradually asymptotes to zero. This behaviour indicates that the interaction of vorticity interfaces, normally associated with the stationary KH mode, does contribute to the growth of the propagating modes, particularly near the transition between the two. This KH contribution is greatest in the low-α portion of the propagat- ing modes where σH < 0 (figure 5.5c). Note that a similar result was observed in §5.2 when the stability diagram of Holmboe (1962) was compared to the profiles of Baines and Mitsudera (1994), consisting of only the upper vorticity interface and density interface (figure 5.4): removal of the lower vorticity interface eliminates the stationary modes but also alters the normal mode growth σ into the propagating region. In other words, despite the propagating modes being primarily an H-type vorticity-internal wave interaction, the presence of additional interfaces will affect changes in σ. This also reflects the Chapter 5. Unstable modes in asymmetric stratified shear layers 82 α J 0.0 60 .0 4 0. 02 0.160 .12 0. 06 0.0 4 0.25 0.5 0.75 1 0 0.1 0.2 0.3 α 0 .1 5 0 .1 0 .0 5 0 -0.05 -0.1 0.25 0.5 0.75 1 α 0.05 0.10.15 0.2 0.25 0.5 0.75 1 (a) σ (b) σKH (c) σH σ self Figure 5.6: Partial growth rates for the symmetric case using smooth profiles. All notation as in figure 5.5 except σself (not shown) must be included to complete the sum. fact, which has been observed in nonlinear studies of the transition (Smyth and Peltier, 1991; Hogg and Ivey, 2003), that the KH mechanism of vorticity wave interaction exerts an influence into the region of propagating modes. The same general behaviour is also seen in the partial growth rates of the smooth profiles shown in figure 5.6. Note however, the differences between the normal mode growth rate σ, in the piecewise and smooth profiles in figures 5.5(a) and 5.6(a). An important difference in obtaining solutions to the TG equation for smooth profiles is that it must be done numerically on a finite domain. A domain height of 10h was chosen since it is large enough to closely approximate solu- tions on an unbounded domain, however, it has the effect of destabilizing the low-α disturbances (Hazel, 1972). In addition, when examining smooth profiles, σself must be included in order to complete the sum in (5.23). Although σself can be a large term, it is always found to be a negative (stabilizing) influence, and does not significantly affect the interpretation of smooth profiles. For this reason we do not show the details of σself in the results. Just as in the piecewise results in figure 5.5 the continuous profiles show that the stationary modes are dominated by the KH component of the partial growth rates. This persists into the propagating region, particularly at small α, where σH continues to be negative. As J is increased the H interaction between shear layer vorticity and the density interface becomes the dominant mechanism of instability growth. We can again conclude that a significant portion of the propa- gating modes experience a KH-type mechanism of inflection point growth near the phase speed transition. Since finite amplitude instabilities in stratified shear layers have generally been found to occur at the wavenumber where the growth rate is a maximum, denoted by αmax, it is these disturbances which we should expect to be the most relevant. It is therefore of interest to compare the relative contributions of σKH and σH to the growth of the instability along the wavenumber of maximum growth curve. This is shown in figure 5.7 where the partial growth rates are plotted as a function of Chapter 5. Unstable modes in asymmetric stratified shear layers 83 (a) Piecewise profiles (b) Smooth profiles c r transition σKHσH σ 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 J g ro w th r at e 0 0.05 0.1 0.15 0.2 0.25 0.3 J (c) (d) (e) c d d e e σH = σKH Figure 5.7: Partial growth rates for the symmetric case along αmax. Piecewise results are shown in (a), while (b) gives the results for smooth profiles along with the numerical results illustrating the nonlinear form of the instabilities. The transition between stationary and propagating modes is given by the vertical grey line (with cr = 0 for small J and cr 6= 0 at larger J). The dashed grey line denotes where σKH = σH . Negative partial growth rates are not shown. J along the αmax-curve. Note that since αmax changes discontinuously across the transition from stationary to propagating modes, figure 5.7 shows that σKH and σH are also discontinuous there. There are two locations of particular interest in these plots: the transition from stationary to prop- agating modes (vertical grey line), and the point at which σKH = σH (vertical dashed grey line). The first location denotes the commonly accepted transition between KH and H modes based on phase speed considerations, and the second shows our diagnostic describing the predominance of either the KH- or H-type growth mechanisms. In the case of the piecewise profiles of figure 5.7(a) these two locations are coincident. However, the smooth profiles show that H-type growth does not predominate until J has surpassed the transition to propagating modes. The partial growth rate diagnostic appears to provide a good description of the nonlinear dynamics for the few simu- lations performed, as shown in figure 5.7(c-e). It is also interesting to note that Smyth and Peltier (1991) have observed nonlinear characteristics of KH instabilities within the propagating modes predicted by linear stability theory in their numerical simulations of the KH-H transition. 5.4.2 Asymmetric profiles Whereas the symmetric profiles show a distinct transition from stationary to propagating instabil- ities, there is no such transition in the case of asymmetric profiles. Instead, two unstable propa- gating (cr 6= 0) modes, each travelling in opposite directions with respect to the mean flow, are Chapter 5. Unstable modes in asymmetric stratified shear layers 84 0 .0 4 (a) σ (b) σKH (c) σH J α 0 .1 6 0 .1 2 0 .0 8 0 .1 2 0 0.25 0.5 0.75 0 0.1 0.2 0.3 α 0.16 0.12 0.08 0 0.25 0.5 0.75 α 0 0 .0 4 0 .0 8 0 0.25 0.5 0.75 Figure 5.8: Partial growth rates for the dominant mode of the asymmetric (a = 0.5) case using piecewise profiles. Details as in figure 5.5. present for all J > 0. However, we note that as J → 0 we recover Rayleigh’s homogeneous shear layer instability, and so KH behaviour is expected at small J . Also, for large J (and large α) it can be seen that both instability regions are in alignment with the resonance condition between the vorticity interfaces and the interfacial waves on the density interface (obtained by equating the phase speeds of each pair of waves, see §5.2). We therefore expect that the unstable modes are of H-type at large J , and some sort of transition between the two linear growth mechanisms must take place. This behaviour is reflected in the partial growth rates of the piecewise profiles shown in figure 5.8. Only the dominant mode consisting of larger growth rates is plotted since it has most often been found to control the development of the instabilities into the nonlinear regime (Lawrence et al., 1991; Haigh, 1995). Figure 5.8 shows that both σKH and σH are positive over the majority of the αJ-plane where the unstable mode is present. Since σKH is a decreasing function of J , and σH is an increasing function of J , we expect a gradual transition from KH- to H-type modes. In other words both KH and H growth mechanisms are present throughout the majority of the unstable region, and they replace one another in a gradual fashion as J is increased. The resulting plots for smooth profiles (figure 5.9) show qualitatively similar results. Some difficulty was encountered in the numerical solution of the TG equation for smooth profiles near the stability boundary due to the near-singular behaviour of the solutions at the critical level (see Alexakis, 2005, for further discussion). Evidence of this can be seen in figure 5.9(c), however, we do not expect this to have significantly affected the results. Focusing on the most amplified wave number, figure 5.10 shows a gradual transition from a growth dominated by σKH to one where σH becomes the stronger influence, as J is increased. This is in accord with previous observations of the nonlinear behaviour of the instabilities (Lawrence et al., 1991, 1998), and can also be seen in the results from the numerical simulations (figure 5.10c- e). A KH-like behaviour can be seen in the instability at J = 0.10, where significant overturning Chapter 5. Unstable modes in asymmetric stratified shear layers 85 J 0 .1 6 0 .1 2 0 .0 8 0 .0 4 0 .1 2 0 .0 8 0 .0 4 0 0.1 0.2 0.3 α 0.1 0.05 0.5 1 1.5 α 0 .1 5 0.2 0 .1 0 .0 50 .1 0 .0 5 0.5 1 1.5 α 0.5 1 1.5 (a) σ (b) σKH (c) σH σ self Figure 5.9: Partial growth rates for the asymmetric (a = 0.5) case using smooth profiles. Details as in figure 5.6. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 J gr ow th ra te (a) Piecewise profiles 0 0.05 0.1 0.15 0.2 0.25 0.3 J (b) Smooth profiles σKH σH σ (c) (d) (e) c c d d e e σH = σKH Figure 5.10: Partial growth rates for the asymmetric case along the wave number of maximum growth. Details as in figure 5.7. is present in the density field, and this is predicted in the larger value of σKH . By J = 0.30 we have a nonlinear form of the instability that is more akin to the H instability, consisting of a cusp-like wave. A result that can be seen directly from figure 5.10 is that the asymmetry extends the KH mechanism to larger values of J than in the symmetric case. 5.5 Application to field profiles We now give an example of how to apply the diagnostic toU and ρ̄ profiles collected from the field. The profiles shown in figure 5.11 were measured in the Fraser River estuary using the methods Chapter 5. Unstable modes in asymmetric stratified shear layers 86 described in chapter 4. Both of the U and ρ̄ profiles can be seen to consist of a number of ‘layers’ which lead to many possible interactions between interfaces. These interfaces can be identified by extrema in the profiles of U ′′ and N2 shown in figure 5.11(a) and 5.11(b), respectively. The presence of multiple inflection points in U , and a number of maxima in N2 indicate that any of the three instability types (KH, H, and T) mentioned in §5.2 may be present. A linear stability analysis of these profiles predicts an unstable mode with a maximum growth rate of σ = 0.024 s−1 at a wavenumber of k = 0.61 m−1 and a displacement eigenfunction amplitude |η̂|, shown in figure 5.11(c). These predictions of linear theory were found in chapter 4 to be in good agreement with observations recorded by an echosounder (figure 5.11d), both in the wavenumber k, and the vertical location given by the maximum of |η̂|. We now determine which interfacial wave interactions are responsible for generating the in- stabilities seen in the echosounding of figure 5.11(d). Due to the many interfaces that are present in the profiles, we shall simplify the analysis by noting that the dominant vertical location of the perturbations is concentrated near the peak of |η̂| at z = 10.9 m. This suggests that the instability is due only to the influence of the interfaces nearest this level. These consist of the two uppermost vorticity interfaces, with extrema at z ≈ 11.5 m and 10 m, as well as the density interface at z ≈ 10.5 m. Assuming this to be the case, we modify the observed profiles so that U ′′ and N2 vanish below, with U ′ remaining constant. The modified profiles will be denoted by an asterisk subscript (e.g. U∗), and are shown in figure 5.11. Performing a linear stability analysis on the modified profiles yields only a 1% change in the peak growth rate, and negligible changes in kmax and |η̂| (see figure 5.11c). We conclude that our choice of modified profiles is justified by these negligible changes in the properties of the unstable mode. The U∗ and N2∗ profiles are of similar form to the asymmetric stratified shear layer examined in the previous section with the exception that the two vorticity interfaces are of different strength, and in the location of the vertical boundaries. Once again we identify a σKH , σH and σself associated with the partial growth rates of the upper vorticity interface and the lower vorticity, density, and upper vorticity interfaces, respectively. Applying the diagnostic formulated in §5.3 leads to {σKH , σH , σself} = {0.011, 0.023,−0.010} s−1. We therefore conclude that the insta- bilities observed in figure 5.11 are primarily of the H-type. This appears to be in agreement with the wave-like features observed in the echosounding image (figure 5.11d) and the conclusions in chapter 4. 5.6 Discussion and conclusions Based on the wave interaction view of stratified shear instability (Holmboe, 1962; Baines and Mit- sudera, 1994; Caulfield, 1994), a diagnostic has been developed to interpret the wave interactions in stratified shear layers that lead to the growth of unstable modes. It can be viewed as an ex- tension of the large-α approximation used by Caulfield (1994) and Baines and Mitsudera (1994), Chapter 5. Unstable modes in asymmetric stratified shear layers 87 0 5 1 0 02468 1 0 1 2 z (m) -1 -0 .5 0 0 .5 0 0 .5 1 x ( m ) 0 2 0 4 0 6 0 (a ) (b ) (c ) (d ) Nρ − ρ 0 2 N *2 U U * U ’’ U *’ ’ |η| , |η ∗ | ^ ^ 1 0 0 1 0 0 Fi gu re 5. 11 : Pr ofi le s ta ke n fr om fie ld m ea su re m en ts in th e Fr as er R iv er es tu ar y (s ee ch ap te r 4 fo r de ta ils ). In (a -c ) th e or ig in al an d m od ifi ed pr ofi le s ar e de no te d by so lid an d da sh ed cu rv es , re sp ec tiv el y. T he de ns ity st ru ct ur e is pl ot te d in (a ) an d th e ve lo ci ty st ru ct ur e in (b ). T he no rm al iz ed am pl itu de of th e di sp la ce m en te ig en fu nc tio n, is sh ow n in (c ) w ith th e so lid ho ri zo nt al lin e in di ca tin g th e ve rt ic al lo ca tio n of th e pe ak di sp la ce m en t. A n ec ho so un di ng im ag e ta ke n at a tim e cl os e to th e pr ofi le s (a ,b ) is sh ow n in (d ), w he re in st ab ili ty ca n be se en ne ar th e pe ak in |η̂| . Chapter 5. Unstable modes in asymmetric stratified shear layers 88 to allow for a general classification of stratified shear instabilities in terms of the KH, H and T types. These types refer to vorticity-vorticity, vorticity-internal wave, and internal wave-internal wave interactions, respectively. The diagnostic has the advantage of quantifying wave growth in the entire αJ-plane, as well as being applicable to smooth profiles. The diagnostic is applied first to the symmetric and asymmetric stratified shear layers. An in- teresting result of the analysis was the presence of the KH growth mechanism in the propagating modes that are normally classified as H-type modes. The asymmetric stratified shear layer, whose region of instability consists entirely of propagating modes, was generally found to be composed of a mix of both KH- and H-type growth mechanisms. The diagnostic suggests that KH-type instabilities are found at larger values of J in the asymmetric case, and appears to be in quali- tative agreement with results from the numerical simulations presented here, as well as those of Carpenter et al. (2007), which cover a range of asymmetries. A drawback of the method is the large number of wave interactions that must be accounted for in profiles with multiple interfaces. In the field profiles examined from the Fraser River estuary, the number of observed interfaces would produce an unwieldy number of interactions that must be accounted for. This issue was dealt with by using a knowledge of the unstable mode to simplify the profiles without significantly affecting the stability properties. This reduced the analysis to that of the stratified shear layers examined earlier. However, this simplification may not always be possible, and there may arise cases in which the instability involves the interaction of numerous interfaces. Extending the use of the diagnostic to smooth profiles allows for the direct interpretation of wave interactions in geophysically relevant profiles measured in the field. These profiles invariably display some type of asymmetry, and the resulting modes of instability are likely to be of a ‘mixed’ or ‘hybrid’ type, i.e. they may involve two or more interaction types as is found in this study. This behaviour has been noted in numerous cases exhibiting asymmetry in a more general sense, for example, the non-Boussinesq flow of Umurhan and Heifetz (2007), and the spatially developing flows of Pawlak and Armi (1998) and Ortiz et al. (2002), in addition to the vertically offset profiles considered in the present study (see Lawrence et al., 1991). In these ‘mixed’ mode circumstances, the partial growth rate diagnostic can provide a useful tool in assessing the characteristics of the instability so that predictions can be made regarding the nonlinear dynamics. 89 Bibliography Alexakis, A. (2005). On Holmboe’s instability for smooth shear and density profiles. Phys. Fluids, 17:084103. Armi, L. and Farmer, D. (1988). The flow of Mediterranean water through the strait of Gibraltar. Prog. Oceanogr., 21:1–98. Baines, P. and Mitsudera, H. (1994). On the mechanism of shear flow instabilities. J. Fluid Mech., 276:327–342. Cairns, R. (1979). The role of negative energy waves in some instabilities of parallel flows. J. Fluid Mech., 92:1–14. Carpenter, J., Lawrence, G., and Smyth, W. (2007). Evolution and mixing of asymmetric Holm- boe instabilities. J. Fluid Mech., 582:103–132. Caulfield, C. (1994). Multiple linear instability of layered stratified shear flow. J. Fluid Mech., 258:255–285. Caulfield, C. and Peltier, W. (2000). The anatomy of the mixing transition in homogenous and stratified free shear layers. J. Fluid Mech., 413:1–47. Haigh, S. (1995). Non-symmetric Holmboe Waves. PhD thesis, University of British Columbia. Hazel, P. (1972). Numerical studies of the stability of inviscid stratified shear flows. J. Fluid Mech., 51:39–61. Helmholtz, H. (1868). On discontinuous movements of fluids. Philos. Mag., 36:337–346. Hogg, A. M. and Ivey, G. (2003). The Kelvin-Helmholtz to Holmboe instability transition in stratified exchange flows. J. Fluid Mech., 477:339–362. Holmboe, J. (1962). On the behavior of symmetric waves in stratified shear layers. Geofys. Publ., 24:67–112. Howard, L. and Maslowe, S. (1973). Stability of stratified shear flows. Boundary-Layer Met., 4:511–523. Kelvin, W. (1871). Hydrokinetic solutions and observations. Philos. Mag., 42:362–377. Koop, C. and Browand, F. (1979). Instability and turbulence in a stratified fluid with shear. J. Fluid Mech., 93:135–159. Bibliography 90 Lawrence, G., Browand, F., and Redekopp, L. (1991). The stability of a sheared density interface. Phys. Fluids, 3(10):2360–2370. Lawrence, G., Haigh, S., and Zhu, Z. (1998). In search of Holmboe’s instability. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 295–304. American Geophysical Union. Michalke, A. (1964). On the inviscid instability of hyperbolic-tangent velocity profile. J. Fluid Mech., 19:543–556. Ortiz, S., Chomaz, J., and Loiseleux, T. (2002). Spatial Holmboe instability. Phys. Fluids, 14:2585–2597. Pawlak, G. and Armi, L. (1998). Vortex dynamics in a spatially accelerating shear layer. J. Fluid Mech., 376:1–35. Rayleigh, J. (1880). On the stability, or instability, of certain fluid motions. Proc. Lond. Math. Soc., 12:57–70. Redekopp, L. (2001). Elements of instability theory for environmental flows. In Environmental Stratified Flows. Kluwer. Smyth, W., Carpenter, J., and Lawrence, G. (2007). Mixing in symmetric Holmboe waves. J. Phys. Oceanogr., 37:1566–1583. Smyth, W., Klaassen, G., and Peltier, W. (1988). Finite amplitude Holmboe waves. Geophys. Astrophys. Fluid Dyn., 43:181–222. Smyth, W. and Peltier, W. (1991). Instability and transition in finite-amplitude Kelvin-Helmholtz and Holmboe waves. J. Fluid Mech., 228:387–415. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Taylor, G. (1931). Effect of variation in density on the stability of superposed streams of fluid. Proc. R. Soc. Lond. A, 132:499–523. Thorpe, S. (1973). Experiments on instability and turbulence in a stratified shear flow. J. Fluid Mech., 61:731–751. Umurhan, O. and Heifetz, E. (2007). Holmboe modes revisited. Phys. Fluids, 19:064102. Wesson, M. and Gregg, M. (1994). Mixing at the Camarinal Sill in the strait of Gibraltar. J. Geophys. Res., 99:9847–9878. Bibliography 91 Winters, K., MacKinnon, J., and Mills, B. (2004). A spectral model for process studies of rotating, density-stratified flows. J. Atmos. Ocean. Tech., 21:69–94. Yoshida, S., Ohtani, M., Nishida, S., and Linden, P. (1998). Mixing processes in a highly strat- ified river. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 389–400. American Geophysical Union. 92 Chapter 6 Conclusions Stratified shear instabilities are an important process in the understanding of turbulent mixing in geophysical flows. This thesis departs from the great majority of previous research on this topic by focusing on the case of a relatively thin density interface embedded in a broader shear layer. Although the case in which the shear and stratification have the same length scale has been studied more extensively, thin sheared interfaces have been observed in a number of previous studies (e.g. Seim and Gregg, 1994; Yoshida et al., 1998; Sharples et al., 2003), all of which also exhibit relatively strong stratifications (J ≈ 0.5). These conditions have a richer dynamics as they are favourable to the generation of instability by the Holmboe mechanism, as well as the Kelvin-Helmholtz. In these stratified shear flows unstable waves are generated on the density interface – and it is these waves that are the topic of this thesis. Each of the four main chapters of the thesis are summarized below, and in each case the original contributions to the study of unstable waves on a sheared density interface are described. In closing, some directions for future research are outlined. 6.1 Summary In the first manuscript chaper, a study of the wave fields generated by the Holmboe instability was undertaken. Two different flows were considered: one in which the mean flow exhibited gradual temporal variations, corresponding to the direct numerical simulations (DNS), and laboratory ex- periments where the flow displayed a gradual spatial variation. This chapter extends the current knowledge of shear instabilities by describing the changes to the wave field that are caused by these gradual variations in the mean flow. In the case of the laboratory experiments, the spatially accelerating mean flow is responsible for ‘stretching’ the waves to longer wavelengths. The appli- cation of gradually varying wave theory has been found to provide an excellent description of this stretching phenomena – despite the constant formation of new waves by the instability process. The effect of the stretching on wave amplitude can also be estimated using a simple application of the conservation of wave action, showing that significant reductions in amplitude can be expected. These experimental results are compared and contrasted to results from multiple wavelength DNS, which have not previously been performed for the Holmboe instability. This led to the identifica- tion of two important processes for the ‘coarsening’ of the wave field, namely, vortex pairing and ejections. Although the vortex pairing process has been described extensively for KH instabili- Chapter 6. Conclusions 93 ties, it has not been observed in symmetric Holmboe instabilities. On the other hand, the ejection process has frequently been observed by numerous investigators (e.g. Thorpe, 1968; Koop and Browand, 1979). However, the importance of ejections as a wave breaking mechanism, and its role in the coarsening of the wave field has been described for the first time. Given the importance of the ejection process, the next chapter focuses on describing the basic mechanism governing its occurrence. This is due to the formation of a stratified vortex couple that consumes the impulse of the shear layer in transporting the ejected fluid against buoyancy forces. The proposed mechanism is able to qualitatively explain observations of ejections in the DNS such as the dependence on interactions between the upper and lower wave modes, and the view of ejections as wave breaking events. This qualitative description of the wave breaking mechanism is a first step in quantifying important wave properties such as maximum wave amplitudes. In chapter 4, shear instabilities are identified in the important geophysical scale flow of the highly stratified Fraser River estuary. Observations of instabilities, made through the use of an echosounder, were verified directly by a stability analysis of the measured profiles. The close agreement between the observations and the linear predictions demonstrates that even in a highly energetic, evolving environment, linear theory can be used to predict features of the instabilities with reasonable accuracy. Previous studies in the Fraser River estuary suggest that it is the gener- ation of shear instabilities that is responsible for the majority of the mixing in the estuary (Geyer and Smith, 1987; MacDonald and Horner-Devine, 2008). These studies have generally focused on the high freshwater discharge conditions, where the salinity intrusion more closely resembles a classic salt wedge. However, the observations presented in chapter 4 show that the structure of the intrusion becomes much more complex, displaying multiple layers, with bottom boundary generated mixing playing a more important role. In addition, the observations show that instabil- ity is generally associated with an asymmetry in the locations of strong shear and strong density gradients. This asymmetry has important consequences for the type of instability that forms, and is the subject of chapter 5. When an asymmetry is present in the profiles of sheared density interfaces it introduces am- biguity in the type of unstable modes that may arise. Focusing on the case of a single density interface allows both Kelvin-Helmholtz and Holmboe-type modes. Since the linear growth mech- anism of these unstable modes can be interpreted as the interaction of two different wave pairs, a diagnostic is developed that quantifies these wave interactions. This diagnostic extends previ- ous methods of mode identification so that it can now be applied to smooth profiles, and for all wavenumbers. The transition from Kelvin-Helmholtz to Holmboe modes in asymmetric stratified shear layers is examined, as well as an instability observed in the Fraser River estuary. Asym- metry is found to lead to ‘mixed’ or ‘hybrid’ modes that are composed of both Kelvin-Helmholtz and Holmboe growth mechanisms, thus illustrating that instability may involve the interaction between numerous waves. Chapter 6. Conclusions 94 6.2 Future work This thesis has identified a number of possible avenues for future research. Since linear theory has been found successful in describing the basic features of the instabilities, it is an appropriate starting point in the study of unstable waves on a sheared density interface. It was found in chapter 2 that predicting the nonlinear evolution of the wave field is intimately linked to both the instability process and the mean flow. Depending on whether the mean flow is varying in time or space, different processes are expected to act on the wave field. However, it remains unclear what overarching principles are responsible for the changes in some of the basic wave parameters such as wavenumber and amplitude. For example, the temporally varying conditions of the simulations show a ‘wave coarsening’ effect which was found to be related to both a vortex pairing mechanism and the ejection process. But what controls the overall rate of coarsening, and how are the instabilities responding to changes in the mean flow? This relationship between the mean flow and the instabilities represents a challenging problem that is perhaps only present in stratified shear layers that exhibit the Holmboe mode of instability. In this case the instability leads to persistent finite amplitude waves that do not immediately break down into turbulent patches, as found in Kelvin-Helmholtz instabilities. However, turbulence is expected to also play a role in the dissipation of the Holmboe wave field, particularly at larger Reynold’s numbers than that employed in this thesis. The greatest area requiring further study is in the application of our current understanding of shear instability to assess its role in geophysically relevant flows. While chapter 4 showed that linear theory could be applied to geophysical observations of shear instabilities with some confidence, it also raised a number of questions regarding the interpretation of the observations and their relevance to the mixing of mass and momentum. Observing instabilities from a moving frame of reference will produce a distorted view due to the Doppler shifting of the wave features. In addition, it is not possible to identify instabilities (particularly with echosounders) after they have broken down and become sufficiently turbulent. Direct numerical simulations of stratified shear flow instabilities have shown that this turbulent phase of the instability life-cycle can persist for a considerable duration (Caulfield and Peltier, 2000; Smyth and Winters, 2003). The turbulent phase is more likely to be the natural state of the high-Reynold’s number instabilities observed in geophysical flows, and makes for difficulties in quantifying the importance of the instability process in the overall dynamics of the flow. Of particular importance is the quantification of mixing due to shear instabilities, with the ultimate goal being a parameterization of shear driven mixing processes. Since mixing character- istics have been found to depend on the type of instability that occurs (Smyth and Winters, 2003; Smyth et al., 2007; Carpenter et al., 2007), this task must account for the instability mechanisms in each particular flow. A method of classifying instabilities by their linear growth mechanisms was described in chapter 5 that aimed at predicting the nonlinear characteristics. This is perhaps Chapter 6. Conclusions 95 a first step towards assessing the type of mixing that may occur at a sheared density interface, and whether the instabilities consist of Kelvin-Helmholtz-like billows or Holmboe-like waves. Future research in this area would link the instability mechanisms of the stratified shear layer with the nonlinear mixing characteristics and wave properties of the resulting turbulent flow. 96 Bibliography Carpenter, J., Lawrence, G., and Smyth, W. (2007). Evolution and mixing of asymmetric Holm- boe instabilities. J. Fluid Mech., 582:103–132. Caulfield, C. and Peltier, W. (2000). The anatomy of the mixing transition in homogenous and stratified free shear layers. J. Fluid Mech., 413:1–47. Geyer, W. and Smith, J. (1987). Shear instability in a highly stratified estuary. J. Phys. Oceanogr., 17:1668–1679. Koop, C. and Browand, F. (1979). Instability and turbulence in a stratified fluid with shear. J. Fluid Mech., 93:135–159. MacDonald, D. and Horner-Devine, A. (2008). Temporal and spatial variability of vertical salt flux in a highly stratified estuary. J. Geophys. Res., 113. Seim, H. and Gregg, M. (1994). Detailed observations of naturally occurring shear instability. J. Geophys. Res., 99:10,049–10,073. Sharples, J., Coates, M., and Sherwood, J. (2003). Quantifying turbulent mixing and oxygen fluxes in a Mediterranean-type, microtidal estuary. Ocean Dyn., 53:126–136. Smyth, W., Carpenter, J., and Lawrence, G. (2007). Mixing in symmetric Holmboe waves. J. Phys. Oceanogr., 37:1566–1583. Smyth, W. and Winters, K. (2003). Turbulence and mixing in Holmboe waves. J. Phys. Oceanogr., 33:694–711. Thorpe, S. (1968). A method of producing a shear flow in a stratified fluid. J. Fluid Mech., 32:693–704. Yoshida, S., Ohtani, M., Nishida, S., and Linden, P. (1998). Mixing processes in a highly strat- ified river. In Physical Processes in Lakes and Oceans, volume 54 of Coastal and Estuarine Studies, pages 389–400. American Geophysical Union.