CRWR Online Report 04-09 Hydrostatic and Non-hydrostatic Internal Wave Models by Bridget M. Wadzuk, B.S., M.S Graduate Research Assistant and Ben R. Hodges, Ph.D. Principal Investigator The University of Texas at Austin December 4, 2004 Final Report to the Office of Naval Research under Contract No. N00014-01-1-0574 This document is available online via World Wide Web at http://www.crwr.utexas.edu/online.html Copyright by Bridget Marie Wadzuk 2004 Abstract Numerical models have become an indispensable tool for ocean and inland flow modeling. Such models typically use the hydrostatic approximation based on the argument that their horizontal length scales are longer than the vertical length scales. There are a wide variety of physical processes in oceans and inland water systems, and many of these processes are adequately modeled with the hydrostatic approximation. However, internal waves contribute to the physics that influence mixing in a density stratified system and have been previously shown to be non- hydrostatic. The neglect of non-hydrostatic pressure in a hydrostatic model is problematic since non-hydrostatic pressure plays a significant role in internal wave evolution balancing nonlinear wave steepening. Where non-hydrostatic pressure is neglected in a model, the governing equations are missing a piece of the physics that control the internal wave evolution, so it should not be surprising that the evolution may be poorly predicted. Despite the knowledge that the non-hydrostatic pressure is necessary for correctly modeling the physics of a steepening internal wave, the high computational cost of solving the non-hydrostatic pressure has limited its use in large-scale systems. Furthermore, the errors associated with hydrostatic modeling of internal waves have not been quantified. This research quantifies the differences between hydrostatic and non-hydrostatic simulations of internal wave evolution and develops a method to a priori determine regions with non-hydrostatic behavior. In quantifying the errors and differences between the two models this research provides the characteristics of model error with grid refinement. Additionally, it is shown that hydrostatic models may develop high wavenumber ?soliton-like? features that are purely a construct of model error, but may seem to mimic physical behaviors of the non-hydrostatic system. Finally, it is shown that regions of significant non-hydrostatic pressure gradients can be identified from a hydrostatic model. This latter finding is a building block towards coupling local non- hydrostatic solutions with global hydrostatic solutions for more efficient computational methods. The work presented here provides the foundations for future non-hydrostatic model development and application. i Chapter 1. Introduction Abstract ............................................................ i Chapter 1. Introduction............................1 1.1 Motivation................................................1 1.2 Objectives ................................................3 1.3 Background..............................................3 1.3.1 Internal Waves.................................. 3 1.3.2 Modeling........................................... 5 1.4 Approach..................................................6 Chapter 2. Methods...................................7 2.1 Numerical Model .....................................7 2.1.1 Governing Equations........................ 7 2.1.2 Solution Methods .............................. 8 2.1.3 Iterative Methods ? Pressure Poisson Equation 10 2.1.4 Boundary Conditions...................... 10 2.2 Analysis Methods ..................................10 2.2.1 Background..................................... 11 2.2.2 Numerical Diffusion ....................... 13 2.2.3 Numerical Dissipation.................... 16 2.2.4 Spectral Analysis ............................ 17 Chapter 3. Verification and validation of the non-hydrostatic solver ............................19 3.1 Verification by the Method of Manufactured Solutions...............................19 3.1.1 Setup ............................................... 19 3.1.2 Error Analysis................................. 19 3.2 Validation and Convergence for an Unsteady Internal Wave...............................22 3.2.1 Setup ............................................... 22 3.2.2 Results............................................. 23 3.3 Validation Against Regime Theory and Laboratory Experiment................................25 3.3.1 Regime Theory................................ 25 3.3.2 Horn?s Laboratory Experiment and Model Simulations.......................................... 26 3.3.3 Results............................................. 28 3.4 Summary................................................32 Chapter 4. Comparison of Hydrostatic and Nonhydrostatic models ..........................33 4.1 Laboratory Scale Comparison ...............33 4.1.1 Setup ............................................... 33 4.1.2 Scenario Simulations ...................... 34 4.1.3 Changing the Grid Resolution ........ 39 4.2 Lake-Scale Comparison.........................47 4.2.1 Setup ............................................... 47 4.2.2 Discussion....................................... 52 4.3 Summary................................................53 Chapter 5. Isolation of non-hydrostatic effects .............. ...............................................59 5.1 Numerical Theory..................................59 5.2 Application of the Non-hydrostatic Pressure Isolation Method...........................62 5.2.1 Setup................................................62 5.2.2 Results.............................................62 5.3 Summary ...............................................64 Chapter 6. Conclusions and Recommendations .........................................71 6.1 Summary Discussion.............................71 6.2 Conclusions ...........................................72 6.3 Recommendations .................................73 Acknowledgements .......................................73 Bibliography ..................................................75 iii iv Chapter 1. Introduction Chapter 1. Introduction 1.1 Motivation Biogeochemical processes in oceans, lakes and estuaries evolve against the background of vertical density stratification. Internal waves are one mechanism which carry energy and momentum through a basin and contribute to mixing events in littoral regions that prompt ecosystem changes. Internal waves are motions that occur beneath the free-surface of a density-stratified waterbody. As an external force (i.e. wind or river inflow) moves the density layers from their equilibrium position, an internal wave is initiated to restore the system to equilibrium (Figure 1.1). Internal waves propagate through the basin and interact with the basin boundaries; this internal wave ? slope interaction is an important source of energy that transports nutrients, biota and contaminants through the water column (Imberger and Ivey, 1993; Javam, et al., 1999). Hydrodynamic processes (e.g. internal waves, tides and eddies) are transport mechanisms in coastal oceans and lakes and their scope and magnitude depends on the processes? speed and length scales. Table 1.1 shows a large range of scales for different hydrodynamic processes and internal waves are considered a middle-scale process. Quantitative assessment of nutrient and constituent transport in a stratified basin requires modeling of internal wave evolution (Imberger, 1994). To have an unambiguous water quality model, the complementary hydrodynamic model, which provides the flow field, should simulate all the significant physics (Gross, et al., 1999). Hydrodynamic models often neglect physics that are irrelevant to the focus of an investigation or can be a priori scaled as small compared to dominant (Marshall, et al., 1997). Large-scale ocean and littoral modeling have traditionally made the hydrostatic approximation, which neglects non- hydrostatic pressure and subsequently vertical momentum (e.g. POM: Blumberg and Mellor, 1987; EFDC: Hamrick, 1992; ROMS: Haidvogel, et al., 2000; Ezer, et al., 2002). The hydrostatic approximation fails at open boundaries (Mahadevan, et al., 1996a) and at steep slopes with strong vertical velocities (Weilbeer and Jankowski, 2000; Horn, et al., 2001; Horn, et al., 2002). Linear waves that are damped by viscous effects may be considered to b)a) c) d) Figure 1.1: Internal wave evolution in a two-layer system. a) System at rest, where the equilibrium position for the free-surface (thick dashed blue line) and pycnocline (thick dashed red line) is flat. b) Wind is applied over the free surface (large red arrow) which moves the free-surface (thin blue line) and pycnocline (thin red line) from their respective equilibrium positions. c) Continued application of wind over the free-surface, so the equilibrium positions for the free-surface and pycnocline are tilted. The free-surface and pycnocline oscillate around the equilibrium position. d) The wind force stops and the free-surface and pycnocline equilibrium position is flat. The free-surface and pycnocline oscillate around the equilibrium position, setting up a basin-scale seiche. 1 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 2 Process/Activity Speed (m/s) Length Oceanic Turbulence 0.1 - 100 1 cm - 100 m Internal Waves 1 - 3 100 m - 1 km Diurnal Tides 0.1 - 10 1 km - 1000 km Mesoscale Eddy 1 10 km - 1000 km Rossby Waves 0.1 - 9 1000 km - 10000 km Surface Gravity Waves 1 - 20 1 m - 100 m Sound waves 1400 200 km behave hydrostatically, while internal waves that steepen nonlinearly are inherently non-hydrostatic (Long, 1972). It is suggested that non-hydrostatic pressure is essential for proper modeling of internal wave development and propagation (Laval, et al., 2003). When numerical models neglect non- hydrostatic pressure, they are neglecting a physical attribute of internal wave evolution. This neglect has repercussions on vertical transport at the topographic boundaries of a basin. Internal waves that nonlinearly steepen and degenerate into solitons can propagate to a basin?s boundary where the wave may shoal and mix, energizing a turbulent benthic boundary layer (Boegman, et al., 2003); the interaction of internal waves and boundary layer mixing is a significant mechanism for nutrient transport and other biogeochemical processes (De Silva, et al., 1997; Nishri, et al., 2000). Modeling the boundary layer is not of interest to the present research, but properly modeling the internal wave evolution which contributes to boundary layer is. During the past decade there has been interest simulating internal waves in numerical models, marked by the development of several non- hydrostatic models (e.g. Mahadevan, et al., 1996a,b; Marshall, et al., 1997; Stansby and Zhou, 1998; Casulli, 1999; Weilbeer and Jankowski, 2000; Fringer and Street, 2003; Daily and Imberger, 2003; Yuan and Wu, 2004). However, there has been little work done on quantifying the differences between hydrostatic and non-hydrostatic models or the scales for which non-hydrostatic pressure impacts internal wave evolution. If the proper scales (i.e. sufficient grid resolutions) are not used in a non-hydrostatic model, then the results of the non-hydrostatic model are no different than a hydrostatic model?s solution, as described in ?4.2 Presently, computational power is insufficient for application of non-hydrostatic models to large domains (i.e. coastal oceans) with sufficient grid resolution. By example, one computational ?work unit? may be considered the processing power required for one operation across the entire domain. Thus, a work unit scales on the number of cells within a domain (e.g. a grid of N ? M size: 1 work unit = N ? M). A hydrostatic model?s basic operations are solving for the three velocity components, salinity transport, temperature transport and the free-surface; the computational effort for each of these operations scales on one work unit, thus totaling six work units for each timestep. A non-hydrostatic model solves for these operations twice (once before the non-hydrostatic solver and once after) and a pressure solution. The work associated with the pressure solution is N ? M ? the number of smoothing iterations per timestep; the number of iterations per timestep for the non- hydrostatic pressure solver may be O(10 2 ) (?3.2 ff.). If the non-hydrostatic model is applied with the same spatial grid and timestep as the hydrostatic model, the non-hydrostatic model requires an additional 10 6 work units, that is, about two orders of magnitude more computational effort. However, as discussed in ?4, the non-hydrostatic model also requires a finer spatial grid and timestep than a hydrostatic model. A typical hydrostatic coastal ocean model applies a spatial grid with an aspect ratio (?z/?x) of O(10 -3 ) (Marshall, et al., 1997; Tartinville, et al., 1998), while ?4 shows that a non-hydrostatic model requires a grid with an aspect ratio of at least O(10 -2 ). This refined spatial grid increases the required work units by another order of magnitude, to ~10 3 . Furthermore, ?4 discusses the need for a non- hydrostatic model to be applied to a smaller timestep, which increases the number of work units by another order of magnitude. Thus, O(10 4 ) work units are needed for the non-hydrostatic model, while the hydrostatic computational effort scales on O(10) work units. According to Moore?s law (Moore, 1965), computational ability doubles about every two years. Assuming that we are reaching the maximum capacity of present computers with the hydrostatic model, the computational power needed to run a coastal ocean non-hydrostatic model with a sufficiently fine spatial grid and timestep will not be available for another 26 years! It is computational power that limits our ability to use non-hydrostatic models. This motivates the present research to quantify the spatial grid and timestep necessary to apply a non-hydrostatic model and develop the theoretical foundations for future methods that may circumvent the constraint of computational power. Table 1.1: Scales of oceanic processes (Gill, 1982; Cushman-Roisin, 1994; Kantha & Clayson, 2000a; Miropol?sky, 2001) Chapter 1. Introduction 1.2 Objectives The present research has two main objectives to address the motivational issues within internal wave modeling: ? quantify the differences between hydrostatic and non-hydrostatic simulations of internal wave evolution, and ? develop a method to a priori determine regions with non-hydrostatic behavior. The first objective is achieved by: 1) quantifying accumulation of error within both hydrostatic and non-hydrostatic models of internal waves; and 2) comparing the space-time evolution of internal wave characteristics in both models. Comparing the models? internal wave evolution and error accumulation provides a means of defining the conditions for which a hydrostatic model is acceptable and those conditions requiring a non- hydrostatic model to capture an internal wave?s propagation. It is shown that a model?s error accumulation and representation of internal wave evolution are dependent on spatial and temporal resolution. Delineation of the spatial and temporal resolution necessary to adequately resolve non- hydrostatic processes is included within the first objective. The second objective is achieved by using a hydrostatic model to estimate non-hydrostatic pressure effects. Using this estimation, local regions where the non-hydrostatic pressure significantly contributes to internal wave evolution are identified. Through reaching this objective, the present work provides the theoretical basis for future model development for a non-hydrostatic solution in isolated, ?non-hydrostatic? areas, concurrent with a hydrostatic solution elsewhere. 3 1.3 Background 1.3.1 Internal Waves Internal waves occur in stably stratified fluids when a water parcel is displaced by some external force (wind, inflow, etc.) and is restored by buoyancy forces; the restoration motion may overshoot the equilibrium position and set up an oscillation thereby forming an internal wave (Turner, 1973; Kantha and Clayson, 2000a). The buoyancy frequency (N) is the upper limit of wave frequencies (?) that can propagate through a system; thus, the only internal waves of interest for the present research are those where N?< , as waves where remain local and do not significantly contribute to the system?s overall transport and energy (Turner, 1973; Lighthill, 1978; Javam, et al., 1999). The lower limit of internal wave frequencies is the Coriolis or inertial wave frequency, f. The inertial frequency is defined as f = 2?sin?, where ? is the angular velocity of the earth?s rotation and ? is the latitude. At f, wave motions become inertial oscillations, where fluid parcels have horizontally circular trajectories (Cushman-Roisin, 1994). While Coriolis forces are known to affect internal wave evolution, their impact is principally three-dimensional. To limit the number of processes under consideration, this project is restricted to vertical stratification in two- dimensional box models where Coriolis effects cannot be represented. Thus, this project is focused primarily on internal waves dominated by gravity, which is typically the case for waves breaking on sloping boundaries (De Silva, et al., 1997). N?> Field studies by Garrett and Munk (1979) indicated that most internal waves follow a consistent energy-wave spectrum (Figure 1.2). There is a full range of wave frequencies and associated energy, where waves of low-frequency are associated with high energy, while those of high- frequency have low energy. Internal waves also evolve over a wide spectrum of spatial scales. An initial low-frequency, long wavelength internal wave may degenerate into a train of high-frequency, short Figure 1.2: Energy-Frequency Spectrum. f is the Coriolis frequency, N is the Brunt-V?is?l? frequency and ? is the frequency, cph (cycles per hour) is a frequency scale and m 2 /cph is an energy scale. (Garrett and Munk, 1979) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models wavelength waves (Kao, et al., 1985; Horn, et al., 2001; Boegman, et al., 2003), as seen in Figure 1.3. This degeneration is connected to the nonlinearity of internal waves. Several field studies (Farmer, 1978; Osborne and Burch, 1980; Apel, 1981; Wiegand and Carmack, 1986) and model results (Segur and Hammack, 1982; Hutter, et al., 1998; Horn, et al., 2002) have verified that an internal wave will behave nonlinearly if the initial wave is sufficiently steep. Steepness is defined as the ratio of a wave?s amplitude to its wavelength. A ?sufficiently? steep wave is one whose timescale for steepening is shorter than its timescale for damping (Farmer, 1978; Horn, et al., 2001). The relationship between the total water depth and thickness of the upper layer can also contribute to a wave?s nonlinearity. The following description of wave evolution is adapted from Kinsman (1965), Whitham (1974) and Lighthill (1978). A wave may be viewed as a recognizable feature of a disturbance that moves with a finite velocity. An infinite number of sinusoidal waves, each with a discrete frequency and discrete wavelength, may be superposed to create the characteristic material surface recognizable as a wave (Figure 1.4 A). Wave speed (c): c k ? = (1.1) where ? is the wave frequency and k is the total wavenumber (i.e. 22 2 kjlm=++; j, l and m are the wavenumbers in the x-, y- and z-directions, respectively) is dependent on the dispersion relationship for an internal wave: ( ) () 22 22 22 2 jl N jlm + ?= ++ (1.2) where the buoyancy frequency, N 2 is: 2 o gd N dz ? = ? (1.3) It follows that wave speed is dependent upon the density gradient of fluid. Waves on different density structures propagate at different speeds (i.e. the crest of the wave, with one density gradient, is moving faster than the trough of the wave, with another density gradient, as seen in Figure 1.4 B&C). This density gradient dependency produces a nonlinearly steepening wave. Wave speed is also dependent on wavelength (the inverse of the wavenumber). As the different sinusoidal waves comprising the material wave move at different speeds, the material wave spreads out over time, or disperses, into a series of waves (Figure 1.4 D). The wave speed dependence on the dispersion relationship shows there is a relationship between Figure 1.3: Laboratory experiment, conducted by Horn, et al. (2001), of a steepening, initial basin-scale wave that degenerates into a train of solitons. The first panel is the initial wave setup. The second panel shows the wave front steepening. The third panel shows the beginning of smaller waves developing behind the lead wave. The last two panels are the train of solitons traveling through the domain. The basin is 6m long and 0.29m high. 4 Chapter 1. Introduction nonlinear acceleration and a dispersive force. The dispersive force that opposes nonlinear wave deformation is the non-hydrostatic pressure gradient. The dispersion of a wave leads to an energy cascade from an initial, low-frequency wave to a resulting, high frequency wave series that is more susceptible to wave breaking (van Haren, 2004). When non- hydrostatic pressure is not modeled while nonlinearities are modeled (i.e. typical of a hydrostatic model), there are two possible results. The first possibility is a wave steepens unabated and creates artificial mixing events through density overturns. The second possibility is a build up of numerical error that balances nonlinear steepening. A The two dominate forms of numerical error are numerical diffusion of mass (from herein called numerical diffusion) and numerical diffusion of momentum. The latter results in dissipation of energy and is henceforth called numerical dissipation. Numerical diffusion artificially weakens sharp gradients across the wave front. Thus, as a wave steepens, numerical diffusion increases and smoothes the wave front. This smoothing impedes wave steepening and reduces the available potential energy in wave, but increases the background potential energy of the system (Figure 1.5 A). The increase in background potential energy is seen in the thickening of the layer that contains the density gradient between the constant density upper and lower layer, which is called the pycnocline (see ?2.2.1 for a description of available and background potential energy). Numerical dissipation has a similar effect in damping wave steepening. As a wave steepens, the velocity shear increases across the wave front, thereby increasing the artificial spreading of momentum through model error. The artificial spreading of momentum reduces the kinetic energy of the wave and inhibits steepening (Figure 1.5 B). B C D Figure 1.4: Cartoon of wave steepening and dispersion. A) An initial wave with a unique wave speed and amplitude. B&C) Wave steepening where the crest is moving faster than trough. D) Wave degeneration into a train of waves. 1.3.2 Modeling Several hydrodynamic models for large-scale applications use the hydrostatic approximation (e.g. POM: Blumberg and Mellor, 1987; EFDC: Hamrick, 1992; ROMS: Haidvogel, et al., 2000). The hydrostatic approximation neglects non-hydrostatic pressure and vertical acceleration, but models B A 5 Figure 1.5: Cartoon of wave steepening. A) Nonlinear steepening wave with numerical diffusion. As the wave steepens, diffusion acts across the wave front and reduces the density gradient, available potential energy and steepness. B) Nonlinear wave steepening with numerical dissipation. The initial wave is primarily moving in a horizontal direction. As the wave steepens, the vertical velocities in front and behind the wave increase, increasing the shear. The increased diffusion of momentum (dissipation) reduces the wave?s kinetic energy and limits the wave steepness. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 6 nonlinear horizontal momentum. The hydrostatic approximation is generally applicable to processes with small aspect ratios (depth: length ? O[10 -3 ]); for example, most large-scale ocean and atmospheric problems. However, the hydrostatic approximation breaks down for mesoscale systems [aspect ratio ~ O(10 -1 ? 10 -2 )] and small-scale events, such as the steepening and breaking of nonlinear internal waves (Kantha and Clayson, 2000a). Hydrostatic models often predict excessive steepening (Wadzuk and Hodges, 2003), which induces artificial breaking and diffusion, thereby developing a weakened pycnocline (Laval, et al., 2003). Hydrostatic models also show increased dissipation of internal wave energy which artificially reduces internal wave amplitude. Internal wave dispersion is not physically modeled in hydrostatic models (as it requires non-hydrostatic pressure), but may still appear as a numerical phenomenon in a simulated evolution (Hodges and Delavan, 2004). The aforementioned errors, which are all numerical, result in an incorrect evolution of an internal wave and its energy. Thus, the need for non-hydrostatic models is evident (Horn, et al., 2001; Laval, et al., 2003). Several quasi-hydrostatic and non-hydrostatic models for large-scale systems have been developed (e.g. Mahadevan, et al., 1996a,b; Casulli and Stelling, 1998; Casulli, 1999; Fringer and Street, 2003). Despite this work, there has been little research on quantifying the timestep and spatial grid scales necessary to model non-hydrostatic behavior. The present research shows that a fine spatial and temporal resolution is needed to capture non- hydrostatic behavior (?4 ff.), thus significantly increasing computational time over hydrostatic models. The aforementioned models applied the non-hydrostatic model over the entire domain, but Stansby and Zhou (1998) recognized that non- hydrostatic effects are local and suggested application of the non-hydrostatic solution only in these limited areas. However, their work provides neither a method to locally solve the pressure, nor a method to locate regions where the non-hydrostatic pressure is important. The present research has developed a method to numerically identify regions with significant non-hydrostatic effect. 1.4 Approach Achieving the objectives of this study required both hydrostatic and non-hydrostatic models. The University of Western Australia Centre for Water Research?s Estuary and Lake Computer Model (CWR-ELCOM: Hodges, 2000) was used as the hydrostatic model. As part of this project, a non- hydrostatic pressure solver was developed as a modification to CWR-ELCOM. CWR-ELCOM is a three-dimensional hydrodynamic model that applies the incompressible, Boussinesq and hydrostatic approximations to solve the Euler equations on an Arakawa C staggered grid with a semi-implicit method and a moving free-surface. The model temporal accuracy is first-order for the flow field, third-order ULTIMATE-QUICKEST (Leonard, 1991) for scalar transport and first-order for the free surface. Surface thermodynamics are solved by bulk transfer models that are not used or investigated in this work. The Euler-Lagrange method for momentum is third-order spatially accurate, while spatial gradients in other terms are second-order. Details of the non-hydrostatic solver developed for the present research are described in ?2 ff. The existing and modified models were used to perform the analysis fulfilling the two research objectives. The first chapter of this document provides the reader with the motivation, objectives and background of this project. The second chapter provides the numerical methods applied in the hydrostatic model and the design of the non- hydrostatic solver, as well as the analytical methods used to quantify the performance of hydrostatic and non-hydrostatic models. The third chapter demonstrates the validation and verification of the non-hydrostatic solver. Chapter four uses the tools developed in chapter two to examine the first objective (comparing the differences between hydrostatic and non-hydrostatic models), including an evaluation of model skill on predicting internal wave evolution and quantification of numerical error. Chapter five develops and demonstrates a method to isolate non-hydrostatic effects. Chapter six concludes this work and makes suggestions for future work in this area of hydrodynamic modeling. Chapter 2. Methods Chapter 2. Methods Several different numerical and analytical methods are used in the present work to develop the models used in analysis and to provide a means to quantify and compare model results. Governing equations (?2.1.1) are used to model the physics of internal wave evolution and are applied to the numerical model with a pressure Poisson solution implemented by the fractional-step method described in ?2.1.2 and ?2.1.3. Along with the governing equations and numerical methods, a set of approximations and conditions are made to simplify the problem and provide boundaries for the solution space. In assessment of model skill, both qualitative and quantitative techniques are used. Section 2.2 provides qualitative methods for assessing model skill by examining the pycnocline displacement and quantitative methods for evaluating model skill through computing numerical dissipation and diffusion. 2.1 Numerical Model This sections discusses a description of the governing equations of CWR-ELCOM and the non- hydrostatic solver and a justification for the inviscid approximation (?2.1.1). The methods used to develop the non-hydrostatics solver (?2.1.2 and ?2.1.3) and the applicable boundary conditions (?2.1.4) are also described. 2.1.1 Governing Equations The governing equations for the present research are the Euler equations, the free-surface equation with the kinematic boundary condition and the linear equation of state for density. In Cartesian tensor form, these are: Incompressible momentum:  () i i o z' nh oi DU 1 g Dt x x P1 x ? ? ?? ? ?? ? ? =?? + ? ??? ? ? ? ?? ? dz ? ? ? (2.1) Incompressible continuity: i i U 0 x ? = ? (2.2) Free-surface ? Kinematic boundary condition integrated over depth with continuity applied: b Udz tx ? ? ? ?? ? =? ?? ? (2.3) Density: ( ) o 1S?=? +? (2.4) In Equations (2.1) - (2.4), i = 1, 2, 3; ? = 1, 2; U i is the velocity; ? is the free-surface elevation; z? is the vertical position in the water column; P nh is the non- hydrostatic pressure; ? is the density; ? o is the reference density; g is gravity. The haline expansion coefficient, ?, is constant for a linear equation of state and S is salinity. The density can be decomposed as:  o ?=? +? (2.5) where  ? is a small departure from the reference density. When the hydrostatic approximation is used, as in CWR-ELCOM, the non-hydrostatic pressure gradient is neglected in Equation (2.1). Within the hydrostatic limit, vertical motions are assumed small (horizontal velocity ~ 10 -1 ms -1 , vertical velocity ~ 10 -4 to 10 -5 ms -1 ); therefore the vertical momentum equation is neglected (Cushman- Roisin, 1994). However, Long (1972) showed that the inclusion of vertical acceleration and non- hydrostatic pressure is necessary to properly model internal wave evolution. Thus, CWR-ELCOM has been modified as described in this chapter to include a non-hydrostatic pressure component based on Equation (2.1). The present research?s focus is on quantifying hydrostatic and non-hydrostatic models? ability to represent internal wave evolution and error accumulation. All analysis for the present research neglects viscous effects and diffusion terms in the momentum and transport equations for two reasons: 1) the effects of viscous damping and viscous effects in wave breaking are at the two extremes of wave phenomena, which are of lesser interest than the evolution of solitons; and 2) by neglecting viscosity and diffusion the dissipation of energy and diffusion of mass can be used to quantify numerical error. Viscous damping acts primarily through the wall boundary layers, which are not modeled in the present research, and the shear between density 7 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models layers. The viscous effect between density layers is immaterial in the present work as the timescale for damping is longer than the timescale for steepening, which is a phenomenon of interest for an evolving nonlinear wave. The modulus of decay provides the timescale for significance of viscous effects. This viscous timescale is the time it takes for a wave?s amplitude to decrease to e -1 of its initial amplitude (Lamb, 1932): ? (m) ? (s) ? (d) 10000 1.267E+12 1.466E+07 100 1.267E+08 1.466E+03 1 1.267E+04 1.466E-01 0.01 1.267E+00 1.466E-05 1.0E-05 1.0E-03 1.0E-01 1.0E+01 1.0E+03 1.0E+05 1.0E+07 1.0E+09 0.01 0.1 1 10 100 1000 10000 Wavelength (m) M odul us of D ecay ( d ) 2 2 8 ? ?= ?? (2.6) where ? is the modulus of decay, ? is wavelength and ? is viscosity. Figure 2.1 shows the modulus of decay for several different wavelengths. This viscous timescale is small for very short waves (e.g. capillary waves, Kinsman, 1965), so viscosity is important in their evolution and attenuation. Conversely, the viscous timescale for gravity waves (e.g. ? > O[1m]) is long compared to the timescales of forcing and wave evolution (see discussion of wave steepening timescale, ?3.3.3 ff.), so viscosity is insignificant for the evolution of longer waves. Viscous and diffusive effects in wave breaking and mixing are not of interest in the present work as large-scale coastal oceans and lake models do not have sufficient grid resolution to capture the details of these processes. The present research is focused on wave behavior prior to breaking, when the wave can be considered a smooth material surface that effectively suppresses turbulent diapycnal diffusion of mass and momentum. Figure 2.1: Modulus of decay for as a function of wavelength, determined from Equation (2.6). Viscosity is 10 -6 m 2 /s. 2.1.2 Solution Methods The fractional step method (Kim and Moin, 1985) was used to incorporate the non-hydrostatic solver into the existing hydrostatic CWR-ELCOM. The fractional step method is mass conservative and has second-order temporal accuracy for 3D incompressible flows with fine spatial discretization and coarse temporal discretization (Casulli and Stelling, 1998; Casulli, 1999; Armfield and Street, 1999). The fractional step method can be outlined as follows: ? Solution of hydrostatic velocity field and free surface ? Solution of non-hydrostatic pressure ? Update of velocity field and free surface to reflect non-hydrostatic pressure. The velocity field is initially approximated via a hydrostatic solution, providing U* and W*, as well as a hydrostatic free surface (?*).  ( ) () 11 22 1 2 1 2 11 22 11 22 ** *n i1 i i,k i,k i,k nn i1,k i,k o i,k b *n i,k i,k *n * * ii i,ki,k i,k i,k b UGUgt x 1 gt z x WGW t UU z x + ++ + ? + + ++ ? +? ?? ??? =???? ? ?? ?? ?? ? ?? ??? ?? ?? = ?? ? ?=?? ? ??? ? ?? ? ? (2.7) where subscripts indicate the cell center (integers) and cell faces (fractions), as seen in Figure 2.2. The x-direction cell length is ?x i , while ?x i+1/2 is the cell center to center distance. The Lagrangian discretization of velocity from the Euler-Lagrange method is represented by GU n , which uses a particle pathline in the time ?n? velocity field to estimate the Lagrangian momentum term at time ?n+1.? This method reduces the artificial damping associated with most low-order methods and can be applied to flows that have a Courant-Friedrich-Lewy (CFL) condition up to two (Hodges, 2000). 8 Chapter 2. Methods 9 The hydrostatic estimated velocity [Equation (2.7)] is subtracted from the discretized non- hydrostatic momentum equation, written as:  ( ) 11 22 1 2 1 2 1 2 11 22 1 2 n1 n1 n1 n i1 i i,k i,k i,k nn i1,k i,k o i,k b n1 n1 nh,i 1,k nh,i,k o i,k n1 n1 nh,i,k 1 nh,i,kn1 n i,k i,k o i,k UGUgt x 1 gt z x PP t x PP t WGW z ++ + + ++ + ? + + ++ + + ++ ++ ++ + ?? ??? =???? ? ?? ?? ? ?? ? ?? ?? ?? ? ? ?? ? ?? ?? ?? ?? ? ? =?? ? ?? ?? ? ? ? (2.8) to obtain velocity correction equations: 11 22 1 2 n1 n1 nh,i 1,k nh,i,kn1 * i,k i,k o i,k PP t UU x ++ ++ ++ + ?? ? ? =?? ?? ?? ? (2.9) 11 22 1 2 n1 n1 nh,i,k 1 nh,i,kn1 * i,k i,k o i,k PP t WW z ++ ++ ++ + ?? ? ? =?? ?? ?? ? (2.10) The discrete form of continuity [Equation (2.2)] is: 11 22 n1 n1 n1 n1 i 1,k i,k i,k 1 i,k i,k i,k UUWW 0 xz ++++ + ++ ?? + ?? = (2.11) which can be applied to the divergence of Equations (2.9) and (2.10) to yield a pressure Poisson equation (in Einstein summation): 2* onh i ii i PU xx t x ??? = ?? ?? (2.12) Equation (2.12) can be written discretely as: i+ 1 / 2 i- 1 / 2 k+ 1 / 2 k- 1 / 2 i,k k+1 k-1 i+1i-1 11 22 11 22 11 1 22 2 n1 n1 n1 n1 nh,i 1,k nh,i,k nh,i,k nh,i 1,k i,k i,k i,k n1 n1 n1 n1 nh,i,k 1 nh,i,k nh,i,k nh,i,k 1 i,k i,k i,k ** * * i,ki,k i,k i,k o i,k i PPPP 1 xx x PPPP 1 zz z UU WW tx z ++++ +? ++++ +? +? + ?? ?? ?? ?? ?? ? ?? ?? ? + ?? ?? ? ?? ? =+ ?? ? ? ,k ?? ?? ?? (2.13) Solution of Equation (2.13) provides a non- hydrostatic pressure, P , which is used to update the horizontal velocity field using Equation (2.9). The updated horizontal velocity field is used to update the free-surface: n1 nh,i,k Figure 2.2: Schematic of grid field + ()11 22 n1 n n1 n1 ,k,k ,k kb t UUz x + ?? ? ++ ??+ ?? ? = ?=? ?? ? ???? ? ?? ? ? (2.14) If Equation (2.10) is used to update the vertical velocity, the resulting velocity field may not be solenoidal due to the residual in the non-hydrostatic pressure solver. To ensure a solenoidal velocity field for mass transport, the vertical velocity is updated diagnostically from the continuity equation: 11 22 11 22 n1 n1 i,k i,k n1 n1 i,ki,k i,k i,k UU WW z x ++ +? ++ +? ??? ?=? ? ? ?? ? (2.15) The free-surface update [Equation (2.14)] has been treated differently in various models. The quasi-hydrostatic models of Mahadevan, et al. (1996a) and Casulli and Stelling (1998) solve only for the hydrostatic free-surface [Equation (2.7)] and do not account for non-hydrostatic effects on surface evolution. Casulli (1999) addressed this problem by correcting the surface elevation for the intermediate step [Equation (2.7)] and after the pressure Poisson equation was applied [Equation (2.14)]. The approach of Casulli (1999) was to advance the velocities, pressure and free-surface at time level ?n+1? while still applying the vertical discretization (?z) of the previous time step (n). His approach removes a nonlinearity that arises when using inconsistent time discretization of ??z? and the free surface, which Hodges (2004) showed exists, but is second-order in time. Thus, the first-order temporal accuracy of CWR-ELCOM makes the inconsistency irrelevant and Casulli?s (1999) treatment appropriate Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models for the updated free-surface. A more complex approach by Chen (2003) applies a double predictor- corrector method, updating the velocity first after the non-hydrostatic pressure is resolved and again after the free-surface is updated. The double predictor- corrector allows the non-hydrostatic change in the free-surface to alter the velocity field. While this method is second-order temporally accurate, Chen (2003) states that a comparison between the results with his method and those of Casulli (1999) show no significant improvement in model performance. 10 2.1.3 Iterative Methods ? Pressure Poisson Equation The pressure Poisson equation [Equation (2.13)] is an elliptic equation, which is inefficient for direct solution by matrix inversion (Mahadevan, et al., 1996a) and is therefore solved with iterative methods. Jacobi, Gauss-Seidel (GS), and successive over-relaxation (SOR) iterative methods were used in the development of the non-hydrostatic solver as a learning exercise for the author. The Jacobi and GS methods were impractical because of long convergence times (i.e. red-black point GS: ~3500 iterations for the L ? norm to converge to 10 -8 in the manufactured solution case described in ?3.1 ff). SOR, with a red-black point iteration scheme, has proved to have adequate convergence times (i.e. red- black point SOR: ~1000 iterations for the L ? norm to converge to 10 -8 ; this reduced the run time by about 70% for the manufactured solution case in ?3.1 ff.). The SOR scheme used is presented in Equation (2.16) with the optimal over-relaxation factor (?) determined by Equation (2.17). () i,kn1 n nh,i,k nh,i,k i,k n1 n1 i 1,k nh,i 1,k i 1,k nh,i 1,k i,k n1 n1 i,k 1 nh,i,k 1 i,k 1 nh,i,k 1 i,k Q P1P a aP aP a aP aP a + ++ ++ ?? ++ ++ ?? ?? =? + ?? ?? ?? ?+ ?? ? ? + ?? ? ? ? ? ? (2.16) g 2 1sin N ?= ?? ? ??+ ?? (2.17) where Q is the source term [i.e. the right hand side of Equation (2.13)], a is the coefficient that represents that spatial discretization and N g is the number of grid cells in the computational domain. Convergence was measured using by the residual (R): (2.18) n1 n1 i,k i,k i 1,k nh,i 1,k i 1,k nh,i 1,k n1 n1 n1 i,k 1 nh,i,k 1 i,k 1 nh,i,k 1 i,k nh,i,k RQaP aP aP aP aP ++ ++ ?? ++ ?? =? ? ?? + When the residual is reduced to a specified convergence criteria (?3.1 ff), the non-hydrostatic pressure is considered resolved and the velocity field and free surface are updated with the new non- hydrostatic pressure. Other methods may be used to solve the pressure Poisson equation more efficiently (e.g. Mahadevan, et al., 1996b; Marshall, et al., 1997; Casulli and Stelling, 1999; Casulli, 1999; Fringer and Street, 2003). Some possibilities include the SOR applied in line-relaxation form, conjugate gradient method or multigrid methods. The multigrid method is well suited to handle stiff problems, such as the pressure Poisson equation, and it may use a Gauss-Seidel iterative scheme, so the non-hydrostatic solver could be adapted into multi- grid form. However, the multi-grid method may be difficult to implement over complex boundaries (He, et al., 1996). 2.1.4 Boundary Conditions The model boundary conditions were assigned so that the vertical solid boundaries perfectly reflect wave propagation and neither vertical nor horizontal boundaries have viscous boundary layers. To achieve these requirements, the velocity is assigned Dirichlet conditions normal to solid boundaries (U = 0 or W = 0) and Neumann conditions tangential to solid boundaries (dU/dx = 0 or dW/dz = 0). Scalar quantities, including the non-hydrostatic pressure, have Neumann conditions (zero gradient) at all solid boundaries. The velocity and scalar conditions are implemented in CWR-ELCOM under the designation of ?free-slip? boundaries. The non- hydrostatic pressure at the free surface is a Dirichlet condition requiring P nh = 0. The boundary conditions for the non-hydrostatic pressure are built directly into the non-hydrostatic solver in the SOR coefficients, ?a,? in Equation (2.16). 2.2 Analysis Methods Both qualitative and quantitative analysis methods are used to compare and contrast hydrostatic and non-hydrostatic model results. The simplest qualitative analysis is visual examination of internal wave evolution over several wave periods in both models and theory. The human eye/brain Chapter 2. Methods combination is capable of distinguishing the overall scope of phase/amplitude errors in an intuitive way that provides the first level of screening for model results. Quantitative analysis methods developed herein are focused on the energy changes within the system. A principle reason for conducting model experiments using the inviscid/diffusionless equations is that any change in the modeled total energy (E T ) is directly attributable to numerical error. Therefore, energy changes are an integrator of numerical dissipation and diffusion errors in the model results. Furthermore, there is an energy shift to smaller wavelengths associated with the degeneration of an internal wave, which is used to identify numerical error. 2.2.1 Background Total energy (E T ) in a system is comprised of kinetic (E K ) and potential energy (E P ), written below for a 2D system with a constant breadth (B): L 2 Ki 0b 1 EUBd 2 ? =? ?? zx P B (2.19) 11 (2.20) L P 0b E g zBdzdx ? =? ?? so that E T may be written as: (2.21) TK EEE=+ The potential energy can further be separated into background potential energy (E B ) and available potential energy (E A ) (Lorenz, 1955). The E B is the potential energy when the system is at its lowest possible energy state. That is, if a system is adiabatically brought to rest and the density field settles without mixing, then the system is considered to have its lowest possible potential energy, which is E B . The E A is the potential energy of a system displaced from its lowest possible energy field, or also described as the E P available to be transferred to E K . Thus, the EA is defined as the difference between the E P and E B : AP EEE? ? (2.22) A simple monochromatic standing wave is an oscillatory exchange of E A (the potential energy in the inclined wave) for E K (kinetic energy when the wave is flat), and back to potential energy again. In an inviscid/diffusionless system, the exchange of E A for E K must be conservative, so it is convenient to define their sum as dynamic energy (E D ): (2.23) DA EEE?+ K The relationship between energies is seen in Figure 2.3 Mass diffusion is the spreading, or weakening, of a stratified density field and is characterized by the increase in the background potential energy (E B ) as heavier particles are diffused upwards and lighter particles diffused downwards. Since E B is the En e r g y Wave Period 01/4 13/41/2 E D E K E A E P E B En e r g y En e r g y En e r g y Figure 2.3: Schematic of energy behavior over a wave period, for a inviscid/diffusionless monochromatic standing wave. The black lines above the graph show the evolution of the wave shape. The purple ( ) line is the potential energy (E P ), the red ( ) line is the background potential energy (E B ), the blue ( ) line is the available potential energy (E A ), the yellow ( ) line is the kinetic energy (E K ) and the green ( ) line is dynamic energy (E D ). The E P is the sum of E B and E A and the E D is the sum of E A and E K . Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models density field when restored to its lowest energy state, it can only change through physical diffusion of the density gradient (Winters, et al., 1995). That is, advection and wave propagation can not directly increase EB. However, in the physical world they may cause shear instabilities and wave breaking which will more readily allow physical diffusion to occur. Resort all water parcels (grid cells) by density, keeping track of the volume of each grid cell. Fill the basin, starting at the bottom with the heaviest parcel, where the volume associated with this parcel is spread evenly over the entire horizontal area. Continue filling the basin from most to least dense, as in Figure 2.5. 12 Numerical diffusion of mass is indistinguishable from physical diffusion of mass in the resulting smoothing and spreading of gradients, which makes it difficult to quantify the numerical diffusion error when the governing equations include physical diffusion. In a model that sets physical diffusion to zero, any change in E B is model error (Laval, et al., 2003), which is commonly referred to as numerical diffusion. Numerical diffusion is dependent on the coarseness of the grid; that is, fine grids reduce numerical diffusion (Hodges, et al., 2000). Thus, as a grid is refined and numerical diffusion is inhibited, modelers typically presume that model representation of E B and internal wave evolution are improved; however, as shall be seen, this is not always true (?4 ff.). Calculate the E P for this resorted system using Equation (2.20), where ??? is the resorted density, ?z? is the height of the layer in the resorted basin, ?dz? is the thickness of the resorted layer and ?dx? is the length of the basin. This is the E B . In the physical world, dynamic energy (E D ) is decreased by: 1) transfer of available potential energy (E A ) and kinetic energy (E K ) to background potential energy (E B ) through turbulent mixing or 2) loss of E K to heat through the action of viscosity. The latter is called ?dissipation,? while the former (mixing) is where waves and advection enhance the rate of physical diffusion. Shear instabilities and wave breaking move a fluid around, (i.e. ?stirring?) which raises some of the heavy fluid, increasing E A and decreasing E K . If physical diffusion is present, a portion of the stirred fluid is irreversibly mixed. Thus, a portion of the E K that went into E A to stir the fluid is, in a sense, converted into E B . Therefore, physical diffusion increases the E B of a system while decreasing the E D of the system. The portion of E K lost due to diffusion is not considered dissipation of energy because the E K was converted to E B rather than heat. Without physical diffusion, stirred fluids cannot irreversibly ?mix? so the heavier fluid would eventually fall back to its lowest energy state, thereby converting the E A back to E K with no impact on E B or E D . Hence, in a diffusionless, inviscid model any decrease in E D represents energy lost in the form of dissipation. As dissipation decreases E K , there is less energy available for conversion to E A and consequently a reduction in the E D of a system. Figure 2.4: Flow chart of sorting algorithm for background potential energy (E B ). The evolution of background potential energy (E B ) and dynamic energy (E D ) provide a measure of numerical diffusion and numerical dissipation. While the total potential energy for any particular density distribution is a relatively simple computation, the calculation of background potential energy is more difficult. Following the example of Winters, et al. (1995), an algorithm was developed to determine E B , as described in Figure 2.4. The model approximates reality by assigning each cell in a system a unique value for density. By example, a 5?5 random domain would result in 25 layers in the resorted system, as in Figure 2.5. Figure 2.5 provides a simple example, but a larger domain results in significantly more layers (i.e. 4380 for a 60?73 grid). The resorted density field yields a vertical resolution that has a much finer scale than the model?s representation of the vertical processes. When the full resorted density field is used in calculations, the results are ?noisy? and are a misleading representation of the model?s resolvable vertical structure. To return the data to the appropriate analysis scale, the resorted density profile is ?binned? into groups of densities to yield a resorted density domain that has the same number of layers as the original vertical resolution (Figure 2.6). Chapter 2. Methods 13 The binned, resorted density field yields a domain that produces smoother calculations (Figure 2.7). 2.2.2 Numerical Diffusion Physical diffusion changes stratification by smoothing density gradients transferring more dense fluid into lower density regions and vice-versa. Typically, changes in the background potential energy provide a measure of the physical diffusion in a system. In a diffusionless system, any changes in the stratification can be attributed to numerical diffusion. Therefore, changes in the background potential energy for a diffusionless system provide an estimate of numerical diffusion. Physical diffusion is modeled by diffusivity for each of the scalars, such as salinity or temperature. Thus, the concept of diffusivity is used to develop a numerical diffusion coefficient based on changes in the background potential energy. The numerical diffusion coefficient acts as a global parameter to compare model performance and error accumulation. The resorted density field allows calculation of the E B per unit area (A): () () B B b Et t,z gz dz A ? =?? ? ? ? ? (2.24) where ? B is the binned, resorted density. The evolution of E B is: () () B B b Et dd t,z gzdz dt A dt ? ?? =? ?? ?? ? (2.25) The free surface for the background density field is constant in a system without inflows or outflows, so the derivative passes through the integral in Equation (2.25). Discretizing Equation (2.25) for ?N? layers yields: ( ) () N B B,k k k k k1 Et dd t,z z gz dt A dt = ?? ? ? ?=?? ?? ? ? ? ? ?? ?? ? ? (2.26) Figure 2.5: Schematic of a random 5?5 density field resorted into 25 layers from most to least dense. The blue indicates least dense and red indicates most dense. The density term for layer ?k? can be represented as: Figure 2.6: Schematic of a random 5?5 density field resorted into 25 layers from most to least dense and then binned into 5 layers. The blue indicates least dense and red indicates most dense. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 0 2 4 6 8 10 12 14 16 18 0 0.05 0.1 0.15 0.2 0.25 t/T E B - E B (0 ) not binned binned E B - E B (0 ) Figure 2.7: Typical change in background potential energy (units: Joules) versus the time normalized by a wave period (T). The green line is before the density field is binned, the blue line is after the density field is binned, with the errorbars representing the standard deviation. () () () () Bk B,k k k k k Bk Bk mt,z dd t,z z z dt dt A z 1d mt,z Adt mt,z A ?? ?? ????= ? ?? ?? ? ?? =? ? ?? =  (2.27) where m?= ? and with ?zAz?= ? k constant in time. The vertical mass flux, , represents the rate of change in the resorted density field; any change in E B m B is due only to the diffusion of mass in the density stratified system. The relationship in Equation (2.27) can be substituted into Equation (2.26) to obtain: () () N BBk k k1 Et mt,z d gz dt A A = ?? = ?? ?? ?  (2.28) Fick?s law states: C qD z ? =? ? (2.29) where D is the diffusion coefficient, q is the mass flux per unit area and C is the mass concentration. The vertical mass flux per unit area in Equation (2.28) may be expressed in the form of Fick?s law: () () () B mz z Dz Az ?? =? ?  (2.30) Making the substitution for the vertical mass flux per unit area into Equation (2.28), yields ( ) () N B kk k1 Et t,z d Dgz dt A z = ???? =? ?? ? ?? ? Bk D (2.31) which can be redefined as N kk k1 F = =? ? (2.32) where ( ) () B Bk kk Et d F dt A t,z gz z ?? = ?? ?? ?? ?=? ? (2.33) F is computed from model results using the rate of change of E B and ? is computed from the resorted density profiles. Local diffusion coefficients (D k ) cannot be directly computed from the global change in E B in Equation (2.32). However, we can estimate an approximate global diffusivity based on a simplified treatment of the density field. In all the present investigations, hyperbolic tangents are used to construct the density field: 14 Chapter 2. Methods 1 py base py py x ez cos z L1 tanh h 2 2 ???? ?? +? ? ? ?????? ???? ?=? + ?? ?? ?? (2.34) where ? base is the density at the center of the pycnocline, ?? py is the density change across the pycnocline, ? is the amplitude, x is the horizontal grid location, L is the length of the basin, z is the vertical grid location and h py is the thickness of the pycnocline. The height at the center of the pycnocline is z py and the initial thickness (h py ) is split equally on each side of the pycnocline center line to designate the initial pycnocline position. Equation (2.34) yields a representative density profile (Figure 2.8). The hyperbolic tangent profile models a continuous stratification in a real world system with three distinct layers; two approximately uniform density layers separated by a region of rapidly changing density. Therefore, dividing a stratified system into three layers: 1 = upper, 2 = pycnocline and 3 = lower (Figure 2.8); Equation (2.32) can be written (2.35) 11 2 2 3 3 DDD?+?+?= 15 F The upper and lower layers have approximately constant density, so ? 1 = ? 2 = 0, which leaves () B 2 2 2 k Et d dt A F D gz z ?? ?? ? ==? ?? ? ? ? (2.36) Equation (2.36) can be approximated as () B 2 13 py py Et d dt A D gz h ?? ?? ? ?? ?? ? ? ?? ?? ?? (2.37) where ? 1 and ? 3 are the densities of the upper and lower layers, respectively. Equation (2.37) yields a global approximation for D across the pycnocline derived from the change in the background potential energy system and changes in the pycnocline thickness and density field. This approximation provides a tool with which numerical diffusion can be estimated for the simple three-layer system. Since this approximation is taken at each time step there is the possibility of oscillations due to small-scale anti- diffusive fluxes. The small-scale anti-diffusive fluxes arise from non-monotonic advection (Leonard, 1991). Oscillations form around sharp wave fronts (i.e. a discontinuity) and are referred to as Gibb?s phenomenon (Kreyszig, 1999). The transport algorithm used in the present research is ULTIMATE QUICKEST; ULTIMATE (Leonard, 1991) smoothes the oscillations in QUICKEST (Leonard, 1979) so that the scalar concentration decreases monotonically in a one-dimensional system. Lin and Falconer (1997) found that ULTIMATE QUICKEST is not monotonic for a two-dimensional system, such as used in the present research. 998 999 1000 1001 1002 0 1 2 3 4 5 6 7 Density (kg/m 3 ) H e i ght ( m ) H e i ght ( m ) Layer 1 Model performance over multiple wave periods is the focus of this work; thus, the change in E B and the pycnocline density gradient are binned over a wave period to smooth anti-diffusive fluxes and provide a clear indication of numerical diffusion. The average form of Equation (2.37) is Layer 2 () B T 2 13 py py T Et 1d dt Tdt A D 1 gz dt Th ?? ?? ? ?? =? ?? ??? ?? ? ?? ? ? (2.38) Layer 3 Figure 2.8: Typical density profile constructed by a hyperbolic-tangent function. The dashed lines separate layers 1, 2 and 3. where T is the wave period. The change in h py is slow if the upper and lower layer densities remain Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models constant and the pycnocline has not diffused over the entire water column (Figure 2.9). With this assumption, Equation (2.38) can be discretized as () () () {} 1 2 1 2 Bn Bn1 n 2 n 13 py py 1 ET ET AT DT gz h ? ? ? ? ? =? ?? ??? ?? ?? (2.39) where the over bar represents the binned data over a wave period. 2.2.3 Numerical Dissipation Energy dissipation in a physical fluid is controlled by molecular viscosity acting at the smallest scales of velocity shear. However, the energy at the smallest scales of motion is fed by the larger scales of motion (Kantha and Clayson, 2000b; Kundu and Cohen, 2002), so the larger scales are considered to be the energy-containing scales and control the rate at which viscosity can work. In modeling, turbulence is commonly discussed in terms of the eddy viscosity associated with a large-scale velocity shear, so that the eddy viscosity characterizes the enhancement of turbulent dissipation for the large-scale features. The concept of eddy viscosity is well-entrenched in our understanding of dissipative phenomenon, so it is useful to develop a measure of a numerical viscosity associated with the dissipation by numerical error. The mechanical energy dissipation rate, typically presented on a unit mass basis (?), (Batchelor, 1967) in a 2D, incompressible fluid is 22 uw u w zx x z ?? ?? ? ?????? ?=? + + + ?? ????? ?? ? ? ????? ?? 2 ? ? ? (2.40) In a diffusionless system, the dissipation rate, integrated over a volume, is equivalent to the rate of change of the E D .  () D 2 o 22 dE d dt uw zx uw 22 xz ? ? =??? ? ???? =?+?? + ? ?? ?? ?? ? ? ? ???? ?? d+ +? ? ?? ?? ?? ?? ?? ? ? ? ? (2.41) Since  ? ? , the density perturbation from the reference density (  ? ) can be neglected, yielding 2 D o 22 dE uw dt z x uw 22 xz ? ? ???? ??? + ? ?? ?? ?? ? ? ? ???? ?? d+ +? ? ?? ?? ?? ?? ? ? ? (2.42) The principle velocity shear (and hence numerical dissipation) occurs across the pycnocline. A characteristic numerical viscosity associated with numerical dissipation is developed that is consistent with the shear characteristics in the pycnocline. Integrating Equation (2.42) over the volume of the pycnocline and substituting the characteristic numerical viscosity of the pycnocline (? py ) for the molecular viscosity (?) yields: 2 D opy 22 py dE uw dt z x uw 22 xz ? ???? ??? +? ?? ?? ??? ? ? ???? ?? + +?? ?? ?? ?? ?? ? ? (2.43) where represents a spatial mean. If the pycnocline volume ( py ? ) is characterized by horizontal area (A) and pycnocline thickness (h py ) 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 t/T h py / H h py / H Figure 2.9: Typical change in pycnocline thickness (normalized by total depth, H) over time (normalized by the wave period, T). This case was scenario 2 for the 30?73 grid (?3.3). 16 Chapter 2. Methods then numerical viscosity, binned over a wave period, may be approximated as () () 1 2 py D n 1 D n n 2 opy 1 22 1 ET ET T uw Ah zx uw 22 xz + ? ? ?? ?= ? ?? ? ? ? ??? ?? ? ?? + ? ?? ?? ? ?? ? ?? ?? ????? ?? ? ++ ? ?? ?? ?? ?? ? ?? (2.44) where an overbar represents the binned data over the wave period. For internal waves, the dominant velocity gradient is ?u/?z, so Equation (2.44) can be reasonably reduced to () () 1 2 py D n 1 D n n 1 2 opy 1 ET ET T u Ah z + ? ? ? ?= ? ? ? ???? ????? ?? ?? ?? ?? ?? ? ?? ???? ? ? (2.45) Equation (2.45) provides an approximation for the numerical viscosity in the model, which is used to compare the relative importance of numerical viscosity with molecular viscosity. 2.2.4 Spectral Analysis As modeled internal waves evolve, dynamic energy is transferred to smaller scale waves or lost to numerical error. The computation of numerical viscosity (?2.2.3) allows quantification of the energy dissipated. In this section we discuss how spectral analysis is used to quantify the transfer of energy into shorter wavelength features. Spectral analysis describes the distribution of signal power over wave frequencies (or wavenumbers). For this work, the spectral analysis decomposes the spatial structure of pycnocline displacement into the power spectral density associated with different wavenumbers (k = 2?/?). A Matlab ? toolbox function (i.e. ?pwelch?) was used to obtain the power spectral density. Welch?s method is used, which is outlined in the following steps: 1) the data is divided into overlapping segments, 2) a Hamming window is used to calculate the modified periodogram for each segment and 3) the modified periodograms are averaged to obtain the estimated power spectral density for the entire data set (Signal Processing Toolbox User?s Guide, 1998). A small data set may reduce the resolution of the power spectral density estimation, thus the data set is replicated to create a longer data set (i.e. in this work, the replicated data set is used to include basin- scale waves). For all analysis in this work, the signal data is divided into eight segments with 50% overlap between adjacent segments and a Hamming window applied on each segment, 32 replicates of the initial half wavelength wave are used (to total 16 wavelengths). The power spectral density is computed at each binned period for all wavelengths, thus the data set is sampled at 2?/?x, which corresponds to wavenumber. As an internal wave steepens and forms solitons, the initial long-wavelength wave evolves into a train of short-wavelength waves. The power spectral analysis provides the power distribution of the internal wave over different wavenumbers. As a wave evolves, an energy shift may occur, transferring energy from the long-wavelength wave to smaller-wavelength waves. For instance, the peak energy at time = 0 coincides with the initial, long- wavelength wave. However, as the wave evolves and dynamic energy (E D ) decreases, there is an energy transfer through different wavelength waves, which appears as a shift in the peak power wavenumber through time. The shift in peak power should coincide with the development of solitons. If a wave does not degenerate into solitons, the peak power should remain with the initial, long- wavelength wave. If a bore develops, then there may be a shift in peak power as the bore provides higher wavenumber components since the wave is no longer sinusoidal. Comparing expected wave behavior with the spectral analysis provides an assessment of model skill, in terms of the model?s ability to represent soliton formation and the evolution of model error. 17 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 18 Chapter 3. Verification and validation Chapter 3. Verification and validation of the non-hydrostatic solver pi/4 pi/2 3pi/4 pi 5pi/4 3pi/2 7pi/4 2pi pi/4 pi/2 3pi/4 pi 5pi/4 3pi/2 7pi/4 2pi 9pi/4 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Verification and validation of a computer code are two necessary and basic steps to demonstrate that the code does what it is intended to do. According to Roache (2002), code verification is an evaluation of error from a known solution, a purely mathematical exercise, while code validation demonstrates the accuracy with which the mathematical model captures the physical phenomena based on theory or measured in the field or laboratory. This chapter verifies the non-hydrostatic pressure solver with an analytical solution, validates the solver?s convergence with a simple test case and validates the model?s solution against theory and a laboratory experiment. Figure 3.1: Solution space for the manufactured solution. The colorbar represents pressure in kg/ms 2 . 3.1 Verification by the Method of Manufactured Solutions The purpose of the manufactured solution is to verify the model?s accuracy and establish the convergence criterion for the numerical model used in this work. 3.1.1 Setup Roache (2002) suggests using the ?Method of Manufactured Solutions? to verify a numerical code. This method uses a continuous mathematical solution independent of the code. The manufactured solution must be non-trivial and should exercise all terms in the numerical code and the corresponding boundary conditions. The manufactured solution is applied to the numerical model?s homogeneous governing equations, producing a non-homogeneous source term that is discretized and added to the numerical model?s governing equations; the resulting discrete non-homogeneous governing equations are used to approximate the exact manufactured solution. This study uses a manufactured two-dimensional analytical solution: (3.1) Pcos(x)cos(z)= Figure 3.1 shows the solution for Equation (3.1). The pressure Poisson equation has the form 22 22 PP Q xz ?? += ?? (3.2) The second derivative of Equation (3.1) is: 2 2 2 2 P cos(x) cos(z) x P cos(x) cos(z) z ? =? ? ? =? ? (3.3) This produces a manufactured source term: (3.4) Q 2 cos(x) cos(z)=? which is used as the discrete source term for the non- hydrostatic solver [e.g. the right hand side of Equation (2.13)] to compute the discrete pressure field of the manufactured solution. Thus, the manufactured solution provides a known solution space which can be compared to the model?s solution driven by the source term, Equation (3.4). This approach allows verification of model performance for different spatial grids and time steps. The verification test domain is a square box, 2? m in length and 7.5 m in height, represented by four different grids, as described in Table 3.1. The boundary conditions of the manufactured solution are the same pressure boundary conditions used in other simulations within the present research (?2.1.4); the domain top (free-surface) has a Dirichlet boundary condition, while all other sides (solid boundaries) have Neumann boundary conditions. 3.1.2 Error Analysis It is necessary to examine the grid convergence error to assess the model?s accuracy at different grid 19 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 20 ) resolutions. Following Roy (2003), grid convergence error is analyzed using the L 1 and L 2 spatial error norms. The spatial error norms consider the entire domain to estimate the order of accuracy of the model. As defined by Roy (2003), the discretization error at grid point (x,z): (3.5) () () ( D grid exact x,z P x,z P x,z?= ? is used in the L 1 and L 2 spatial error norms: () () grid grid N D n1 1 grid 1/2 N 2 D n1 2 grid x,z Lnorm N x,z Lnorm N = ? = ? ? = ? ? ?? = ?? ?? ? ? ? (3.6) Table 3.1: Grids used in the manufactured solution (?3.1). where N grid is the number of grid points in the computational domain. Figure 3.2 shows the L 1 and L 2 spatial error norms. The three different grids used (Table 3.1 cases A - C) were refined in both the horizontal and vertical direction. As the grid is refined, the spatial error norms decrease with second-order behavior, verifying that the non- hydrostatic solver is spatially second-order accurate. In real world problems, analytical solutions typically do not exist. Another way to look at error is the residual [from Equation (3.2)]: 22 grid grid 2 P (x,z) P (x,z) R(x,z) Q(x,z) xz ?? ?? ? 2 Case Grid (nx ? nz) A 10 ? 10 B 20 ? 20 C 40 ? 40 D 41 ? 30 (3.7) Figure 3.2: The L 1 and L 2 spatial error norms of the manufactured solution for grid refinement (Table 3.1 cases A ? C). 10 1 10 2 10 -3 10 -2 10 -1 # grid points in one direction No r m L 1 norm L 2 norm 1 st order slope 2 nd order slope No r m No r m Chapter 3. Verification and validation The L 1 , L 2 and L ? residual norms are then defined: () () () grid grid N grid n1 1R grid 1/2 N 2 grid n1 2R grid Rgrid Rx,z Lnorm N Rx,z Lnorm N L norm max R x, z = = ? = ?? ?? ? = ?? ?? = ? ? ? (3.8) Figure 3.3 shows the L 1 , L 2 and L ? residual norms for increasing number of iterations that the non-hydrostatic solver performs for the manufactured solution. This iteration count is equivalent to the number of iterations per timestep of the pressure solution in an unsteady problem (?3.2). The change in residual norms with number of iterations measures the model?s evolution towards convergence for a specific amount of computational effort. The norms continue to decrease with increasing number of iterations until the residual calculation [Equation (3.7)] reaches machine accuracy. However, the solution at machine accuracy will differ from the analytical solution by the truncation error in the discrete equations (Hirsch, 1988). Convergence is defined when the norms decrease by a specified order of magnitude. Hirsch (1988) gives the example of, ?a second-order accurate space discretization, with ?x = 10 -2 , will produce an error of the order of 10 -4 on the solution, which cannot be reduced further even if the residual equals 10 -14 .? Thus, for second-order discretization of the non-hydrostatic solver, the solution is considered converged when the L ? residual norm decreases by four orders of magnitude. The 1E-12 1E-11 1E-10 1E-09 1E-08 1E-07 1E-06 1E-05 0.0001 0.001 0.01 0.1 1 10 100 1 10 100 1000 10000 # Iterations Nor m norm inf Norm 2 Norm 1 Nor m Nor m Nor m Figure 3.3: Residual norms for manufactured solution (case D, Table 3.1). Solid blue line is L ? norm, dotted blue line is L 2 norm and dashed red line is L 1 norm. The number of iterations are the number of iterations the non-hydrostatic solver performs per timestep. 21 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models manufactured solution drops by four orders of magnitude (Figure 3.3) in about 50 iterations. Figure 3.4 shows the comparison of L ? spatial error and error norms. The L ? spatial error norm flattens at the same number of iterations (~50) where the L ? residual norm has decreased by four orders of magnitude. Thus, the L ? residual norm can be used to reasonably represent the error when it is not possible to calculate the L ? spatial error norm (i.e. where there is no analytical solution). In summary, applying the non-hydrostatic solver to the manufactured solution shows that as the grid is refined, the error decreases and the computed solution converges to the manufactured solution. The spatial error norms verify that the method is second-order accurate. The residual norms are a good proxy for the spatial error norms when an analytical solution is not available. Thus, the L ? residual norm is used to establish the convergence criterion in the above verification of the non- hydrostatic solver; which is the reduction of the L ? residual norm by four orders of magnitude. This convergence criterion is used in all further model simulations. 3.2 Validation and Convergence for an Unsteady Internal Wave The purpose of this test case is to determine the number of iterations to reach a converged solution (i.e. based on ?3.1) in an unsteady internal wave. This test case is sufficiently simple, yet exemplifies many of the nuances of internal wave modeling. 3.2.1 Setup The validation test basin used was 10 m long and 7.5 m deep. The initial wave was a cosine wave with amplitude of 1.125 m, where the upper layer depth to total depth ratio was 0.3. A hyperbolic tangent function [Equation (2.34)] was used to construct the 0.0000001 0.000001 0.00001 0.0001 0.001 0.01 0.1 1 10 100 1 10 100 1000 # iterations No rm residual spatial error No rm Figure 3.4: Spatial error and residual L ? norms for manufactured solution for case D. The red dotted line is where the L ? residual norm decreases by four orders of magnitude (i.e. the convergence criterion). The number of iterations are the number of iterations the non-hydrostatic solver performs per timestep. 22 Chapter 3. Verification and validation density profile (Figure 2.8). The test basin?s density profile has an initial pycnocline that is 5m thick with a density change of 4 kg/m 3 . A square grid of 0.25 m * 0.25 m was applied to this domain. 3.2.2 Results Akin to the manufactured solution in ?3.1, the L 1 , L 2 and L ? residual norms for the internal wave test case decrease with increasing number of iterations that the non-hydrostatic solver performs within each timestep of the simulation (Figure 3.5). However, the decrease in convergence is slower than the manufactured solution (Figure 3.6). The L ? residual norm decreases by four orders of magnitude in about 50 iterations for the manufactured solution case, while the L ? residual norm decreases by the same amount in about 180 iterations for the internal wave test case. The longer convergence time for the test case compared to that of the manufactured solution is discussed in ?6. The estimated root mean square error is used to evaluate the density field for the different number of iterations that the non-hydrostatic solver performs per timestep. The estimated RMS error compares the difference in density at each grid cell within the pycnocline between a simulation with 1000 iterations, which is considered the converged case, against simulations with a lesser number of iterations (i.e. 10, 50, 75, 100, 200 and 500). The estimated root mean square error (? RMS ) of the density field is: () grid 1/ 2 N 2 k,n 1000,n RMS n1 py grid k N = ?? ??? ?? ? ? = ?? ?? ?? ? ? rapid change in the density profile (e.g. Figure 2.8, (3.9) where k is different numbers of iterations the non- hydrostatic solver performs per timestep and ?? py is the density change across the pycnocline. The estimated RMS error is computed only within the pycnocline (e.g. layer 2 in Figure 2.8) because layers 1 and 3 are essentially uniform density. The pycnocline is defined by the region where there is 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1 10 100 1000 # Iterations No rm Norm 1 Norm 2 Norm inf Figure 3.5: Residual norms for the internal wave test case. The solid purple line is the L ? norm, the dotted blue line is the L 2 norm and the dashed red line is the L 1 norm. The number of iterations are the number of iterations the non-hydrostatic solver performs per timestep. 23 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 24 n 00 or s r by in ror drostatic pressure solution requires more n- 1.E-08 1.E-07 1.E-06 1.E-05 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1 10 100 1000 # Iterations No rm Internal wave test case Manufactured solution No rm No rm No rm No rm No rm Figure 3.6: L ? norm for the manufactured solution (case D, dashed blue line) and the internal wave test case (solid red line). The number of iterations are the number of iterations the non- hydrostatic solver performs per timestep. the pycnocline is demarcated between 998.1 kg/m 3 and 1001.9 kg/m 3 ). The RMS error is normalized by the density change across the pycnocline (e.g. in this case, 3.8 kg/m 3 ). Figure 3.7 shows the RMS error values for the six different non-hydrostatic solver iterations per timestep (i.e. 10, 50, 75, 100, 200 and 500). Using 200 iterations per timestep, the RMS error is one or more order of magnitude smaller tha the 10, 50, 75 and 100 iteration cases. After the initial growth in the first 200 timesteps, the error growth rate for the 200 iteration case [O(10 -5 ) in 8 timesteps] is slower than the growth for cases with less iterations [i.e. the 100 iteration case grows to O(10 -3 ) in 800 timesteps]. Of course, the RMS err at 200 iterations is larger than the 500 iteration case, but using 500 iterations increases the computational time (200 iterations: 0.61 s/timestep; 500 iterations: 0.75 s/timestep) without necessarily improving the model results. That is, as discussed in ?3.1.2, the solution cannot reduce the error below truncation a iterations are increased past the convergence criterion (decrease in L ? residual norm by fou orders of magnitude) (Hirsch, 1988). At 200 iterations, the L ? residual norm has decreased four orders of magnitude, so any further decrease the L ? residual norm provided by increasing the number of non-hydrostatic solver iterations per timestep is unnecessary. This estimated RMS er analysis supports the findings from the norm analysis in Figure 3.5. The non-hy computational time. For example, on a Dell Workstation PWS530 Xeon Processor with a 2.4 GHz CPU and 3.5 GB of RAM, the hydrostatic model takes 0.31 seconds/timestep, while the no hydrostatic model takes 0.66 seconds/timestep for Chapter 3. Verification and validation 25 , e ame l, on is for on equation estimates the press ls d s 3.3 Validation Against Regime Theory and This section compares the numerical model?s ves 3.3.1 Regime Theory in referred to as Horn) ry , 2) that oretical imescales prov ning awn the same timestep and grid resolution where the convergence criteria of the non-hydrostatic solve reducing the L r is ? norm by four orders of magnitude; for a grid with 60 horizontal grid cells ? 73 vertical grid cells (0.01 m horizontal resolution ? 0.004 m vertical resolution) and a timestep of 0.012 seconds simulating 10 wave periods requires 90000 timestep iterations, the hydrostatic model takes 7.9 hours, while the non-hydrostatic model takes 16.4 hours. However, the above mechanics understates the scop of the problem: as discussed in ?4 ff., the non- hydrostatic model requires a finer timestep and spatial resolution than the hydrostatic model to capture the resolvable processes, therefore a more meaningful comparison is the computational time required for a coarse-resolution hydrostatic model to a fine resolution non-hydrostatic model. While the choice of the coarse-resolution grid is somewhat arbitrary, for illustration we might consider a hydrostatic solution on a 30 ? 73 grid on the s problem, with a 0.2 s time step and 1800 timestep iterations. The resulting solution time is 0.09 hours (5.5 minutes). With this example, the difference in computational expense between the hydrostatic and non-hydrostatic models is evident; the spatial and temporal requirements of the non-hydrostatic mode combined with the pressure Poisson solution, render it a huge computational undertaking. The solution of the pressure Poisson equation accounts for 89% of the increase in simulation time. Updating the velocity and free-surface field increases the simulation time by 10%. The above discussi a two-dimensional solution; a three-dimensional solution may further increase the difference in computational effort. The pressure Poiss 0 200 400 600 800 1000 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 # Timesteps 10 50 75 100 200 500 ? RMS ure at a point using the pressure from the surrounding cells. Other non-hydrostatic mode (e.g. Marshall, et al., 1997; Casulli, 1999) use the pressure from time ?n? as the initial estimate of the pressure field for time ?n+1,? however the present research found that the time ?n? pressure may be a poor approximation for the time ?n+1? pressure fiel due to the internal wave propagation (?6.1 ff.). Setting the pressure field to zero at each timestep provides a pressure Poisson solution that converge faster than solution with the previous timestep?s pressure field (Figure 3.8). Figure 3.7: Root mean square error (? RMS ) for density field in the test case. Each line represents the ? RMS for a different number of iterations per timestep. The density range is 3.8 kg/m 3 . Laboratory Experiment representation of internal wave evolution to the evolution of several laboratory-scale internal wa and theoretical wave evolution. Horn, et al. (2001; here compared internal wave evolution in a laborato experiment with four theoretical timescales of different internal wave phenomena: 1) damping steepening, 3) shear instability (Kelvin-Helmholtz billow formation) and 4) internal bore formations. The theoretical timescales are a function of the basin and internal wave dimensions, such as length, height and amplitude; Figure 3.9 shows a schematic of a basin. From the theoretical timescales, Horn developed a regime diagram (e.g. Figure 3.10) uses the theoretical timescales to determine the prevailing wave phenomena under different conditions. Horn tested and validated the the regimes in a laboratory experiment. Regime diagrams based on the t ide the boundaries between damping, steepe into solitons, and formation of Kelvin-Helmholtz billows or undular bores. A single ?universal? regime diagram (i.e. for all waves) cannot be dr as the regime boundaries are functions of wave amplitude, upper layer thickness, total depth, basin Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models Figure 3.8: L ? norm for the non-hydrostatic solver when the pressure at the present timestep is approximated from zero (solid blue line) and the previous timesteps solution (dashed red line). 26 ata in rn, Regime 1 in Figure 3.10 is the ar n rnal d, while the wave evolution in regimes 1 and 2 ct the d for n?s Laboratory Experiment and Model Simulations series of model simulations: length l length, pycnocline thickness, and the density change across the pycnocline. While there are many possible ways to graph this multidimensional d set, Horn showed that a graph of wave amplitude to upper layer depth ratio (a/h) and upper layer depth to total depth ratio (h/H) provides a clear delineation of the timescales and regimes. Figure 3.10 is an example of a regime diagram for the basin used Horn?s experiments. As discussed in Ho region where viscous damping dominates nonlinear steepening, typically for internal waves with small amplitudes. Regime 2 is where nonline steepening is dominant and solitons are formed. Regime 2 flows have been observed many times i different field studies (e.g. Boegman, et al., 2003; Farmer, 1978; Wiegand and Carmack, 1986). Regime 3 is supercritical flow, with large wave amplitudes and a shallow pycnocline, which are associated with internal hydraulic jumps and inte bores. Regime 4 physics are dominated by Kelvin- Helmholtz billows, which predominantly occur in shallow systems with deep pycnoclines; they are a result of shear instabilities at the interface between layers. Regime 5 physics are characterized by large amplitude internal waves that develop into bores and may include Kelvin-Helmholtz billows. Waves developing in the high-mixing regimes 3, 4 and 5 (supercritical flow, Kelvin-Helmholtz billows and bores) have the onset of principle phenomena occurring within one-quarter of the wave perio (damping and soliton formation) typically requires one or more wave periods to significantly affe wave characteristics. Horn?s regime diagram provides a framework to validate the non-hydrostatic solver. Simulation experiments were performe several scenarios corresponding to the different regimes. 3.3.2 Hor The dimensions of Horn?s laboratory experiment were used for a (L) = 6 m, height (H) = 0.29 m and width = 0.3 m (Figure 3.9). The regime diagram in Figure 3.10 is defined for the above scales. Horn?s experimental basin is a closed system, with solid boundaries on al sides and the top. The present model is designed 0.0001 0.001 0.01 0.1 1 10 100 1 10 100 1000 # iterations Nor m P = Previous Pressure P = 0 H-h h H a ? L ? 2 ? 1 Figure 3.9: Schematic of basin. Chapter 3. Verification and validation 27 1 m) inic f rn?s tank diffe 3 hown in Table 3.2; these a / h h / H Figure 3.10: The regime boundaries from laboratory experiments; ordinate: amplitude-upper layer depth ratio, abscissa: upper layer depth-total depth ratio. A typical interface thickness, 1 cm, was used to determine the timescales, T KH and Td. The star represents Kelvin-Helmholtz billows and bore, the diamond is the broken undular bore, the triangle is the solitons, the square is steepening, and the circle is damped linear waves. (from Horn, et al, 2001) only for free-surface simulations, so this is one area where the model and experimental conditions diverge. However, the free-surface remains essentially flat (maximum displacement of 0.00 as the initial conditions are an entirely barocl flow, which has no significant coupling to the free surface (Kantha and Clayson, 2000b). The model boundary conditions are as described in ?2.1.4. The model viscosity is set to zero, which is a second area where the model and experimental conditions diverge. However, as discussed in ?2.1.1, this only effects the damping rate of modeled waves at fine grid scales when the numerical dissipation is less than the physical dissipation. Thus, it can be expected that the model may show development o solitons (regime 2) under conditions when Ho experiments show damped linear waves. Modeling the viscosity in Horn?s experiment is impractical as the principle viscous damping in the experiment is not the shear across the wave interface, but is the boundary layers on the sides and top of the tank. Modeling the drag of these boundary layers is impractical with the present model and of little relevance to the present subject of interest. Horn?s experiments were set up in a tilted with two fluids of different density. The density rence between the two layers was 20 kg/m across a pycnocline thickness of 0.01 to 0.02 m. The tilted tank was quickly moved to a horizontal position to initialize the wave. This impulsive start creates an initial condition that is approximately a linear tilt. Model simulations of Horn?s experiments are initialized with a cosine wave having an amplitude that corresponds to the tank?s initial angle of tilt in Horn?s experiment. The difference in initial wave setup is the third area of divergence between the model simulation and the laboratory experiment. However, this difference is not expected to affect the wave evolution. A three layer hyperbolic tangent salinity profile [Equation (2.34)] was used to establish the initial density profile that is close to Horn?s condition (Figure 3.11). Nine scenarios were modeled to reproduce some of Horn?s experiments, as s scenarios were chosen based on Horn?s reported results (see Horn, et al., 2001 Figures 5 and 6). These scenarios exemplify the influence of nonlinearity, and the subsequent development non- Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 28 o ?s not available for numerical qualitative comparisons are made to f n aboratory essary to review the principal diffe osity in e ary el the e hydrostatic pressure gradients, on internal wave evolution. Group A maintains the same depth ratio (h/H) and varies the amplitude ratio (a/h) to illust the relationship between nonlinearity and initial wave amplitude. Group B changes the depth ratio and amplitude ratio. Changing the depth ratio als influences the nonlinearity of the system; that is, as the depth ratio decreases, the nonlinearity increases. Scenarios 6 and 9 have a depth ratio of 0.5. At this depth ratio, the wave is considered weakly nonlinear. Therefore nonlinear steepening and soliton formation are not significant wave evolutions. Scenario 6 has a small amplitude ratio (a/h = 0.18), so the expected phenomenon is a damped wave. Scenario 9 has a large amplitude ratio (a/h = 1.0) and Kelvin- Helmholtz billows are expected to form. 3.3.3 Results rate Horn raw data was comparisons, so figures from Horn, et al. (2001). Specifically, the model simulations are compared qualitatively against Horn?s laboratory experiment to demonstrate the model?s ability to capture soliton development in Figures 3.12 - 3.14. The model is compared to the experimental results by comparing the evolution o pycnocline displacement at a point in the tank; the amplitude of the wave and solitons, time of emergence of solitons and number of solitons withi a wave series are the three major points of comparison. Before a discussion of the model and l results, it is nec rence between the model simulations and Horn?s experiments; that is the neglect of visc the model simulations. The experimental tank is narrow and has a lid, so there is significant area (i.e. the surface area of the interface is 1.8 m 2 and the surface area of the tank is 3.8 m 2 ) that creates a viscous boundary layer which will dominate the wave damping. The computational domain in th model simulations is the same size, but the bound layers are not modeled. Without viscosity, the mod simulations have no physical damping to restrain wave steepening. The model does have numerical dissipation which may cause some damping, however, the numerical viscosity due to numerical dissipation is less than molecular viscosity for wave periods up to the emergence of solitons (?4.1 ff.). Therefore, as steepening is occurring, the wav experiences negligible damping, so the modeled wave tends to develop solitons sooner and with greater amplitudes than in the experiments with viscous effects. Wave Evolution Interfacial displace wave evolution. T ment can be used as a measure of he interfacial vertical Table 3.2: Scenarios used in the laboratory-scale case (?3.3) Scenario h/H a/h Group A 10.300.15 20 3 30.300.45 40 6 50.300.9 Group B 60.500.18 70423 20.300. 80245 90.501.0 10 15 20 25 30 0 5 10 15 20 25 He i g h t ( c m) Model Horn He i g h t ( c m) He i g h t ( c m) He i g h t ( c m) He i g h t ( c m) He i g h t ( c m) He i g h t ( c m) He i g h t ( c m) Figure 3.11: Initial density profile for Horn?s laboratory experiment (red, dashed line) and the modeled hyperbolic tangent density profile. Density anomaly 3 Chapter 3. Verification and validation 29 tory data to l ped waves. Phys e ent; ing the le each . teepen and develop solitons, however Horn ing ed orted o rved by Horn: the n ut the e billo is 13 elvin- the displacement recorded by the central wavegauge in Horn?s tank is used to compare the labora the isopycnal displacement in the center of the mode domain. Figure 3.12 shows the model and Horn?s results for group A (?3.3.2., Table 3.2). In general, the model matches the qualitative evolution of wave phenomena and the period of the waves, but over- predicts the wave amplitude. This is consistent with comparing an inviscid model to a laboratory experiment with viscous fluid. There tends to be one less soliton in the wave series in the model simulations than in the laboratory experiments; this will be discussed later in this section. Scenario 1 is within regime 1 where the dominating wave behavior is viscous dam ical viscous damping is not present and th damping due to numerical dissipation is less than physical damping. Therefore, it is expected that modeled waves within regime 1 will have less damping (i.e. these wave will have slightly larger ?bumps? develop) than in the laboratory experim this is seen in a comparison of scenario 1 between the model simulation and Horn?s experiment. Scenarios 2-4 are in regime 2 with soliton development, which the model represents, show ading wave with the largest amplitude and successive wave in the train with smaller amplitudes The model results qualitatively follow Horn?s trend of more rapid soliton formation with increasing nonlinearity. Scenario 5 is in regime 2 and should theoretically s observed a broken undular bore in the lead wave (Figure 3.12). The wave evolution, as described by the model simulation, shows bore development (Figure 3.15), but the breaking described by Horn is not observed. While detail data on the breaking phenomena were not rep by Horn, it is likely their spatial scales are too fine t be captured in the present model. Behind the bore, the model shows soliton formation, which is also indicated in Horn?s results. The model simulation for scenario 5 shows slightly fewer solitons in the wave series than Horn?s experiment. The model results for group B (?3.3.2., Table 3.2) show the same characteristics obse onlinearity of internal wave evolution increases as the amplitude ratio (a/h) increases and as the depth ratio (h/H) decreases (Figure 3.13). With increased nonlinearity, non-hydrostatic pressure effects increase, which allows the non-hydrostatic model simulations to depict soliton formation. Witho non-hydrostatic pressure, the modeled evolution of the wave is quite different (?4 ff.). Similar to the results for group A, the model simulations of group B generally match the qualitative evolution of wav phenomena and the period of the waves seen in Horn?s laboratory experiment. The wave amplitude is over-predicted in the model simulations. Scenario 9 has a large amplitude ratio (1.0) and is on the boundary between Kelvin-Helmholtz ws (regime 4) and bore formation (regime 5). The timescale for Kelvin-Helmholtz formation s and the timescale for bore formation is exactly one quarter of the wave period (25 s). Bores and Kelvin- Helmholtz billows form when their respective timescales are less than one quarter of the wave period (Horn, et al. 2001). For this scenario, K Helmholtz billows were observed by Horn, while model simulation showed the wave degenerate into a system of higher mode waves (Figure 3.16). Kelvin- Helmholtz billows are a fine-scale phenomenon, which require a small grid resolution to model them. The grid scale (60?73) used may be too coarse horizontally to capture this event. Emergence of Solitons A critical point of c of solitons. Horn define omparison is the emergence s the emergence of solitons as, ?? s oliton the time when the waves are sufficiently well separated that the depth (measured from the crest of the leading wave) of the trough between the leading solitons (measured from the crest of the leading wave) is 25% of the amplitude of the leading wave? (Horn, et al., 2001). Horn used three wavegauge spaced approximately 1.5 m apart. As solitons may have emerged between wavegauges, Horn reports the time it took for the wave to move between the wavegauge that showed solitons and the upstream wavegauge. Horn?s definition and method of s emergence was used to identify the emergence of solitons for the model simulations. The period of the emergence of solitons for Horn?s experiments and model simulations is seen in Figures 3.12 and 3.13. Horn defines the timescale of steepening as: s L T = ?? where L is the basin length, ? is the amplitude and (3.10) 12 hh3 c 12 2hh ? ?= (3.11) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 30 Scenario 1: h/H = 0.30, a/h = 0.15 Scenario 2: h/H = 0.30, a/h = 0.30 Scenario 3: h/H = 0.30, a/h = 0.45 Scenario 4: h/H = 0.30, a/h = 0.60 Scenario 5: h/H = 0.30, a/h = 0.90 0 50 100 150 200 250 300 350 400 -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) -5 5 ? (c m ) -5 5 ? (c m ) Time (s) Figure 3.12: Interfacial displacements for Group A at the center of the tank. The figures on the left are Horn?s experiment (Horn, et al., 2001). The figures on the right are the model simulations. In both sets of figures, the dotted lines indicate the range in which solitons emerged. Scenario 8: h/H = 0.20, a/h = 0.45 Scenario 2: h/H = 0.30, a/h = 0.30 Scenario 7: h/H = 0.40, a/h = 0.23 Scenario 6: h/H = 0.50, a/h = 0.18 -2 2 ? (c m) -2 2 ? (c m ) -2 2 ? (c m ) 0 50 100 150 200 250 300 350 400 -2 2 ? (c m) Time (s) Figure 3.13: Interfacial displacements for Group B at the center of the tank. The figures on the left are Horn?s experiment (Horn, et al., 2001). The figures on the right are the model simulations. In both sets of figures, the dotted lines indicate the range in which solitons emerged. Chapter 3. Verification and validation 31 where h 1 is the upper layer thickness, h 2 is the lower layer thickness, H is the total depth, c is the wave speed: 12 hh cg' H = (3.12) g? is reduced gravity. T ng was calcu rted in Table 3.3. The model simulations have soliton formation slightly earlier than the timescale of steepening and Horn?s observations (Figures 3.14). This discr scid approxim Numb he timescale of steepeni lated for each scenario and are repo epancy can be attributed to the invi ation in the model. er of Solitons Scenarios 3-5 had less solitons emerge in the model simulation than in Horn?s laboratory experiment. The reason for this is unknown. The number of solito the difference in wave am time of soliton be .1 ff.). ?73) ngth (12 m) by 120 grid cells. The author?s hypothesis that this resolution would be sufficient is incorrect. There needs to be more investigation into the lower limit of resolution needed to capture all of the solitons in the wave series. However, the grid resolution must be balanced with viscosity since horizontal grid refinement decreases the numerical viscosity to the same order of magnitude or smaller than molecular viscosity (?4.1 ff). With decreased numerical Figure 3.14: Comparison of timescale of steepening with Horn?s observation of the emergence of solitons (blue lines with circle) and model simulation observation of the emergence of solitons (red lines with square). ns is not a function of viscosity, like plitude and emergence. The number of solitons appears to dependent on the horizontal grid resolution (?4 The grid resolution examined in this section (60 represents the initial wavele 0 50 100 150 200 250 300 350 400 450 0 50 100 150 200 250 300 350 400 450 T S (s) So l i ton em ergence (s ) 0 50 100 150 200 250 300 350 400 -5 5 ? (c m ) ? Figure 3.16: Isopycnal displacement for scenario 9 at the center of th (c m ) e tank. h/H = 0.5 and a/h = 1.0. Table 3.3: Timescale of steepening for scenarios in internal wave test case (?3.3). Scenario Ts (s) 1 427 2 213 3 142 7 446 8 124 9-- 4 107 571 6-- 1 2 3 4 5 6 0.29 Length (m) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) enario 5, bore development at time 60s. Figure 3.15: Sc Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 32 scosity, the modeled wave amplitude and time of Summary The manufactured solution (?3.1) verified that the in e solver iterations per timestep needed for co s ex on- h hysics of internal wave evolution. Non-hydrostatic pressure has a dispersive effect on internal wave evolution. The non-hydrostatic model clearly demonstrates the emergence of solitons for a wave that lies within regime 2. The major difference between Horn?s laboratory experiment and the model simulations is the viscous orn?s experiment are not included in the model simulations. Horn?s experiment has a large boundary layer area compared to volume of the basin, thus viscous effects are prevalent in the laboratory results. The model simulations did not model viscosity, so solitons emerged earlier and soliton amplitudes were larger. However, the basic characteristics of the wave train were captured by the model. The model was orrectly represent the number of solitons in the wave series at t on- h tested a e n of i vi emergence of solitons may become significantly different than what is predicted by theory and shown in laboratory experiments. 3.4 effects in H non-hydrostatic solver is accurate and convergent. The convergence criterion was established to be when the L ? residual norm is decreased by four orders of magnitude. The internal wave test bas (?3.2) validated the convergence criterion established by the manufactured solution and determined the number of non-hydrostatic pressur not able to c he grid resolution used. In summary, the n ydrostatic model is shown to be valid when gainst theory and a laboratory experiment. Th on-hydrostatic model captures the physics nternal wave evolution. nvergence. Comparison between Horn? periments and model results validates the n ydrostatic solver by showing it captures the p Chapter 4. Comparison of models 33 odel the physics of (a non-hydrostatic process), yet van (2004) observed significant id c model nd ares the laboratory expe rostatic and non-hydrostatic model simulations. One laboratory-scale internal wave ase is examined on different grid meshes to examine the effect of grid resolution on model skill for the hydrostatic and non-hydrostatic models. The lake-scale simulation is also analyzed on different grid meshes. The analytical techniques in ?2.2 are used to quantitatively evaluate changes in energy, which is used to compare differences between hydrostatic and non-hydrostatic models on the different grid meshes. 4.1 Laboratory Scale Comparison Results of hydrostatic and non-hydrostatic models have been compared with the same laboratory experiments (Horn, et al., 2001) used in ?3.3. Non- hydrostatic pressure is a small-scale effect; thus, it is advantageous to model a laboratory-scale internal wave because the small-scale grid (e.g. 0.1 m ? 0.004 m) allows the non-hydrostatic pressure to affect wave evolution and not be damped out due to the size of the grid. The hydrostatic and non- hydrostatic models used eight different grid resolutions to compare model skill for one scenario (scenario 2, Table 3.2). The grid aspect ratios (?z/?x) for all resolutions were of O(10 -2 ) or greater. 4.1.1 Setup Basin dimensions and the scenarios are identical to those listed in ?3.3.1 and Table 3.2. Each scenario was separately simulated using the both hydrostatic The timestep is dependent on the physics of the wave and the grid resolution. Due to the characteristics of the hydrostatic and non-hydrostatic models, the timestep required to yield a stable solution is different. For both models, the timestep selection is based on the Courant-Friedrich-Lewy (CFL) condition for baroclinic motions. The hydrostatic model uses the form of the CFL condition: Chapter 4. Comparison of Hydrostatic and Nonhydrostatic models Hydrostatic models do not m wave dispersion Hodges and Dela wave dispersion in a hydrostatic model which produces soliton-like formations. Such soliton-like formations may mimic the expected internal wave evolution, however this coincidence is a false positive. Without non-hydrostatic pressure, the soliton-like formations can only be an artifact of the numerical method, and is thus dependent on gr resolution and numerical truncation error. Previous literature has not examined the effect of grid resolution on hydrostatic and non-hydrostati performance. This chapter examines laboratory a lake scale simulations. The laboratory-scale simulations examine nine different internal waves (?3.3) and qualitatively comp riments to the hyd c and non-hydrostatic models (?4.1.2). For all scenarios, a model grid of 60?73 (?x = 0.1 m and ?z = 0.004 m) was used. The pycnocline is represented by five cells (h py /?z = 5). Scenario 2 is modeled on several different grid meshes, as described in Table 4.1, to compare model error (?4.1.3). Each simulation for the changing grid resolution was run for ten wave periods as this time allows the internal wave to evolve and develop into solitons. All simulations used the inviscid/diffusionless approximation. x t CFL c x ? ? ? # Cells in Z-dir 15 29 73 ? o /?x 30 x x x 60 # Cells in X-dir 60 x x x 120 600 x x 1200 h py /?z 125 Table 4.1: Grid meshes for scenario 2 simulations Those meshes with an ?x? were performed. . (1.13) where ?t is the timestep for a specific, initial internal wave speed (c) and grid resolution (?x). The CFL condition is a non-dimensional number which limits the timestep and grid discretization needed to stabilize the conditionally stable explicit Euler scheme (Ferziger and Peri?, 2002) used in CWR- ELCOM. The baroclinic CFL condition is the most restrictive condition, as opposed to the barotropic or advective CFL condition, for density-stratified flows (Hodges, 2000). In CWR-ELCOM, the maximum allowed baroclinic CFL x = ?2 (Hodges, 2000), but this work uses a conservative CFL x condition of 1/3. The hydrostatic model is dissipative (?4.1.3 ff.), so strong vertical motions are inhibited. The CFL x condition in Equation (1.13) provides a sufficient timestep for the hydrostatic model. However, the non-hydrostatic model is generally less dissipative, especially for finer horizontal grid resolutions (?4.1.3 ff.), so strong vertical motions may occur. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 34 hus, the timestep that was appropriate for the hydrostatic model is not suitable for the non- e odel l FL T hydrostatic model (i.e. strong vertical motions caus instabilities, Figure 4.1) and the timestep must be smaller. The timestep for the non-hydrostatic m is determined by the vertical velocity and vertica grid resolution, yielding a new definition for the C condition: z ac t CFL xz ? ? ?? (1.14) where a is the wave amplitude and ac/?x is a measure of vertical velocity of a bore. The maximum CFL the theoretical number of solitons (N) that will evolve in z condition allowed is 1/3. The smaller timestep simulation does show a marked improvement in the internal wave evolution in non-hydrostatic model (Figure 4.1), and is thus used for further analysis. The hydrostatic model is not significantly altered by the different timesteps. 4.1.2 Scenario Simulations This section uses qualitative comparisons between the laboratory experiment and the hydrostatic and non-hydrostatic model simulations. The number of solitons that emerge in a wave train for the laboratory experiment and the model simulations is compared to the theoretical number of solitons that will evolve. Kao, et al. (1985) defined the a train as: M N1? + ? (1.15) where M is defined for a two-layer system in a model with the Boussinesq approximation: 1 2 12 o22 12 2 hh hh3 ML ? ?? =? ? ? ? ? The length of the basin is L, h (1.16) and the 1 is the upper layer thickness, h 2 is the lower layer thickness and ? o is the wave amplitude. The results are separated into the simulations for group A, group B, scenario 9 number of solitons that emerge in a wave train. Group A Group A has the same interface depth ratio 0.3) and varies the a (h/H = mplitude of the wave for each scena olitons ic model ormation. Scenario 5 has the greatest initial amplitude and forms an undular bore. The boratory experiment shows signal, indicating that the bore was initially turbulent (Horn, et al., 2001). The non-hydrostatic model iment he rio (Figure 4.2). At small amplitudes (scenario 1), there is some small steepening but no s emerge; both models show similar evolution. For scenarios with larger amplitudes (regime 2), the nonlinearity of the wave increases causing the wave to steepen and evolve into a train of solitons with decreasing amplitude. The non-hydrostatic model represents this evolution, while the hydrostat showed bore f la a high-frequency shows the bore formation, but does not show a high- frequency signal or any other indication of turbulence in the wave evolution. The non- hydrostatic model was able to characterize the soliton train that formed in the laboratory exper after the bore was damped by turbulent mixing. T hydrostatic model developed the initial bore with a slightly smaller amplitude than the non-hydrostatic model. As time progresses, the hydrostatic model retains the bore shape in the leading wave and develops soliton-like features behind it. Group B Group B varies the interface depth ratio for each scenario (Figure 4.3). For scenarios with larger depth ratios, that is a thick upper layer, the nonlinearity of the wave is small. When the depth ratio is at mid-depth (h/H = 0.5, scenario 6), the Figure 4.1: Scenario 6 pycnocline displacement simulated with a coarse timestep, determined by CFL , and a fine timestep, determined by CFL . x z 0 50 100 150 200 250 300 350 400 -4 -2 0 2 4 Time (s) ? (c m ) coarse ?t fine ?t ? (c m ) ? (c m ) ? (c m ) Chapter 4. Comparison of models Figur l disp he lef erim 2001). The f ddle hydro ulations. Th d lines in Horn nd the d. 35 lacements for Group A measured at the center of the tank. The figures on t are the non-hydrostatic model simulations. The figures on the right are the non-hydrostatic model indicate the range of time in which solitons emerge e 4.2: Interfacia igures in the mi ?s experiment a t are Horn?s exp static model sim ent (Horn, et al., e dotte Scenario 1: h/H = 0.30, a/h = 0.15 2: h/H = 0.30Scenario , a/h = 0.30 3: h/H = 0.30Scenario , a/h = 0.45 4: h/H = 0.30Scenario , a/h = 0.60 0.Scenario 5: h/H = 30, a/h = 0.90 0 50 100 150 200 250 300 350 400 -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) -5 5 ? (c m ) -5 5 ? (c m ) Time -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) ? (c m ) -5 5 ? (c m ) ? (c m ) ? (c m ) -5 5 ? (c m ) ? (c m ) ? (c m ) ? (c m ) 0 50 100 150 200 250 300 350 400 -5 5 ? (c m ) ? (c m ) ? (c m ) ? (c m ) ? (c m ) ? (c Tim Horn?s experiment Non-hy ydrostatic drostatic model H model m ) e (s) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models system has no appreciable nonlinearity and the hydrostatic and non-hydrostatic model produce similar results. A thin upper layer (i.e scenario 8) causes the wave to steepen quickly and results in a train of solitary waves with larger amplitudes (a ~ 2.5 cm) than a thicker upper layer (i.e. scenario 7, a ~ 1.7 cm), 36 h is seen in the non-hydrostatic model. The hydrostatic model is able to predict the onset of steepenin ut the wave evolution is entirely unphysic Scenario whic g, b al. 9 Horn, et al. (2001) observed shear instabilities developing into Kelvin-Helmholtz billows in only one experiment (scenario 9). There is no figure of laboratory observations of this experiment, and it is not specified how the wave evolves after the billowing event, so no direct comparison can be made between the laboratory experiments and model simulations. Neither the hydrostatic, nor the non- hydrostatic models were able to capture Kelvin- Helmholtz billows. Horn, et al. (2001) observed two bores form at both ends of the tank and propagate towards the center, which the non-hydrostatic model did simu The non-hydrostatic model shows the bore pro tion, with a series of solitons behind the bore (Fig 4.4). The hydrostatic model has neither bore form n, nor Kelvin-Helmholtz billows. Number olitons in Wave Train late. paga ure atio of S In gener e non-hydrostatic model underestimates the number of observed solitons in the wave series from the laboratory experiment (?3.3.3). The theoretical number of solitons, as defined by Kao, et al. (1985), for each scenario is in Table 4.2. Figure 4.5 shows how the number of solitons in a train for the model simulations and laboratory experiments compare with the theoretical number of solitons. The theoretical number of solitons is always larger than what was observed in the laboratory experiment and the model simulations. Summa al, th ry -hyd rator s (Fi hows ry ex The non tic mode oduces results of the labo eriment o differen scenario s 4.2 and . The n ydrostati model s tons emerge sooner th n the laborato ments an slightly grea plitude isopycnal dis ement. T can be attrib to the application o viscid approximation as discusse in ?3.3.3. The hydrostatic model captures the evolution for weakly nonlinear waves (scenarios 1 and 6), but presents a poor ntation of nonlinear wave evolution ( 2-5, 7 and 8; Figures 4.2 and 4.3). In sc 2-5, 7 and 8, the hydrostatic model exactly m the non- hydrostatic model ste causes the hydrostatic model vel a bore while the non-hydrostatic model evol o a train of solitons. Bore dev ment in the hydrostatic model is typically followed by so on-like formations, but these are d from observations in the laboratory experimen non-hydrostatic model simulations. The hy ic soliton-like formations are near the cre bore front and their amplitudes are small c d to the amplitudes of the soliton tr the laboratory experiments. A compariso ernal wave evolution for scenario 2 is ure 4.6. This confirms observations of h and non- hydrostatic internal wave e om the isopycnal displacement. Discussion represe scenarios enarios atches epening op into ves int me solit ifferent ts and drostat st of the ompare ains in n of int seen in Fig ydrostatic volution fr until to de elop rosta y exp gure soli peri ter am plac uted f the in l repr quite well f 4.3) d in his d the r the on-h an i t c Observations of the isopycnal displacement demonstrate that the non-h model performs better than hydrostatic mod spect reproducing the int wa on of the laboratory experim by . (2001). The non-hydrostatic model sho fferences in internal wave evolution am fferent scenarios, while the hydros l either shows a damped wave (scenario 1 a bore with small soliton-like formations (sce 7 and 8). The soliton-like formations beh re do not disperse or change in ampl e progresses; this observation is the same narios where bore formation occurs. Th ostatic model shows the train of solitons nd the leading wave and decrease de with time, as in the laboratory experimen drostatic pressu a dispersive pro g, 1972); the non-h static model?s ab ulate disper verifies that the static model captures the overarching ph onlinear wave evolu Tab 2: Theoretical number olitons (N) that w p for a 2-layer, Bou sq mode Simus12345 89 N 6.2 8.3 9.9 11.3 13.6 6 12.7 1.0 ydrostatic el, with re ve evoluti Horn, et al ws clear di ong the di tatic mode nd 6) or a nario 2-6, ind the bo itude as tim for all sce e non-hydr disperses behi in amplitu t. Non-hy perty (Lon ility to sim non-hydro ysics of n ill develo 67 1.0 5. to ernal ents re is ydro sion tion. of sle 4. ssine lation l. Chapter 4. Comparison of models 37 Figure 4.3: Int ial displ figu The figures in t middle ar ght a Horn?s experi and the tons erfac he ment acements for Group B measured at the center of the tank. The e the non-hydrostatic model simulations. The figures on the ri non-hydrostatic model indicate the range of time in which soli res on the left are Horn?s experiment (Horn, et al., 2001). re the hydrostatic model simulations. The dotted lines in emerged. Scenario 8 = 0.2 Scenari = 0. , a/h = 0.30 Scenari 0.23 Scenari 0.18 : h/H o 2: h/H o 7: h/H 0, a/h = 0.45 30 = 0.40, a/h = o 6: h/H = 0.50, a/h = -2 2 ? (c m) -2 2 ? (c m ) -2 2 ? (c m ) 0 50 100 150 200 250 300 350 400 -2 2 ? (c m) Time (s) Horn?s experiment ro drostatiNon-hyd static model Hy c model -2 2 ? (c m ) ? (c m ) ? (c m ) ? (c m ) -2 2 ? (c m ) ? (c m ) ? (c m ) -2 2 ? (c m ) Time 0 50 100 150 200 250 300 350 400 -2 2 ? (c m ) ? (c m ) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 38 0 50 100 150 200 250 300 35 0 2 4 6 8 10 12 14 16 12345678 Simulation # o f S o li t o n s Theory Laboratory Model # o f S o li t o n s Figure 4.5: Number of solitons that develop according to theory (Ka al., 198 d in the laboratory experiment and model simulations. 1 2 3 4 5 6 0.05 0.15 0.25 Length (m) He i g h t ( m ) Figure 4.6: a) Hydrostat d b) non rostatic el of sce 2 internal wa on. 1 2 3 4 5 6 0.05 0.15 0.25 Le (m) He i g h t ( m ) a) b) ngth ve evolutinario mod-hydic an 5) an et o, 0 400 -5 5 ? (c m ) ? (c m ) Figure 4.4: Isopycna placement for scenario 9 at the center of the tank. h/H = 0.5 and a/h = 1.0. l dis Chapter 4. Comparison of models Wave evolution in scenario 9 includes the appearance of shear instabilities in a Kelvin- Helmholtz billow (Horn, et al, 2001). While neither the hydrostatic nor the non-hydrostatic models were able to capture the billows. However, the non- hydrostatic model was able to simulate the bore forma 39 Horn, et al. (2 evolves after the Kelvin-Helmholtz billow collapses, but it is speculated that the wave will degenerate into a train of solitons. The non-hydrostatic model shows a series of solitons develop behind the leading bore. If the speculation of the wave evolution is correct, than the non-hydrostatic model may be a reasonable estimation of the wave evolution after the initial Kelvin-Helmholtz billow formation and collapse. 4.1.3 Changing the Grid Resolution dynamic energy ev the power spectral Isopycnal Displacement tion propagating from the end of the tank. 001) does not specify how the wave Scenario 2 was examined on several different grids (Table 4.1) to asses the effect of grid resolution on the models? performance. Scenario 2 lies within regime 2 and is expected to have steepening and soliton formation. The analysis techniques discussed in ?2.2 are used to quantify differences between the models and grid resolutions. Specifically, the different analysis methods used to quantify model results are the isopycnal displacement, background potential energy evolution, numerical diffusivity, olution, numerical viscosity and density. The isopycnal displacement of scenario 2 for the different grid resolutions is seen in Figures 4.7 and 4.8. The wave evolution for a grid refined vertically (Figure 4.7) shows that the leading wave moves slightly faster for finer vertical grids. The amplitudes of all the solitons within the wave train are smaller for finer vertical grids (~ 1 cm difference 7). hydrostatic models. Horizontal refinement substantially changes the pycnocline displacement (Figure 4.8). Horizontal grid refinement also develops differences in the wave period. As the grid is horizontally refined, the hydrostatic and non- hydrostatic models predict a smaller wave period; that is in Figure 4.8, the finer horizontal grid has an average wave period of 120 s, while the coarser grid has an average wave period of 230 s. Furthermore, the effect of grid scale on the wave influences the formation of soliton trains. The non-hydrostatic model has more solitons in the wave train for the finer grid than the coarser grid. Finally, the leading soliton has approximately the same wavelength (~1 m) through time and for the different grids in the non-hydrostatic model. In the hydrostatic model, the wave evolves into a bore followed by soliton-like formations; the wavelengths of the soliton-like formations change in size as the grid is refined. between the coarsest and finest grids in Figure 4. This is seen in both the hydrostatic and non- Background Potential Energy Changes in the background potential energy (E B ) represent numerical diffusion in a diffusionless model (?2.2). The non-hydrostatic and hydrostatic models show an increase in the background potential 0 100 200 300 400 500 600 700 800 -10 10 ? (c m ) Time (s) 30x15 30x29 30x73 0 100 200 300 400 500 600 700 10 800 -10 Time (s) ? (c m ) Non-hydrostatic a) Figure 4.7: Isopycnal displacement of scenario 2 for vertical grid refinement, measured at the center of basin. Hydrostatic b) a) Non-hydrostatic model, b) hydrostatic model. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 40 at (Table 4.3). In the odels, refinement ows less change in E 73 and 60?73) is about an order of magnitude less than rse 15 and 60?15). The system can be considered unphysical when the normalized background potential energy (Figure 4.9) goes above unity. In the non-hydrostatic model, the E B goes above unity around three wave periods he ve ins physical for only nt vertical grid has more effective numerical diffusion than a fine vertical grid, so grid refinement reduces the numerical diffusion. Numerical Diffusivity energy (E B ), with the hydrostatic model increasing a slightly lower rate (Figure 4.9). In both models, grid refinement in the horizontal direction (Figure 4.9 c and d) shows more change in E B over time. In the non-hydrostatic model, the finest horizontal grid (600?29) has a change in E B of about 92% over nine wave periods, while the coarsest grid (30?29) has a change in E B by 80% (Table 4.3). The hydrostatic model has a change in E B by 89% for the finest grid and 77% for the coarsest grid hydrostatic and non-hydrostatic m in the vertical direction sh B (Figure 4.9 a and b). The finest vertical grid (30? for the coarsest grid (30?15). Refining the grid horizontally for this vertical resolution decreases t time when the E the coa st vertical grid (30? B goes unphysical to 2.2 wave periods. For the 600?29 grid, the E B increases abo one around five wave periods, while the 30?29 grid is physical until about nine wave periods. The hydrostatic model has slightly less diffusion than the non-hydrostatic model and rema a short time longer. Conversely, vertical refineme significantly improves the numerical diffusion in both models (Figure 4.9); at ten wave periods the E B is well below unity for the finest vertical resolution, irrespective of the horizontal resolution. A coarse Non-hydrostatic Hydrostatic 0 50 100 150 The numerical diffusion coefficient (Figure 4.10) derived in ?2.2 confirms what was observed in the background potential energy (E B ); vertical refinement reduces the system?s diffusivity, while refinement in the horizontal increases the system?s diffusivity (Figure 4.10). The hydrostatic model has lower diffusivities than the non-hydrostatic model for all grids, except 60?73. As time progresses, there is an initial growth in the diffusivity, followed by a continual decline. This trend matches the behavior of E B , which increases rapidly at first and flattens with time (Figure 4.9). The numerical diffusivity is orders of magnitude larger than the usion is significant in the models and physical diffusion (in diffusivity of salt, indicating that numerical diff 200 250 300 350 400 -10 10 Time (s) ? (c m ) 0 50 100 150 200 250 300 350 400 -10 10 Time (s) ? (c m ) 30x73 60x73 a) b) Figure 4.8: Isopycnal displacement of scenario 2 for h ic mod Grid Hydrostatic Non-hydrostatic 30*15 90 91 60*15 79 90 600*15 77 80 30*29 89 92 % Increase orizontal grid refinement, measured at the center el. of basin. a) Non-hydrostatic model, b) hydrostat 60*29 78 85 600*29 77 80 Table 4.3: Percent increase in background potential energy for scenario 2. Chapter 4. Comparison of models e 4.9: Background potential energy (E 41 Figur cenario al grid , B) Horizontal grid = 60, C) Vertical grid = 15, esent hydrosta rostati odel results (h). The difference between E B and t B i mic ener B ) of s D) Vertical grid = 29. Lines repr non- he initial E is normalized by the init al dyna 2 for different grid resolutions. A) Horizont tic model results (nh), markers represent hyd gy (E = 30 c m D ). 0 2 4 6 8 10 -2 10 -1 10 0 10 1 t/T 3 3 3 3 3 3 0x15 nh 0x15 hyd 0x29 nh 0x29 hyd 0x73 nh 0x73 hyd 0 2 4 6 8 10 -1 10 0 10 1 t/T 60x15 nh 60x15 hyd 60x29 nh 60x29 hyd 60x73 nh 60x73 hyd a) b) [E B ? E B (0 ) ] / E D (0) 0 2 4 6 8 10 -1 10 0 10 1 t/T 30x29 nh 30x29 hyd 60x29 nh 60x29 hyd 600x29 nh 600x29 hyd d) [E B ? E B (0 ) ] / E D (0) [E B ? E B (0 ) ] / E D (0) 0 2 4 6 8 10 10 -1 10 0 10 1 30x15 nh 30x15 hyd 60x15 nh 60x15 nh 600x15 hyd 600x15 nh t/T c) [E B ? E B (0) ] / E D (0 ) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 42 Fig 10: Numerical diffusivity (?) of scenario 2 ifferent grid resolutions. a) Horizontal grid = 3 b) horizontal grid = 60, c) vertical grid = 15, d) vertical grid = 29. Lines represent non-hydrostatic mod results (nh), markers represent hydrostatic model results (h). ? s = 10 -9 m 2 /s. ure 4. for d 0, el 0 2 4 6 8 10 0 10 1 10 2 10 3 10 4 t/T ? / ? s ? / ? s a) 30x15 nh 30x15 hyd 30x29 nh 30x29 hyd 30x73 nh 30x73 hyd 0 2 4 6 8 10 0 10 1 10 2 10 3 10 4 t/T ? / ? s ? / ? s 30x15 nh 30x15 hyd 60x15 nh 60x15 hyd 600x15 hyd 600x15 nh c) 0 2 4 6 8 10 0 10 1 10 2 10 3 10 4 t/T ? / ? s ? / ? s b) 60x15 nh 60x15 hyd 60x29 nh 60x29 hyd 60x73 nh 60x73 hyd ? / ? s ? / ? s 0 2 4 6 8 10 0 10 1 10 2 10 3 10 4 30x29 nh 30x29 hyd 60x29 nh 60x29 hyd 600x29 nh 600x29 hyd t/T d) Chapter 4. Comparison of models 43 e absence of turbulence and wave breaking) would ve a negligible effect on internal wave evolution. Dynamic Energy th ha Decreases in dynamic energy (E D ) represent the dissipation by numerical error (?2.2). The non- hydrostatic model has lower dissipation rates than the hydrostatic model for all grid resolutions (Figure 4.11). Horizontal grid refinement (Figure 4.11 c and d) decreases the rate at which E D is dissipated. For the finest horizontal grid resolution (600?29), the non-hydrostatic model dissipates E D by 8% over nine wave periods, whereas the hydrostatic model dissipates E D by 85%. The coarsest grid (30?29) dissipates E D by 80% in the non-hydrostatic model and the hydrostatic model dissipates E D by 90% (Table 4.4). Vertical grid refinement (Figure 4.11 a and b) generally increases the dissipation rate. However, the 600?29 grid has less dissipation than the coarser 600?15 grid, which is counter to other observations of vertical grid refinement (Figure 4.12). Numerical Viscosity Numerical viscosity (?2.2) tends to be smaller for r resolution grids (Figure 4.13); the non- rostatic model produces numerical viscosities that are smaller than those produced by the hydrostatic model. Generally, the numerical viscosity increases and then decreases as time progresses. The decrease in numerical viscosity with time corresponds to an overall decrease in dynamic energy (Figure 4.11). Once the dynamic energy begins to decrease, the velocity shears that drive the numerical viscosity are being reduced. That is, any system with a lower dynamic energy should generally lose dynamic energy at a lower rate than a system with a higher dynamic energy, and thus a lower numerical viscosity is expected, which is seen in Figures 4.11 a and b and 4.13 a and b. The numerical viscosity is slightly larger than molecular viscosity; as the grid is refined, the numerical viscosity decreases and nears the molecular viscosity. Numerical viscosities computed for the hydrostatic model are slightly greater than those for the non-hydrostatic model, which is expected as the hydrostatic model has greater numerical dissipation than the non-hydrostatic model. Power Spectral Density fine hyd evolution is associated with the wavelength that has the maximum pycnocline displacement (i.e. amplitude). For all grid resolutions, the hydrostatic model had no shift in peak power from the initial wavelength within ten wave periods (Figure 4.15). However, there is still some energy transferred to different wavelength waves. Timeslices of the power spectral density show as the wavelength decreases there is energy associated with certain wavelengths that follow a harmonic pattern (Figure 4.16). The non-hydrostatic model does show a shift in peak power between four and seven wave periods depending on the grid resolution (Figure 4.15). Finer horizontal grids shift the peak power sooner than coarser horizontal grids. The peak power shifts to different wavelengths, although the range is small (~ .12 ? 0.17 wavelengths). Once the peak power shift has taken place, all grid resolutions for the non- hydrostatic model show similar peak power/wavelength evolution. For example, as seen in Figure 4.15d, the 600?29 grid shifts the peak power at t/T = 4 to ?/? o = 0.12. At t/T = 4.5, the 600?29 mesh shifts to ?/? o = 0.15, which coincides with the peak power shift of the 60?29 grid. Prior to the peak power shift, the non-hydrostatic model shows a harmonic energy shift to ?/? o = 0.5 and some smaller wavelengths (Figure 4.17). However, after the peak power shift, the peak power is located at a wavelength ratio (?/? o ) of 0.2 for all grid resolutions. The spectral analysis shows that the peak power shift is to 20% of the initial basin-scale wavelength, which is consistent with the degeneration of the wave into a train of five solitons (Figure 4.2). The increase in power at small wavelengths shows that energy is being directly transferred from the basin-scale wave to the train of solitons. Grid Hydrostatic Non-hydrostatic 30*15 85 77 60*15 85 47 600*15 80 20 30*29 91 80 60*29 90 52 600*29 89 8 % Increase Table 4.4: Percent decrease in dynamic energy for scenario 2. 0 The power spectral density (?2.2) quantifies the energy transfer of different wavelength waves over time (Figure 4.14). The peak power in a wave Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 44 Figure 4.11: Dy = 60, c) Vertical grid = 15, d) V l results (h). E D is normalized b namic energy (E D ) of scenario 2 for different grid resolutions, a) Horizontal grid = 30, b) Horizontal grid ertical grid = 29. Lines represent non-hydrostatic model results (nh), markers represent hydrostatic mode y the initial E D . 0 1 2 3 4 5 6 7 8 9 10 -1 10 0 30x15 n 30x15 h 60x15 n 60x15 h 600x15 h 600x15 n h yd h yd yd h t/T 0 1 2 3 4 5 6 7 8 9 10 10 -1 10 0 t/T 60x15 60x15 60x29 60x29 60x73 60x73 nh hyd nh hyd nh hyd 0 1 2 3 4 5 6 7 8 9 10 -1 10 0 30x15 nh 30x15 hyd 29 nh 29 hyd 73 nh 73 hyd 30x 30x 30x 30x t/T c) E D / E D (0) a) b) E D / E D (0) E D / E D (0) d) E D / E D (0) 0 1 2 3 4 5 6 7 8 9 10 10 -1 0 t/T 30x29 nh 30x29 hyd 60x29 nh 60x29 hyd 600x29 nh 600x29 hyd Chapter 4. Comparison of models Discussion 45 The key findings of the laboratory-scale model simulation are: 1) in both the hydrostatic an n- hydrostatic models horizontal grid refinement always leads to more growth in the background potential energy and larger numerical diffusivities, indicating the horizontal grid refinement increases the numerical diffusion of the system, 2) the hydrostatic model is more dissipative and somewhat less diffusive than the non-hydrostatic model and 3) the non-hydrostatic model is nearly free of numerical dissipation for a sufficiently refined horizontal grid and the controlling mechanism for model skill is numerical diffusion. The increase in background potential energy with horizontal grid refinement is opposite of the conventional wisdom that reducing grid size must improve model skill and reduce error. This finding corroborates the finding of Hodges and Delavan (2004) and may be explained by examining the concept of energy exchange from ?2.2. Dissipation decreases the dynamic energy, which is the kinetic and available potential energy. The available potential energy is the energy that activates diffusion (when mixing occurs). Thus, if there is less available potential energy to be mixed, then there is less possible diffusion. As the grid is horizontally refined, the dissipation rate decreases. Re ing merical dissipation allows the backgrou ntial energy to continue to grow, thereby increasing numerical diffusion. The non-hydrostatic model has less dissipation than the hydrostatic model (Figure 4.11), and thus has more diffusion (Figure 4.9). Vertical grid refinement generally has less numerical diffusion and more numerical dissipation in both models, which is tied to bove idea that the more dynamic energy is dam the less energy available to be diffused. The results indicate that non-hydrostatic model numerical diffusion is s dependent on the vertical grid resolution (F 9 a and b and Figure 4.10 a and b) and num on is strongly dependent on the hor grid (Figure 4.11 c and d and Figure 4.13 c . That is, as a grid is vertically refined, num iffusion decreases and horizontal refin decreases numerical dissipation. Indeed merical dissipation in fine horizontal such that the numerical viscosity is on the o molecular viscosity. It is an expensive t a nearly numerically dissipative free m .g. a grid with 1200 cells per wavelength; fo sin a horizontal grid of 600). However, consi urbulent eddy viscosities where it is suppose ddy viscosity is 100 times the molecular viscosity, all of the grids examined would have smaller ical viscosities than the eddy viscosity (Figur . Ev the eddy viscosity is 10 times the lar viscosity, several grids (30?73, 60?73, 00?15 and 600?29) would have numerical viscosities less than the eddy viscosity. This shows that coarser grid resolutions may still produce l that performs well with respect to numerical ity when compared with eddy viscosity peculated that with numerical viscosities low eddy viscosity, then the numerical dissipation e 4.11) could be considered small. Following, l with low numerical dissipation would h ? skill controlled by numerical diffusion. The non-hydrostatic mo ntation of internal wave evolution closel the evolution seen in the laborato ent. The wavelength and amplitude is tely the same for different resolutions. The odel also shows a shift in peak pow haracteristic wavelength of the solitons. T ccurs at all grid resolutions and is an indi he non- hydrostatic model?s ability to e physics of internal wave evolution. The model was not able to always captur of observed solitons occurring in train. This difference is most likely due t with numerical dispersion. Numerical dispersion is Figure 4.12: Dynamic energy of scenario 2 for the non-hydrostatic model. Six grids are displayed where the horizontal refinement (horizontal grid = 30, 60 and 600) is the same for two different vertical grids (vertical grid = 15 and 29). E D is normalized by the initial E D . 0 2 4 6 8 10 -1 10 0 D D () 30x15 60x15 600x15 600x29 60x29 30x29 t/T D D () E D / E D (0) the a ped, for the trongly igure 4. erical dissipati izontal and d) erical d ement , the nu grids is rder of ask to have odel (e r this ba dering t d the e numer e 4.13) molecu 60?15, 6 a mode viscos . It is s er than (Figur a mode ave its del?s represe y matched ry experim approxima non-hydrostatic m er to the c his shift o cation of t capture th non-hydrostatic e the numbers the wave o an issue en if d no duc nd nu pote Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 46 0 2 4 6 8 10 -1 10 0 10 1 10 2 t/T ? / ? m ? / ? m 30x15 nh 30x15 hyd 30x29 nh 30x29 hyd 30x73 nh 30x73 hyd 0 2 4 6 8 10 -1 10 0 10 1 10 2 t/T ? / ? m ? / ? m 60x15 nh 60x15 hyd 60x29 nh 60x29 hyd 60x73 nh 60x73 hyd 0 2 4 6 8 10 -2 10 -1 10 0 10 1 10 2 t/T m ? / ? m ? / ? 30x15 nh 30x15 hyd 60x15 nh 60x15 hyd 600x15 hyd 600x15 nh 0 5 10 -2 0 2 4 6 8 10 12 14 16 18 20 t/T ? / ? m ? / ? m 30x29 nh 30x29 hyd 60x29 nh 60x29 hyd 600x29 nh 600x29 hyd Figure 4.13: Numerical viscosity (?) of scenario 2 grid = 60, c) vertical grid = 15, d) vertical grid = markers represent hydrostatic model results (h). viscosity includes negative values. ? for 29. L , Note, 2/s. a) different grids. a) Horizontal grid = 30, b) horizontal ines represent non-hydrostatic model results (nh) d) the ordinate is not log scale because the numerical m = 10 -6 m b) c) d) Chapter 4. Comparison of models affected by discretization methods in time, which were not analyzed in the present research; it remains a po 10 -1 10 0 10 0 normalized wavelength (? / 2L) n o r m al i z ed t i m e ( t / T ) -7 -6 -5 -4 -3 -2 Figure 4.14: Representative power spectral density for non-hydrostatic model of scenario 2. The colorbar represents the log 10 relative PSD. 47 interest for future research. The hydr odel did not show the soliton development observed in the laboratory experiment. Analysis of the power spectral density shows the hydrostatic model retains most of the wave?s energy within the bore (represented as the initial wavelength wave in Figure 4.15). Small soliton-like formations were observed behind the bore, but the peak power is not shifted to them (Figures 4.15 and 4.16). The soliton-like formations vary for different horizontal grid resolutions in both amplitude and wavelength. The soliton-like formations are considered an artifact of the model and grid resolution because the soliton- like formations differ significantly from the laboratory experiment, vary with the grid resolution and the peak power is not transferred to them. It is thought that the hydrostatic model will continue to change the soliton train characteristics with grid refinement. 4.2 Lake-Scale Comparison The hydrostatic and non-hydrostatic models were app mod labo dev Zhou scale problems. Those that do (Marshall, et al., 1997; Casulli, 1999), do not model internal waves. Th testing the behavior of a non-hydrostatic m on a larger-scale internal wave has not previously been done. Similar to ?4.1.3, the grid is varied and the analytical methods of ?2.2 are used to quantify the effect of grid resolution on the model?s skill. 4.2.1 Setup A 2D closed basin with length (L) of 12500 m and height (H) of 50 m was chosen as representative of a medium-sized lake (Delavan, 2003). The initial wave shape was a cosine and the density field was constructed with a three-layer hyperbolic tangent function [Equation (2.34)]. The density difference across the pycnocline was varied to assess the effect of density gradient on model performance. The amplitude ratio (a/h) was 0.3 and the depth ratio (h/H) was 0.3. The wave used in this simulation has amplitude and depth ratios that indicate steepening and soliton formation should occur (i.e. regime 2). Several different grid meshes were examined, as listed in Table 4.5. The timestep used was based on Equation (1.13); this will be discussed further in int of ostatic m us, odel lied to a lake-scale comparison to demonstrate el performance on a larger scale than the ratory experiment. Most of the previously eloped non-hydrostatic models (e.g. Stansby and , 1998; Chen, 2003) do not attempt real-world ?4.2.2 ff. Each simulation was conducted over 100 wave periods. The inviscid/diffusionless approximation was used. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 48 0 0.5 1 0 1 2 3 4 5 6 7 8 ?/? o t/ T t/ T 30x15 nh 30x15 hyd 30x29 nh 30x29 hyd 30x73 nh 30x73 hyd 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 ?/? o t/ T t/ T 60x15 nh 60x15 hyd 60x29 nh yd60x29 h 0 0.5 1 0 1 2 3 4 5 6 7 8 9 ?/? o t/ T t/ T 30x15 nh 30x15 hyd 60x15 nh 60x15 hyd 600x15 hyd 600x15 nh 0 0.5 1 0 1 2 3 4 5 6 7 8 9 10 ?/? o t / T t / T 30x29 nh 30x29 hyd 60x29 nh 60x29 hyd 600x29 nh 600x29 hyd Figure 4.15: Wavelength with peak Horizontal grid = 30, b) horizon grid = 29. Lines represent non-hydro represent hydrostatic model results (h PSD o tal grid = 6 static s yd). f scenario 2 for different grids. a) 0, c) vertical grid = 15, d) vertical model results (nh), marker Chapter 4. Comparison of models 49 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 PS D ?/? o 30x29, t/T=2.5 60x29, t/T=2.5 30x29, t/T=7.5 60x29, t/T=7.5 600x29, t/T=2.4 600x29, t/T=7.3 PS D Figure 4. and t/T = Figure ulation of scenario 2 at t/T = 2.5 and t/T = 7.5. 16: Typical timeslices of PSD for hydrostatic model simulation of scenario 2 at t/T = 2.5 7.5. 4. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 17: Typical timeslices of PSD for non-hydrostatic sim ?/? o PS D 30x29, t/T=2.5 60x29, t/T=2.5 30x29, t/T=7.5 60x29, t/T=7.5 600x29, t/T=2.4 600x29, t/T=7.3 PS D Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models Results Results are analyzed for nine simulations using the loss of dynamic energy (Figure 4.19), computed numerical viscosity (Figure 4.20), change in background potential energy (Figure 4.21), computed numerical diffusivity (Figure 4.22) and the peak power spectral density (Figure 4.23). As a note to the results, the finest horizontal grid went completely unstable for all vertical resolutions (81?9, 81?27, 81?81), as seen in the E B , E D , numerical diffusivity and adve 50 viscosity in Figures 4.19 - 4.22. The vertical ctive CFL: advective z t CFL w z ? ? ? ? (1.17) was greater than one, indicating the presence of strong vertical motions. The Euler-Lagrange Method is stable for hydrostatic models (Hodges, 2000), but the stability has not been investigated in this research or elsewhere in literature for non-hydrostatic models. Reducing the timestep to adhere to Equation (1.14) did not mitigate the instability. The addition of viscosity had no effect on the instability either. This issue is discussed further in ?4.2.3. The internal wave initially steepens with small, soliton-like features developing (Figure 4.18) and is damped out after ten wave periods. The system is dissipative, losing an order of magnitude of dynamic energy in less than ten wave periods for all grid resolutions (Figure 4.19). The large numerical dissipation accounts for the damping effect observed on the wave?s evolution. For vertical and horizontal grid refinement, there is nearly the same loss in dynamic energy for the first three wave periods. After this time, grid refinement causes more dynamic energy to be dissipated. The most significant differences in dynamic energy amongst grid resolutions were due to refinement in the vertical tory- scale hydrostatic simulations (i.e. the dynamic energy decreases about an order of magnitude in ten wave periods). However, in the laboratory-scale # Cells in Z-dir 92781 ? o /?x 9xxx18 # Cells in X-dir 27 x x x 54 81 nh nh x 160 h py /?z 139 grid. The decrease in dynamic energy in the lake- scale simulations is comparable to the labora Table 4.5: Grid meshes for lake scale comparison. Those meshes with an ?x? were performed; nh indicates that only the non-hydrostatic model was used. 2000 4000 6000 8000 10000 12000 10 20 30 40 Length (m) He i g h t ( m ) t/T = 0.9 He i g h t ( m ) 2000 4000 6000 8000 10000 12000 10 20 30 40 Length (m) He i g h t ( m ) t/T = 2.9 Figure 4.18: Internal wave evolution for lake-scale case at t/T = 0.9 and t/T = 2.9. Chapter 4. Comparison of models 51 Figure 4.19: Dynamic energy of lake-scale for different grid resolutions. a)Vertical grid id rid = 9, b) vertical gr 81, d)horizontal g 27. Lines represent n results (nh), markers model results (h). = 27, c) vertical grid = = 9, e) horizontal grid = on-hydrostatic model represent hydrostatic 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (0 ) t/T 9 x 81 nh 9 x 81 h 27 x 81 nh 27 x 81 h 81 x 81 nh 81 x 81 h E D / E D (0 ) E D / E D (0 ) 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (0 ) t/T 9 x 9 9 x 9 9 x 2 9 x 2 9 x 8 9 x 8 nh h 7 nh 7 h 1 nh 1 h E D / E D (0 ) E D / E D (0 ) 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (0 ) t/T 27 x 9 nh 27 x 9 h 27 x 27 nh 27 x 27 h 27 x 81 nh 27 x 81 h E D / E D (0 ) E D / E D (0 ) 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (0 ) t/T 9 x 9 x 27 x 27 x 81 x 9 nh 9 h 9 nh 9 h 9 nh E D / E D (0 ) E D / E D (0 ) 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (0 ) Tt/ 9 x 27 nh 9 x 27 h 27 x 27 nh 27 x 27 h 81 x 27 nh E D / E D (0 ) E D / E D (0 ) d) e) c) b) a) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models non-hydrostatic model simulations had less numerical dissipation than the hydrostatic model and horizontal grid refinement dominated changes in numerical dissipation, which was not seen in the lake-scale simulations. After about ten wave periods, the numerical viscosity (Figure 4.20) becomes scattered. At ten wave periods the wave is damped, so there are small velocities and subsequently low numerical dissipation. Therefore, after ten wave periods, small fluctuations in the dynamic energy are amplified resulting in scatter. This scatter is irrelevant to the present research as we are interested in wave evolution before damping. All grids show the numerical viscosity initially increase and then decrease after two or three wave periods (Figure 4.20). This is similar to the laboratory-scale results (?4.1.3). Vertically refined grids (Figure 4.20 d 52 e) showed lower values for the numerical viscosity, similar to the laboratory-scale results; there was not a significant difference in the numerical viscosity for horizontally refined grids, which was not seen in the laboratory-scale results. One other difference between the lake-scale and laboratory- scale results for numerical viscosity is that the numerical viscosity for the lake-scale simulations were over 100 times greater than molecular viscosity, while the laboratory-scale simulations? numerical viscosity never exceeded 30 times molecular viscosity. The growth in background potential energy is greater for horizontal grid refinement and less for vertical grid refinement (Figure 4.21) similar to the laboratory-scale. Vertical grid refinement dominates changes in numerical diffusion over different grid resolutions. The normalized change in E B increases above unity and is considered unphysical around three wave periods for the coarsest vertical grid. The finest vertical grid has a change in E B from ~ 0.01 to 0.4 over 100 wave periods (Figure 4.21 d and e). Changes in the horizontal grid have a relatively small effect on the growth of E B . Likewise, the numerical diffusivities (Figure 4.22) decrease by an order of magnitude for vertical grid refinement, while horizontal refinement produces relatively similar numerical diffusivities (Figure 4.22). This confirms the vertical grid as the controlling mechanism of numerical diffusion, which is the same result in the laboratory-scale simulations (?4.1.3). The power spectral density (Figure 4.23) shows a shift in peak power around ten wave periods; the wavelength with peak power varies with the grid resolution. A coarse horizontal grid (e.g. Figure 4.23a) shows all power shifting to the same wavele (?/? o ~ 0.5) for different vertical finer horizontal resolution, as the efined (Figure 4.23 b), the peak aller wavelengths (?/? o = 0.2) and ack towards larger wavelengths e progresses. However, unlike the imulations where the peak power e wavelengths irrespective of the e peak power shifts to different the grid is horizontally refined the non-hydrostatic model tical results to the hydrostatic o grid resolutions that showed en the hydrostatic and non- ere the 27?9 and 9?27 grids. Figures 4.21 and 4.22 show that the 27?9 and 9?27 grids had less diffusion in the non-hydrostatic model than the hydrostatic model. After ten wave periods, the background potential energy in the non- hydrostatic model converges to the background potential energy in the hydrostatic model. The numerical diffusivities in the non-hydrostatic model are less than the numerical diffusivities for the hydrostatic model for the entire 100 wave period simulation. The density gradient across the pycnocline was varied (?? = 0.1, 1.0 and 10.0 kg/m 3 ) to asses its effect on the modeling of internal waves. There is no significant difference between hydrostatic and non- hydrostatic simulations for the identical grid resolution (Figure 2.24). 4.2.2 Discussion Long (1972) states that there is a relationship between nonlinearity and non-hydrostatic pressure; a sufficiently steep wave has significant non- hydrostatic effects affecting the internal wave evolution. The wave simulated in this section is nonlinear and lies within regime 2 (a/h = 0.3, h/H = 0.3). Thus, non-hydrostatic pressure effects should disperse the steepened wave and have different numerical error behavior than the hydrostatic model. However, the results show the non-hydrostatic model is nearly identical in all grid resolutions to the hydrostatic model. The similarity between the hydrostatic and non-hydrostatic models can be explained by grid resolution. Grid resolution plays a key role in how the non-hydrostatic pressure affects internal wave evolution. The model grid can only represent a single, average non-hydrostatic pressure for a single cell. Thus, for cells with small grid aspect ratios (i.e. 9?9: aspect ratio = 4[10 -3 ]; 9?81: ngth resolutions. For a grid is vertically r power shifts to sm gradually moves b (?/? o = 0.5) as tim laboratory-scale s shifted to the sam grid resolution, th wavelengths when (Figure 4.23 c-e). Application of generates nearly iden model. The only tw any difference betwe hydrostatic models w and Chapter 4. Comparison of models 53 essure n ion (in his author, y ith ted stability may be due to the d r l., If only grid mesh. Stabi vant to real-world scale basins. nt and sure tic rse linearly steepening initial wave odel model is numerically e lake tio the is )] aspect ratio = 4[10 -4 ]), the non-hydrostatic pr effect is muted. The 27?9 grid has the largest grid aspect ratio (0.012) of all the grids successfully applied to the lake-scale wave. This grid resolutio did show small differences in numerical diffus terms of the background potential energy and numerical diffusivity) between the hydrostatic and non-hydrostatic models. It is believed by t if finer grid resolutions with larger aspect ratios had remained stable (i.e. the 81?9 and 81?27 cases), these cases would have shown effects from the inclusion of non-hydrostatic pressure. The mechanism causing the model instability for fine horizontal grid resolutions is presently unknown. The timestep is not controlling vertical motions, unlike the laboratory case (?4.1.1), so this is not a simple CFL issue. There are several possible sources of error for this instability. There may be an error in the non-hydrostatic model that was only detected in this situation. The most likely possibilit is the use of the explicit Euler discretization for the baroclinic term in CWR-ELCOM. While explicit Euler discretizations are unstable, they have been successfully applied as part of models, such as CWR-ELCOM where truncation error serves as stabilization. Stability is clearly only a problem w the non-hydrostatic pressure solution as the hydrostatic model remains stable under all tes conditions. The issue with ecrease in numerical dissipation seen in the laboratory-scale simulation. Another possibility fo this instability could be a more global problem with the application of the fractional-step method for internal wave evolution. While several other non- hydrostatic models have used the fractional-step method, there has been no reported work on grid refinement and very little work on internal wave modeling. Most of the known non-hydrostatic models (Mahadevan, et al. 1996a, b; Marshall, et a 1997; Stansby and Zhou, 1998; Casulli, 1999; Chen, 2003) do not explicitly model internal waves and present results only for laboratory scale problems. lake or ocean scale problems are investigated, the general circulation patterns and free-surface elevations are reported. Daily and Imberger (2003) model internal waves, however only on the laboratory scale. None of the aforementioned models present results for more than one lity at fine resolutions in a large domain remains an open issue for future research. 4.3 Summary The non-hydrostatic model laboratory-scale simulations in the present work compared well with the theory and laboratory experiments of Horn. Horn, et al. (2001) examined data from several published field studies and confirmed the applicability of the theoretical timescales used to develop a regime diagram. Therefore, it is inferred that the results of the laboratory-scale model simulations are rele The non-hydrostatic model provides a substantially better representation of internal wave evolution than the hydrostatic model for laboratory scales. The features of internal wave developme (such as soliton formation) are physical in the non- hydrostatic model, while a function of the grid resolution in the hydrostatic model. There is no appreciable difference between the hydrostatic and non-hydrostatic model for weakly nonlinear, subsequently weakly non-hydrostatic waves. In cases with a strongly nonlinear, non-hydrostatic wave on a sufficiently fine grid, such as in the laboratory simulations, the non-hydrostatic pres effect is significant. The addition of non-hydrosta pressure dispe s the non into a train of solitons, as observed in the laboratory, and captures the energy cascade through the wave train. Adding non-hydrostatic pressure does not eliminate numerical error, but it does change the characteristics of model error; numerical diffusion is the mechanism that determines m skill, whereas the hydrostatic dissipative. The largest grid aspect ratio applied to th scale case for a stable non-hydrostatic simulation is 0.012 (27?9 grid). The largest grid aspect ra applied to scenario 2 is 2.0 (600?15 grid), and smallest grid aspect ratio is 0.02 (30?73 grid). It deduced from the results of the laboratory scale and lake scale cases that a non-hydrostatic model needs a sufficiently large grid aspect ratio [?z/?x > O(10 -2 to affect internal wave evolution. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 54 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 t/T ? / ? m 9 x 9 nh 9 x 9 h 27 x 9 nh 27 x 9 h ? / ? m 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 t/T ? / ? m 9 x 9 nh 9 x 9 h 9 x 27 nh 9 x 27 h 9 x 81 nh 9 x 81 h ? / ? m ? / ? m 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 t/T ? / ? m 27 x 9 nh 27 x 9 h 27 x 27 nh 27 x 27 h 27 x 81 nh 27 x 81 h ? / ? m ? / ? m 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 t/T ? / ? m 9 x 81 nh 9 x 81 h 27 x 81 nh 27 81 h 81 h x 81 x ? / ? m ? / ? m 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 t/T ? / ? m 9 x 27 nh 9 x 27 h 27 27 nh 27 h x 27 x ? / ? m ? / ? m a) d) Nu a d 1, 7 l is eft b) c) e) Figure 4.20: for different gri 9, b) vertical gr d)horizontal gri Lines represent (nh), markers r results (h). The hydrostatic cas 81 makes the fi out. merical viscosity of lake-sc d resolutions. a)Vertical gri id = 27, c) vertical grid = 8 d = 9, e) horizontal grid = 2 non-hydrostatic model resu epresent hydrostatic model instability in all non- es where the horizontal grid gures unclear and so were l le = . ts Chapter 4. Comparison of models 55 10 0 10 1 10 2 10 -1 10 0 10 1 E B - E B (0 ) / E D (0 ) t/T 9 x 9 nh 9 x 9 h 27 x 9 nh 27 x 9 h 81 x 9 nh E B - E B (0 ) / E D (0 ) E B - E B (0 ) / E D (0 ) 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 E B - E B (0 ) / E D (0 ) t/T 9 x 27 nh 9 x 27 h 27 x 27 nh 27 x 27 h 81 x 27 nh E B - E B (0 ) / E D (0 ) E B - E B (0 ) / E D (0 ) 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 E B - E B (0 ) / E D (0 ) t/T 27 x 9 nh 27 x 9 h 27 x 27 nh 27 x 27 h 27 x 81 nh 27 x 81 h E B - E B (0 ) / E D (0 ) 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 E B - E B (0 ) / E D (0 ) t/T 9 x 9 nh 9 x 9 h 9 x 27 nh 9 x 27 h 9 x 81 nh 9 x 81 h E B - E B (0 ) / E D (0 ) E B - E B (0 ) / E D (0 ) a) b) d) e) Figure 4.21: Background potential energy of lake-scale for different grid resolutions. rtical grid = 9, b) vertical grid = 27, c) cal grid = 81, d)horizontal grid = 9, e) zontal grid = 27. Lines represent non- ostatic model results (nh), markers esent hydrostatic model results (h). a)Ve verti hori hydr repr 10 0 10 1 10 2 -1.5 -1 -0.5 0 0.5 1 1.5 2 E B - E B (0 ) / E D (0 ) t/T 9 x 81 nh 9 x 81 h 27 x 81 nh 27 x 81 h 81 x 81 nh 81 x 81 h E B - E B (0 ) / E D (0 ) E B - E B (0 ) / E D (0 ) c) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 56 10 0 10 1 10 2 -8 -6 -4 -2 0 2 4 6 8 x 10 4 ? / ? s 9 x 81 nh 9 x 81 h 27 x 81 nh 27 x 81 h 81 x 81 nh 81 x 81 h t/T ? / ? s ? / ? s 10 0 10 1 10 2 10 0 10 1 10 2 10 3 10 4 10 5 10 6 ? / ? s t/T 27 x 9 nh 27 x 9 h 27 x 27 nh 27 x 27 h 27 x 81 nh 27 x 81 h ? / ? s ? / ? s 10 0 10 1 10 2 10 0 10 1 10 2 10 3 10 4 10 5 10 6 ? / ? s t/T 9 x 27 nh 9 x 27 h 27 x 27 nh 27 x 27 h 81 x 27 nh ? / ? s ? / ? s 10 0 10 1 10 2 10 0 10 1 10 2 10 3 10 4 10 5 10 6 ? / ? s t/T 9 x 9 nh 9 x 9 h 9 x 27 nh 9 x 27 h 9 x 81 nh 9 x 81 h ? / ? s ? / ? s 10 0 10 1 10 2 10 0 10 1 10 2 10 3 10 4 10 5 10 6 ? / ? s t/T 9 x 9 nh 9 x 9 h 27 x 9 nh 27 x 9 h 81 x 9 nh ? / ? s ? / ? s a) b) c) e) d) Figure 4.22: Numerical diffusivity of lake scale for different grid resolutions. a)Vertical grid = 9, b) vertical grid = 27, c) vertical 81, d)horizontal grid = 9, e) horizontal gr 27. Lines represent non-hydrostatic mode results (nh), markers represent hydros model results (h). grid = id = l tatic Chapter 4. Comparison of models 57 ) vertical grid = 27, e) vertical grid = 81. Lines represent non-hydrostatic model results (nh), markers represent ults (h). Figure 4.23: Wavelengths with peak PSD of lake-scale for different grid resolutions. a)Horizontal grid = 9, b) horizontal grid = 27m c)vertical grid = 9, d hydrostatic model res 10 -2 10 -1 10 0 10 -1 10 0 10 1 10 2 t/ T ?/? o t/ T t/ T 9x81 nh 9x81 h 27x81 nh 27x81 h 81x81 h 10 -2 10 -1 10 0 10 -1 10 0 10 1 10 2 t/T ?/? o t/T t/T 9x27 nh 9x27 h 27x27 nh 27x27 h 10 -2 10 -1 10 0 10 -1 10 0 10 1 10 2 t / T ?/? o t / T t / T 9x9 nh 9x9 h 27x9 nh 27x9 h 10 -2 10 -1 10 0 10 -1 10 0 10 1 10 2 t/T ?/? o t/T t/T 27x9 nh 27x9 h 27x27 nh 27x27 h 27x81 nh 27x81 h 10 -2 10 -1 10 0 10 -1 10 0 10 1 10 2 t/T ?/? o t/T t/T 9x9 nh 9x9 h 9x27 nh 9x27 h 9x81 nh 9x81 h a) b) c) d) e) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 58 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 E B - E B (1 ) / E D (1 ) t/T 0.1 nh 0.1 h 1.0 nh 1.0 h 10.0 nh 10.0 h E B - E B (1 ) / E D (1 ) E B - E B (1 ) / E D (1 ) 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 E D / E D (1 ) t/T 0.1 nh 0.1 h 1.0 nh 1.0 h 10.0 nh 10.0 h E D / E D (1 ) E D / E D (1 ) 10 0 10 1 10 2 10 1 10 2 10 3 10 4 10 5 ? / ? s t/T 0.1 nh 0.1 h 1.0 nh 1.0 h 10.0 nh 10.0 h ? / ? s ? / ? s 0 10 20 30 40 50 60 70 80 90 100 -1000 -500 0 500 1000 1500 2000 2500 t/T 0.1 nh 0.1 h 1.0 nh 1.0 h 10.0 nh 10.0 h ? / ? m ? / ? m ? / ? m a) b) c) d) Figure 4.24: Varying the density gradient of the lake-scale simulation. a)Background potential energy, b) numerical diffusivity, c) dynamic energy, d) numerical viscosity. Lines represent non-hydrostatic model results (nh), markers represent hydrostatic model results (h). Chapter 5. Isolation of non-hydrostatic effects 59 effe In en oceans, est wave issue a non- neither nor pr press fou non- simu beha identif press advan hydr using could hydr simp me applicatio devel rema prese of fou wave (Fig The st beha non-hy is used to identify local regions with non-hydrostatic behavior. The new non-hydrostatic pressure isolation approach is based on the separation of hydrostatic parameter obtained from hydrostatic pressure isolation method terizes a cell as hydrostatic or non-hydrostatic. chapter provides the numerical theory used to p the non-hydrostatic pressure isolation me d an application of the method for two diff aves. Further, the non-hydrostatic pressure n method is applied to one wave on thre rent grid resolutions to assess the effect of grid ion on the non-hydrostatic pressure isolation . 5.1 Numerical Theory The governing equations (see ?2 eld a vertical momentum equation for invis : Chapter 5. Isolation of non-hydrostatic cts vironments with hydrostatic aspect ratios (e.g. uaries, lakes), non-hydrostatic pressure gradients are only locally significant, such as at steep fronts. Stansby and Zhou (1998) discussed this nd suggested isolating regions of significant hydrostatic pressure gradients; however, they identified regions of non-hydrostatic behavior ovided a local solution for the non-hydrostatic ure. The present research builds the ndations for linking hydrostatic models with local hydrostatic solutions to yield practical lations of non-hydrostatic internal wave vior. This chapter presents a new method for ying local regions where the non-hydrostatic ure is expected to be significant. The key ce in the state-of-the-art is that the non- ostatic region is demonstrated to be identifiable only a hydrostatic solution. Thus, this method be used to allow the application of a local non- hydrostatic solver where necessary, while ?hydrostatic? areas are more simply computed with a ostatic model. The ability to ?boot-strap? from ler hydrostatic methods to non-hydrostatic thods should allow faster solution than global n of non-hydrostatic methods. The opment of a local non-hydrostatic solver ins an issue for further research (?6 ff.), so the nt work can be seen as a necessary development ndational methods. As a wave propagates, the fluid around the front is pushed away from the wave front ure 5.1); this is where the strongest vertical velocities are located (Daily and Imberger, 2003). rong vertical velocities cause the wave to ve more nonlinearly, and subsequently, more drostatically. Thus, the vertical acceleration hydrostatic and non-hydrostatic pressure solutions in the fractional-step method, so that the regions with non-hydrostatic behavior are a priori identified from the hydrostatic solution. The hydrostatic solution provides the vertical acceleration used in the non- hydrostatic pressure isolation method. A non- the non- charac This develo thod an erent w isolatio e diffe resolut method .1.1) yi cid flow nh o PDW 1 Dt z ? =? ? ? (5.1) This equation is a relationship betwee al acceleration and the non-hydrostatic p vertical gradient, which ordinarily exists only on- hydrostatic model. The vertical veloc is representative of a flow field that accounts for non- hydrostatic pressure. In a hydrostatic the non-hydrostatic pressure is neglected uation (5.1) is not modeled. The hydrostatic l velocity is computed diagnostically b ing the continuity equation (see ?2.1.1). The hydrostatic pressure isolation method e hydrostatic vertical velocity and the relationship in Equation (5.1) to approximate the non static pressure. n vertic ressure for a n ity (W) model, and Eq vertica y enforc non- uses th -hydro * nh o PDW 1 Dt z ? =? ? ? (5.2) W* is the hydrostatic vertical velocity; this is the same vertical velocity used in Equation (2.7). Therefore, the left-hand side of Equat ) is known after the hydrostatic model is e right-hand side is unknown and is app with the non-hydrostatic pressure isolation comparison of the hydrostatic vertical acceleration (DW * /Dt) and the modeled non-hydrostatic pressure vertical gradient is seen in Figure 5.2. The hydrostatic vertical acceleration is the approximated ion (5.2 applied; th roximated method. A Wave Direction Figure 5.1: Direction of vertical velocities for a horizontally propagating wave front. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 60 , c e ion ak dients. The approximated and modeled non-hydrostatic pressure vertical gradients are symmetrical around the wave front, although the modeled non- hydrostatic pressure vertical gradient shows a slightly wider range of influence around the steep wave front. Aside from the principle wave front, there is some steepening associated with a higher mode wave at the right end of the basin in Figure 5.2. Both the approximated and modeled non-hydrostatic pressure vertical gradients identify this location. . non-hydrostatic pressure vertical gradient. Note ?modeled? refers to anywhere the non-hydrostati model is used and ?approximated? refers to when th hydrostatic model is used to make an approximat of the non-hydrostatic pressure. The location of pe values for the approximated non-hydrostatic pressure vertical gradient is a reasonable match for the peak modeled non-hydrostatic pressure vertical gra Horizontal flow dominates internal wave propagation and is prognostically determined in hydrostatic models, thus it is practical to use the horizontal momentum equation in the development of the non-hydrostatic pressure isolation method The non-hydrostatic horizontal momentum equation contains the non-hydrostatic pressure horizontal gradient:  nh x ? ? oo z' PDU 1 1 gdz Dt x x ? ?? ??? ? ? =? + ? ? ?? ?? ? ut an f the n-h accel h ? ? (5.3) The non-hydrostatic pressure horizontal gradient is not known in a hydrostatic model, b approximation of the non-hydrostatic pressure horizontal gradient effect is required to identify regions of significant non-hydrostatic behavior. The approximation for the non-hydrostatic vertical gradient can be used to estimate the non-hydrostatic pressure horizontal gradient. Vertically integrating Equation (5.2) produces an approximation o non-hydrostatic pressure: Figure 5.2: Comparison of hydrostatic/no Density profile; b) hydrostatic vertical vertical gradient. The main non-hydrostatic also some non-hydrostatic effect on the rig ydrostatic vertical momentum terms. a) eration; c) modeled non-hydrostatic pressure influence is at the steep wave front. There is t side of b and c. 1 2 3 4 5 6 7 8 9 10 2 4 6 Length (m) He i g h t ( m ) a) Density b) DW * /Dt c) dP nh /dz He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) Chapter 5. Isolation of non-hydrostatic effects 61 a) Density b) Approximate dP nh /dx 1 2 3 4 5 6 7 8 9 10 2 4 6 c) Model dP nh /dx He i g h t ( m ) * nh o DW Pdz Dt =?? ? (5.4) The horizontal spatial derivative of the approximated non-hydrostatic pressure [Equation (5.4)] is: * nh o PDW 1 dz xDt x ?? =? ?? ? (5.5) Equation (5.5) approximates the non-hydrostatic pressure horizontal gradient, as seen in Figure 5.3. Similar to the non-hydrostatic pressure vertical gradient, the approximated non-hydrostatic pressure horizontal gradient matches the occurrence of the modeled non-hydrostatic pressure horizontal gradien t at the steep wa e tatic effect is large. The screening parameter relates the non-hydrostatic pressure orizontal gradient to other terms within th orizontal momentum equation that influe internal wave evolution. Where the screening parameter is large, that location is considered to have signi the Fi pr t; c) model non- hydrostatic pressure horizontal gradient. ve front. A screening parameter is developed to quantify the relevance of the approximated non-hydrostatic pressure with regards to its influence on internal wave evolution, thereby identifying regions wher the non-hydros h e h nce ficant non-hydrostatic behavior. Non-hydrostatic pressure effects are linked to the nonlinearity of a wave and are propagated at speed of the wave (?6 ff.). As the non-hydrostatic pressure gradients are strongest at the wave front, their relative importance on internal wave evolution is dependent on how the non-hydrostatic pressure gradient scales with the horizontal acceleration. Therefore, the screening parameter that isolates gure 5.3: Comparison of non-hydrostatic pressure horizontal gradient. a) Density ofile; b) approximate non-hydrostatic pressure horizontal gradien Length (m) b) Approxi ate dP nh /dx 1 2 3 4 5 6 7 8 9 10 c) odel dP nh /dx He i g h t ( m ) Length (m) b) Approximate dP nh /dx 1 2 3 4 5 6 7 8 9 10 2 4 6 c) Model dP nh /dx He i g h t ( m ) Length (m) b) Approxi ate dP nh /dx 1 2 3 4 5 6 7 8 9 10 c) odel dP nh /dx He i g h t ( m ) Length (m) b) Approximate dP nh /dx 1 2 3 4 5 6 7 8 9 10 2 4 6 c) Model dP nh /dx He i g h t ( m ) Length (m) b) Approxi ate dP nh /dx 1 2 3 4 5 6 7 8 9 10 c) odel dP nh /dx He i g h t ( m ) Length (m) b) Approximate dP nh /dx 1 2 3 4 5 6 7 8 9 10 2 4 6 c) Model dP nh /dx He i g h t ( m ) Length (m) b) Approxi ate dP nh /dx 1 2 3 4 5 6 7 8 9 10 c) Model dP nh /dx He i g h t ( m ) Length (m) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 62 regions of non-hydrostatic behavior compares the non-hydrostatic horizontal gradient to the hydrostatic horizontal acceleration (DU * /Dt): nh o ** max,n x DU DU Dt Dt ? ?? ?? ?? ?? ?? ??? ?? (5.6) The difference between the maximum horizontal acceleration at a time step (n) and the horizontal acceleration at each cell is used to tigate the effect of cells with very small horizontal accelerations that may amplify the screening eter erroneously. The total horizontal acceleration is used to include the horizontal wave agation and the nonlinear steepening. The total horizontal acceleration can be decomposed into the horizontal wave propagation, which is characterized ocal horizontal acceleration: P1??? mi param prop by the l * U t ? ? (5.7) nonlinear steepening, which is characterized dvection acceleration: and the by the a ** ** UU UW xz?? s of the approximated non-hydrostatic dient, the screening parameter is: ?? + (5.8) In term pressure gra * ** max,n Dt Dt ??? ?? DW dz xDt DU DU ??? ?? ? ?? ?? ?? ? (5.9) Figure 5.4 shows the comparison between the approximated and modeled screening parameter. The approximated screening parameter slightly over- predicts non-hydrostatic regions through the water column, whereas the modeled screening parameter is located only at the steep wave front. 5.2 Application of the Non-hydrostatic Pressure Isolation Method The non-hydrostatic pressure isolation method is applied to two different internal waves to demonstrate the method?s ability to isolate regions with significant non-hydrostatic pressure effects. 5.2.1 Setup The test basin is a 2D rectangular basin with an initial internal wave setup as specified in Table 5.1. The two different internal waves (internal wave 1 and internal wave 2) were initialized as a cosine wave with a three-layer hyperbolic-tangent density profile. The amplitude of internal wave 1 is 1.125 m, yielding an amplitude ratio (a/h) of 0.5. The amplitude of internal wave 2 is 0.3 m, yielding a/h = 0.1. Internal wave 1 is steep (a/L = 0.11), while internal wave 2 is much less steep (a/L = 0.04). Internal wave 1 was simulated with three different square grids (20?15, 40?30 and 80?60, where each grid cell has dimensions of 0.5 m ? 0.5 125 m, sen to avoid damping the non-hydrostatic effect over a grid with a small aspect ratio (?4.3). Internal wave 2 was simulated ternal wav 7 2. 0 ess h py Density gradient ?? m, 0.25 m ? 0.25 m and 0.125 m ? 0. respectively). Square grids were cho with the 40?30 grid. 5.2.2 Results The non-hydrostatic pressure isolation method is an effective way to screen for and target areas that require a non-hydrostatic solution. However, the s of test basin. Value Units Table 5.1: Dimension In Length L Height H Upper layer depth h Depth ratio h/H Pycnocline thickn Characteristic Symbol e 1 Internal wave 2 10 10 m .5 7.5 m 25 3.0 m .3 0.4 -- 55 44kg/m 3 Chapter 5. Isolation of non-hydrostatic effects 63 Figure 5.4: Comparison of screening parameter, ?. a) Density profile; b,d) ? calculated with the approximate non-hydrostatic pressure horizontal gradient; c,e) ? calculated with the model non- hydrostatic pressure horizontal gradient. In b and c, as cells scale from dark blue to dark red, the cell is considered more non-hydrostatic. 1 2 3 4 5 6 7 8 9 10 2 4 6 Length (m) He i g h t ( m ) d) Approximate ? e) Model ? He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) i He i g h t ( m ) t t He i g h t ( m ) l ( ) He i g h t ( m ) ? 2 4 6 He i g h t ( m ) a) Density c) Model ? Length (m) He i g h t ( m ) 1 2 3 4 5 6 7 8 9 10 0 2 4 6 b) Approximate ? He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 64 reening parameter requires a criterion for scriminating hydrostatic and non-hydrostatic gions. Where the screening parameter exceeds the riterion, that cell is considered non-hydrostatic, hile a cell that has a screening parameter below the riterion is considered hydrostatic. Quantifying the riterion as a function of density gradients, wave haracteristics and grid resolution is beyond the ope of this work. However, to demonstrate the avior of the non-hydrostatic pressure isolation method, three different criteria were chosen (? = 1, 5 and 10) and applied a posteriori to model results. Figures 5.5 and 5.6 show that as the non-hydrostatic screening parameter is increased, the region that is considered non-hydrostatic decreases. Figure 5.6 shows spikes in the number of cells that exceed the criterion about every 400 timesteps. This time coincides with when the steep wave fronts reflect off the basin wall, where the wave has large vertical velocities and therefore the non-hydrostatic effect is strong. Three grid resolutions were applied to internal wave 1. The finest grid resolution (80?60) identified the most cells as non-hydrostatic, while the coarsest grid (20?15) identified the least (Figure 5.7). This does not imply that there are more non-hydrostatic regions with the finer grid, but rather there are more cells representing this region. Figures 5.4, 5.8 and 5.9 show that all grids have similar ability to identify the region of non-hydrostatic behavior. Internal wave 1 is strongly nonlinear and so the non-hydrostatic pressure effect is strong at the wave front. Internal wave 2 is weakly nonlinear, nonetheless the non-hydrostatic pressure isolation method identifies the steepest part of the wave, as seen in Figure 5.10. Figure 5.11 shows a comparison between the number of cells with a screening parameter greater than one for internal wave 1 and internal wave 2. Internal wave 2 shows large rises in the number of cells considered non-hydrostatic when the wave is traveling along the basin and has the least number of non-hydrostatic cells when the wave is at the basin walls. When the wave is reflecting off the basin walls, the vertical and horizontal accelerations are greater than anywhere else in the basin, so the screening parameter is able to properly detect these regions with non-hydrostatic behavior. When internal wave 2 is traveling along the basin, it should not behave non-hydrostatically, as it is a near- linear wave. However, the horizontal accelerations are very small so the screening parameter is amplified. The filter used in the screening parameter to eliminate small accelerations does not work in the limit of very small horizontal accelerations. Thus, the screening parameter delineates regions as non- hydrostatic erroneously. 5.3 Summary The new non-hydrostatic pressure isolation method identifies regions where the non-hydrostatic pressure effect is significant. For a sufficiently steep wave, the method is able to predict these regions quite well. For weakly nonlinear waves, the non-hydrostatic pressure isolation method tends to over-predict non- hydrostatic regions, due to the nature of a weakly nonlinear system. The method of the screening parameter takes the difference between the local horizontal acceleration and the maximum horizontal acceleration in the field to ameliorate the problem of small accelerations. Despite this construction, a weakly nonlinear system has such small accelerations that the non-hydrostatic screening parameter is still amplified. Thus, the present method is limited by the need to identify a characteristic horizontal acceleration scale, which herein is taken as the maximum value of the domain. If a wave is weakly nonlinear, and therefore hydrostatic, the non-hydrostatic pressure isolation method is not applicable. However, if a wave has significant steepening, then the non-hydrostatic pressure isolation method identifies non-hydrostatic regions. The results from the non-hydrostatic pressure isolation method must not be viewed alone, but with a general understanding of the nature of the modeled wave. sc di re c w c c c sc beh Chapter 5. Isolation of non-hydrostatic effects 65 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 800 1000 # Timesteps # Ce l l s > ? ? = 1 ? = 5 ? = 10 Figure 5.6: Number of cells considered non-hydrostatic for internal wave 1; ? > criteria for a 40?30 grid -hydrostatic and Figure 5.5: Internal wave 1 domains with dif white are regions considered non ferent non-hydrostatic screening parameters. The the black are regions considered hydrostatic. Leng De p t h (m th (m) ) 1 2 3 4 5 6 7 8 9 10 0 2 4 6 ? = 5 ? = 10 = 5 ? = 1 De p t h (m ) Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 66 0 200 400 600 800 1000 1200 1400 1600 0 1000 2000 3000 4000 # timesteps # c e l l s ? > 1 40x30 20x15 80x60 # c e l l s ? > 1 # c e l l s ? > 1 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 800 1000 # Timesteps # C e l l s ? > 1 0 40x30 20x15 80x60 # C e l l s ? > 1 0 # C e l l s ? > 1 0 0 200 400 600 800 1000 1200 1400 1600 0 500 1000 1500 # timesteps # C e lls ? > 5 40x30 20x20 80x60 # C e lls ? > 5 # C e lls ? > 5 # C e lls ? > 5 # C e lls ? > 5 # C e lls ? > 5 # C e lls ? > 5 Figure 5.7: Number of cells considered non-hydrostatic for internal wave 1. Each panel has three different grids (blue line is 40?30 grid, red line is 20?15 grid and green line is 80?60 grid). The top panel is for a screening parameter, ? > 1; the middle panel has ? > 5 and the bottom panel has ? > 10. Chapter 5. Isolation of non-hydrostatic effects Figure 5.8: Comparison of non-hydrostatic screening parameter, gamma, on the 20?15 grid for internal wave 1. a) Density profile; b) ? calculated with the approximate non-hydrostatic pressure horizontal gradient; c) ? calculated with the model non-hydrostatic pressure horizontal gradient. 1 2 3 4 5 6 7 8 9 10 2 4 6 Length (m) He i g h t (m) b) Approximate ? e) Model ? a) Density He i g h t (m) 67 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 68 th 1 2 3 4 5 6 7 8 9 10 2 4 6 Length (m) He i g h t (m) c) Model ? b) Approximate ? a) Density He i g h t (m) He i g h t (m) Figure 5.9: Comparison of non-hydrostatic screening parameter, gamma, on the 80?60 grid for internal wave 1. a) Density profile; b) ? calculated with the approximate non-hydrostatic pressure horizontal gradient; c) ? calculated with e model non-hydrostatic pressure horizontal gradient. Chapter 5. Isolation of non-hydros tatic effects 69 1 2 3 4 5 6 7 8 9 10 2 4 6 Length (m) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) He i g h t ( m ) a) Density b) Approximate dP nh /dx c) Model dP nh /dx Figure 5.10: Comparison of non-hydrostatic pressure horizontal gradient for internal wave 2. a) Density profile; b) approximate non-hydrostatic pressure horizontal gradient; c) model non-hydrostatic pressure horizontal Figure 5.11: Number of cells considered non-hydrostatic for internal wave 1 (blue line) and internal wave 2 (red line). ? > 1. gradient. 0 0.5 1 1.5 2 2.5 3 3.5 4 0 500 1000 1500 t/T # C e l l s ? > 1 # C e l l s ? > 1 # C e l l s ? > 1 # C e l l s ? > 1 internal wave 1 internal wave 2 Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 70 Chapter 6. Conclusions and Recommendations 71 Chapter 6. Conclusions and Recommendations The two main objectives of the present research are: ? quantify the differences between hydrostatic and non-hydrostatic simulations of internal wave evolution, and ? develop a method to a priori determine regions with non-hydrostatic behavior. Both of these objectives have been accomplished, and the results provide insight into the mechanisms that control internal wave modeling. This chapter: 1) summarizes the work and conclusions, 2) discusses a numerical issue that limits the rapidity of the pressure Poisson solution?s convergence and 3) provides recommendations for future work. 6.1 Summary Discussion Differences between the hydrostatic and non- hydrostatic representation of internal waves are considerable (?4.1). The hydrostatic model simulates half of the physics necessary to capture internal wave evolution (i.e. nonlinear acceleration). This allows the wave to nonlinearly steepen, but there is no physical term to model dispersion. Thus, the hydrostatic model is limited in capability to modeling hydrostatic shallow waves that remain essentially linear damped waves (?4.1, scenarios 1 and 6). However, the hydrostatic model is incapable of physically capturing the dispersion of a steepened wave into a train of solitons, as the hydrostatic model neglects the physics that control this process. The soliton-like formations observed in the hydrostatic model (in this work (?4) and Hodges and Delavan, 2004) change when the grid resolution for a hydrostatically-modeled wave is varied. This indicates that the size and shape of the soliton-like formations are dependent on the grid (?4.1.3) and therefore are a model fabrication. This finding is further confirmed by the power spectral density distributed over wavelengths; the peak power is not shifted to the wavelengths of the soliton-like formations, and any power that is seen in these waves is a result of harmonics (?4.1.3). Numerical error within the hydrostatic model accumulates to reduce the energy in the wave and eventually damp out the wave?s steepness. Numerical diffusion can be reduced by vertically refining the grid, but the hydrostatic system is numerically dissipative, irrespective of the grid resolution. Unlike the hydrostatic model, the non- hydrostatic model includes the non-hydrostatic pressure term, which acts as a dispersive force. Thus, as a wave steepens nonlinearly, the non- hydrostatic pressure applies an effect that disperses the wave into a train of solitons. Grid resolution affects the occurrence of the train of solitons (i.e. finer horizontal grids simulate more emerging soliton trains than a coarser grid, Figure 4.8). However, the characteristic wavelength of the solitons is relatively unaffected by grid resolution, compared to the hydrostatic model. The peak power in a wave system shifts to the characteristic soliton wavelength after solitons have emerged; the shift in power is indicative of an energy transfer from basin-scale wavelengths to soliton wavelengths. It is thus concluded that the non-hydrostatic model is indeed modeling the physics of internal wave evolution. Non-hydrostatic models do not eliminate numerical error, however, the characteristics of numerical error are different in a non-hydrostatic model than in a hydrostatic model. Like the hydrostatic model, vertical grid refinement reduces numerical diffusion. In contrast to the hydrostatic model, horizontal grid refinement greatly reduces numerical dissipation in the laboratory-scale case. The laboratory-scale case (?4.1.3) decreases dynamic energy by only 8% over nine wave periods for the finest grid resolution compared to the coarsest grid resolution which decreases dynamic energy by 85%. It appears that the neglect of non-hydrostatic pressure effectively acts as a dissipative term in the hydrostatic equations. Perhaps an alternative viewpoint is that the non-hydrostatic pressure redistribution of momentum reduces the numerical dissipation that is associated with shallow water conservation of momentum. Grid resolution and the subsequent grid aspect ratio are the two model parameters that control the non-hydrostatic model?s ability to capture internal wave evolution. Most of the lake-scale cases had very small grid aspect ratios [O(10 -3 ) ? O(10 -4 )], and showed no significant difference between the hydrostatic and non-hydrostatic models. The one lake-scale case non-hydrostatic model that differed (in diffusion only) from the hydrostatic model had a grid aspect ratio of O(10 -2 ); the laboratory scale case, where results were significantly different between the hydrostatic and non-hydrostatic models, had grid aspect ratios of at least O(10 -2 ). Due to the successful use of the non-hydrostatic model for larger grid aspect ratios, the present research suggests that the grid aspect ratio must be at least of O(10 -2 ) to include the effect of non-hydrostatic Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 72 ial resolution error non- poral resolution must also condition. he and is pr gradients ve. The anal ca son equation nce ves d For example, flow .g. Su, et al., 1999) e step tatic the me e id needed to effectively apply the non-hydrostatic model renders the s ture w c ler e n- th n he a n- - er nt r -2 ) pressure. In addition to the fine spat needed for a physical, low numerical hydrostatic model, the tem be fine enough to maintain a small vertical CFL T pressure Poisson solution is inherently stiff oblematic for the motion of stiff pressure at the front of an internal wa yti l solution (?3.1) provides a pressure Pois that is independent of a flow field. However, the pressure Poisson solution for hydrodynamic problems (e.g. internal wave evolution) depends upon the gradient of velocity. The gradient of velocity generally does not match the speed of propagation of the non-hydrostatic pressure gradient. Indeed, this appears to be a key differe between non-hydrostatic solutions of internal wa and flows that have been more commonly studie using non-hydrostatic models. over a backwards facing step (e has an unsteady recirculation region behind th due to non-hydrostatic pressure. The non-hydros pressure gradients in the recirculation oscillate at same timescale as the motion of the reattachment point, which is much longer than the timescale of advection. Thus, the advection velocity limits the model time step, from which it follows that local changes in the non-hydrostatic pressure gradient from one time step to another are small. As the ti ?n+1? pressure field is a good initial guess at the tim ?n? pressure field in a backwards-facing step, a Poisson solver should be able to rapidly converge. In contrast, when modeling an internal wave, for common subcritical flows the wave is faster than the fluid velocity, so the baroclinic wave speed limits the model time step. As the internal wave propagates rapidly through the model grid, non-hydrostatic pressure gradients at the wave front are propagated at the wave speed. Therefore, the non-hydrostatic pressure field may change quite dramatically between time steps, such that the time ?n? pressure field is a poor approximation of the time ?n+1? pressure field. This effect is demonstrated by simulating a monochromatic wave (?3.2) under two different conditions for the Poisson solution: 1) with the time ?n? pressure is used as the starting point for the time ?n+1? pressure solution; and 2) with zero pressure used as the starting point for all pressure solutions. Figure 3.9 shows that using the time ?n? pressure always requires a larger number of iterations to meet the convergence criterion. The fine timestep and spatial gr olution of real-world scale problems a fu capability. The present research investigated a ne method to identify regions with significant non- hydrostatic behavior. Using the non-hydrostati pressure isolation method, a model may be developed which would ?boot-strap? from a simp hydrostatic model to a non-hydrostatic model in th regions identified as ?non-hydrostatic.? Such a hydrostatic/non-hydrostatic hybrid model should allow faster solution of internal wave evolution. However, there is a limit to the ability of the no hydrostatic pressure isolation method. The method over-predicts non-hydrostatic regions in shallow waves. These shallow waves are ably modeled wi the hydrostatic model, since steepening and solito formation is not the dominant phenomenon in t wave?s evolution. Thus, the non-hydrostatic pressure isolation method should only be applied to waves which are expected to steepen and form solitons (i.e. regime 2, Figure 3.10). It is crucial to water quality modeling to have hydrodynamic model to skillfully simulate the internal wave evolution and the flow field in a basin (Gross, et al, 1999). For damped linear internal waves, a hydrostatic model may adequately model the behavior as long as viscous damping is modeled. For internal waves where nonlinear steepening and soliton formation are the prevalent features, a no hydrostatic model provides a physical simulation of internal wave evolution on a sufficiently fine grid, whereas the hydrostatic model does not. Thus, non hydrostatic models should be used when nonlinear internal waves and solitons are present in the flow field; this will in turn create a more realistic representation of all the system processes that depend on the flow field, such as those in a wat quality model. 6.2 Conclusions The conclusions that may be drawn from the prese work are: ? The present hydrostatic model of the Eule equations shows soliton-like formations for nonlinear/non-hydrostatic waves, which are an artifact of the numerical model and grid resolution. ? There is some evidence that the grid aspect ratio (?z/?x) must be greater than O(10 to capture the effect of non-hydrostatic pressure in the present non-hydrostatic model. Chapter 6. Conclusions and Recommendations 73 y mendations st be ulation n d- g eed, ce c hydr the gion to r e on- e handled by buffe to t e c o finer grids. Additionally, the essure solver should be conditioned conjugate gradient Ack o This rese Naval R 2001 of 01-1-057 tre for Water Research at the Uni ELC research t of Texas ? A hydrostatic model can be used to identif internal-wave regions with non-hydrostatic behavior. 6.3 Recom The research presented here investigates emerging areas within non-hydrostatic modeling. This work began to chip away at some of the issues challenging non-hydrostatic modeling, but much remains to be investigated. The most notable issue that mu examined in non-hydrostatic modeling is sim time. To diminish simulation times, this author believes there are two main approaches: 1) address the issue of stiffness, and 2) construct a hybrid model that uses the non-hydrostatic pressure isolation method. The pressure Poisson equation is inherently stiff, but the stiffness is exacerbated by the advectio methods used. Presently, the velocity at the previous timestep is used in the source term of the pressure Poisson equation; this velocity moves much slower than the wave and the non-hydrostatic pressure gradients. One suggestion to improve the source term of the pressure Poisson equation would be to use a predicted advection velocity field from the previous two timesteps. For instance, use a secon order Adams-Bashforth which predicts the ?n+1? velocity from the ?n-1? and ?n? velocities. Using the past two timesteps may better estimate the strong velocity gradients near the propagating wave front. A second approach to improve the source term in the pressure Poisson equation is to use the non- hydrostatic pressure isolation method to identify the point at an instant in time with the peak non- hydrostatic pressure gradient. This point should coincide with the propagating wave front and the strongest velocity gradients. Identifying the location with peak non-hydrostatic pressure gradients, alon with the known timestep, will yield the velocity that of the peak non-hydrostatic pressure gradient. This velocity should be characteristic of the wave sp which may be used to provide a better estimation of the non-hydrostatic pressure field. The non-hydrostatic pressure isolation method developed here only provides the relative importan of the non-hydrostatic pressure effect in the internal wave?s evolution. In application, the non-hydrostati screening parameter needs some criterion to delineate a region as ?hydrostatic? or ?non- ostatic.? For example, where the screening parameter is above a preset criterion, a cell is designated as ?non-hydrostatic.? At this cell, the non-hydrostatic solver may be called to calculate non-hydrostatic pressure and the updated velocity field and free-surface. The criterion for a re be considered non-hydrostatic must be established. The author does not know if there will be a global value for the non-hydrostatic screening paramete criterion, or if it will be dependent on each wave characteristics. Another issue associated with the implementation of the non-hydrostatic pressure isolation method is that any change in the velocity field propagates instantaneously over the entire domain. Thus, changes within the interior of the domain may cause perturbations on the free-surfac that could feed back into the flow in the form of non- conservation in mass and volume. There must be some transition between the ?hydrostatic? and n hydrostatic? cells; the transition may b r cells. The buffer cells will account for any free-surface perturbation and make a correction the system to ensure conservation of mass and volume. The size of the buffer region must also be investigated to determine if it is a function of the magnitude of the non-hydrostatic pressure gradients or the size of the non-hydrostatic region. Finally, the non-hydrostatic model must be examined on more real-world scale cases with different grid resolutions to determine if the instability found in the lake-scale case (?4.2) is a model error or an issue with the fractional-step method. To verify the present research?s assessmen that a grid aspect ratio of at least O(10 -2 ) must b used to observe differences between the hydrostati and non-hydrostatic model, the lake-scale case should be applied t present non-hydrostatic pr adapted to use a pre method or a multigrid method to improve convergence time. n wledgements arch was supported under the Office of esearch Young Investigator Program Award Ben R. Hodges, ONR Contract No. N00014- 4. The Cen versity of Western Australia provided the OM numerical model used in this work under a partnership agreement with the Universi y at Austin. Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 74 Bibliography 75 Bibliography Apel, J.R., (1981). Satellite sensing of ocean surface dynamics. Ann. Rev, Earth & Planetary Sci., 8:303-342. Armfield, S. and R. Street, (1999). The fractional-step method for the Navier-Stokes equations: the accuracy of three variations. J. Comp. Phys., 153:2:660-665. Batchelor, G.K., (1967). An Introduction to Fluid Dynamics. Cambridge Univeristy Press, New York, 615 pgs. Blumberg, A.F., & G.L. Mellor, (1987). A description of a three-dimensional coastal ocean circulation model, in Three-Dimensional Coastal Ocean Models, N. S. Heaps (ed.), American Geophysical Union, Washington D.C. pp. 1-16. Boegman, L., J. Imberger, J.N. Ivey, & J.P. Antenucci, (2003). High frequency internal waves in large stratified lakes. Limnol. & Oceanogr., 48:12:895-919. Casulli, V., (1999). A semi-implicit finite difference method for non-hydrostatic, free-surface flows. Int. J. for Numer. Methods in Fluids, 30:425- 440. Casulli, V. & G.S. Stelling, (1998). Numerical simulation of 3D quasi-hydrostatic free-surface flows. J. Hydraul. Engng. 124:678-686. Chen, X., (2003). A fully hydrodynamic model for three-dimensional free-surface flows. Int. J. for Numer. Methods in Fluids, 42:929-952. Cushman-Roisin, B., (1994). Introduction to Geophysical Fluid Dynamics. Prentice Hall, New Jersey, 320 pgs. Daily, C. & J. Imberger, (2003). Modeling solitons under the hydrostatic and Boussinesq approximations. Int. J. for Numer. Methods in Fluids, 43:231-252 Delavan, S.K., (2003). Evolution Accumulation of Numerical Errors in Hydrostatic Models of Internal Waves, M.S. Thesis, Department of Civil Engineering, University of Texas, Austin, Aug, 2003. De Silva, I.P.D., J. Imberger, & G.N. Ivey, (1997). Localized mixing due to a breaking internal wave ray at a sloping bed. J. Fluid Mech. 350:1-27. Ezer, T., H, Arango, A.F. Shchepetkin (2002). Developments in terrain-following ocean models: intercomparisons of numerical aspects. Ocean Modelling, 4:249-267 Farmer, D.M., (1978). Observations of long nonlinear internal waves in a lake. J. Phys. Oceanogr., 8:63-73. Ferziger, J.H. & M. Peri?, (2002). Computational Methods for Fluid Mechanics. Springer, Berlin, Germany, 423 pgs. Fringer, O. B., & R. L. Street. 2003. The dynamics of breaking progressive interfacial waves. J. Fluid Mech. 494: 319-353. Garrett, C. and W. Munk, (1979). Internal waves in the ocean. Ann. Review Fluid Mech., 11:339-369. Gill, A.E., (1982). Atmosphere-Ocean Dynamics, Academic Press: New York, 1982. Gross, E.S., J.R. Koseff, & S.G. Monismith, (1999). Three-dimensional salinity simulations of South San Francisco Bay, J. Hyd. Engr., 125:11:1199-1209. Haidvogel, D.B., H.G. Arango, K. Hedstrom, A. Beckmann, P. Malanotte-Rizzoli, A.F. Shchepetkin, (2000). Model evaluation experiments in the North Atlantic Basin: simulations in nonlinear terrain-following coordinates. Dynamics of Atmospheres and Oceans. 32:239?281. Hamrick, J.M., (1992). A three-dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. The College of William and Mary, Virginia Institute of Marine Science, Special Report, 317, 63 pgs. He, P., M. Salcudean, I.S. Gartshore, and P. Nowak, (1996). Multigrid calculation of fluid flows in complex 3D geometries using curvilinear grids. Computers and Fluids, 25:4:395-419. Hirsch, C., (1988). Numerical Computation of Internal and External Flows, Volume 1: Fundamentals of Numerical Discretization. John Wiley & Sons, Inc., New York, 515 pgs. Hodges, B.R., (2000). Numerical Techniques in CWR-ELCOM (code release v.1). CWR Manuscript WP 1422 BH, Centre Water Res., Nedlands, Western Australia, Australia. Hodges, B.R., (2004). Accuracy of semi- implicit shallow-water Crank-Nicolson discretization. ASCE J. Eng. Mech., 130:8:904-910. Hodges, B.R. & S.K. Delavan, (2004). Numerical diffusion and dissipation in hydrostatic Wadzuk & Hodges: Hydrostatic and Nonhydrostatic Internal Wave Models 76 waves. ASCE Engineering ence Proceedings, June 13-16, 2004 00). Modeling basin-scale internal wave ge-scale interfacial gravity wave 2). A weakly nonlinear model of long h., Imbe d.], paradigm of planetary problems. Elsev 77-491. ion 396:183- 201. c Process, Acad Geophysical Fluid Flows, Acad ). g and breaking over a slope. J. Flu d ress, New York, 730 pgs. ts with a ydraul. Engrgr., 129:215-224. Dover ges, and R. Stock s: ic rpolation. Computer Meth. in Appl. Mech cons 8), Waves in Fluids, Cam dal flow ng, . Street, (199 l. (199 Adcroft, (199 Z. (2001). Dynamics of Inter onto models of internal Mechanics Confer . Hodges, B.R., J. Imberger, A. Saggio and K. Winters, (20 s in a stratified lake, Limnology and Oceanography, 45:7:1603-1620. Horn, D.A., J. Imberger, & G.N. Ivey, (2001). The degeneration of lar s in lakes. J. Fluid Mech., 434:181-207. Horn, D. A., J. Imberger, G. N. Ivey, and L. G. Redekopp. (200 internal waves in closed basins. J. Fluid Mec 467:269-287. Hutter, K., G. Bauer, Y. Wang & P. Guting, (1998). Forced motion response in enclosed lakes. In Physical Processes in Lakes and Oceans (ed. J. rger), AGU. Imberger, J., (1994). Transport processes in lakes: A review, p. 79-193. In R. Margalef [e Limnology now: A ier Science. Imberger, J. & G.N. Ivey, (1993). Boundary mixing in stratified reservoirs. J. Fluid Mech., 248:4 Javam A., J. Imberger & S.W. Armfield, (1999). Numerical study of internal wave reflect from sloping boundaries. J. Fluid Mech., Kantha, L.H. & C.A. Clayson, (2000a). Numerical Models of Oceans and Oceani emic Press, New York, 940 pgs. Kantha, L.H. & C.A. Clayson, (2000b). Small Scale Processes in emic Press, New York, 888 pgs. Kao, T.W., F. Pan, & D. Renouard, (1985 Internal solitons on the pycnocline: generation, propagation, and shoalin id Mech., 159:19-53. Kim, J. and P. Moin, (1985) Application of a fractional-step method to incompressibleNavier- Stokes equations. J. Comp. Phys., 59:308-323. Kinsman, B., (1965). Wind Waves: Their Generation and Propagation on the Ocean Surface. Dover Publications, Inc., New York, 676 pgs. Kreyszig, E, (1999). Advanced Engineering Mathematics. John Wiley & Sons, Inc., New York, 1156 pgs. Kundu, P.K. and I.M. Cohen, (2002). Flui Mechanics, Academic P Laval, B., B.R. Hodges, & J. Imberger, (2003a). Reducing numerical diffusion effec pycnocline filter. J. H Lamb, H., (1932). Hydrodynamics. Publications, Inc., New York, 738 pgs. Laval, B., J. Imberger, B.R. Hod er. (2003b). Modeling Circulation in Lake Spatial and Temporal Variations. Limnology and Oceanography 48:3:983-994. Leonard, B.P. (1979). Stable and accurate convective modeling procedure based on quadrat upstream inte . & Engr., 19:1:59-98. Leonard, B.P. (1991). The ultimate ervative difference scheme applied to unsteady one-dimensional advection. Computer Meth. in Appl. Mech. & Engr., 88:17-74. Lighthill, J. (197 bridge University Press, Cambridge, 504 pgs. Lin, B. and R.A. Falconer, (1997). Ti and transport modeling using ULTIMATE QUICKEST scheme. J. Hyd. Engr., 123:4:303-314. Long, R.R., (1972). The steepening of lo internal waves. Tellus, 24:88-99. Lorenz, E.N., (1955). Available potential energy and the maintenance of the general circulation. Tellus, 7:157-167. Mahadevan, A., J. Oliger, & R. L 6a). A non-hydrostatic mesoscale ocean mode Part I: Wellposedness and scaling. J. Phys. Oceanogr., 26:1868-1880. Mahadevan, A., J. Oliger, & R. L. Street, 6b). A non-hydrostatic mesoscale ocean model. Part II: Numerical implementation. J. Phys. Oceanogr., 26:1881-1900. Marshall, J., C. Hill, L. Perlman, & A. 7). Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling, J. Geophys. Res., 102:5733-5752. Miropol?sky, Yu. nal Gravity Waves in the Ocean, Kluwer, Dordrecht, 406 pgs. Moore, G. (1965). Cramming more components integrated circuits. Electronics, 38:8. Bibliography 77 and the respective biogeochemical processes in th . Internal nce, 208:4443:451- 460. ng. 124:1 autics, 41:4:595-604. sing Toolbox User?s Guide, (199 lems. Int. J. Numer. Methods in Fluid ). simulation by the gas-kinetic schem ). A coastal ocean model intercomparison study for a fluids. nonli - ng -18, 2003. astal Modeling: Proc rence, ar . 86). The lake. y & E.A. mixi ith nt. J. 5. Nishri, A., J. Imberger, W. Eckert, I. Ostrovsky, & Y. Geifman, (2000). The physical regime e lower water mass of Lake Kinneret, Limnol. Oceanogr., 45:4:972-981. Osborne A.R. and T.L. Burch, (1980) solitons in the Andaman sea. Scie Roache, P. J., (2002). Code verification by the method of manufactured solutions. J. Fluids E :4-10. Roy, C.J., (2003). Grid convergence error analysis for mixed-order numerical schemes. Amer. Inst. Aeronautics Astron Segur, H. & J.L. Hammack, (1982). Soliton models of long internal waves. J. Fluid Mech., 118:285-304. Signal Proces 8). The MathWorks, Inc. Stansby, P.K. & J.G. Zhou, (1998). Shallow- water flow solver with non-hydrostatic pressure: 2D vertical plane prob s, 28:542-563. Su, M. D., K. Xu, and M. S. Ghidaoui. (1999 Low-speed flow e. J. Comput. Phys., 150:1:17-39. Tartinville, B., E. Deleersnijder, P. Lazure, R. Proctor, K.G. Ruddick, and R.E. Uittenbogaard, (1998 three-dimensional idealized test case. App. Math. Modelling., 22:3:165-182. Turner, J. S. (1973). Buoyancy effects in Cambridge Univ. Press, NY USA. van Haren, H., (2004). Some observations of nearly modified internal wave spectra. J. Geophys. Res., 109(C03045). Wadzuk, B.M. & B.R. Hodges, (2003). Comparing hydrostatic and non-hydrostatic Navier Stokes models of internal waves. ASCE Engineeri Mechanics Conference, July 16 Weilbeer, H. and J.A. Jankowski, (2000). A three-dimensional non-hydrostatic model for free surface flows ? development, verification and limitations. Estuarine and Co eedings of the 6th International Confe Nov. 3-5, 1999. Whitham, G.B., (1974). Linear and nonline waves, Wiley, NY, USA, 636 pgs Wiegand, R.C. and E.C. Carmack, (19 climatology of internal waves in a deep temperate J. Geophys. Res., 91:C3:3951-3958. Winters, K.B., P.N. Lombard, J.J. Rile D?Asaro, (1995). Available potential energy and ng in density-stratified flows. J. Fluid Mech., 289:115-128. Yuan, H. & C.H. Wu, (2004). A two- dimensional vertical non-hydrostatic ? model w an implicit method for free-surface flows. I Numer. Meth. Fluids, 44:811-83