Wave induced hydrodynamic complexity and transport in the nearshore

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Request accessible transcript

Transcript (if available)

Content WAVE INDUCED HYDRODYNAMIC COMPLEXITY AND TRANSPORT IN THE NEARSHORE by Sangyoung Son A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL ENGINEERING) August 2012 Copyright 2012 Sangyoung Son Dedication To my lovely wife Bohyang Kim, and our beautiful children, Jaiden Jiho Son and Irene Jiin Son ii Acknowledgements I would like to express my deepest appreciation to Dr. Patrick Lynett for his guidance, support, patience and trust throughout my doctoral studies. It has truly been an honor and a pleasure to work with you for past four years. I also would like to thank my committee members, Dr. Jiin-Jen Lee, Dr. Larry Redekopp, Dr. Carter Wellford, Dr. Hung Leung Wong for their valuable comments on my works. I would not forget to thank Dr. Scott Socolofsky, Dr. James Kaihatu, Dr. Jennifer Irish and Dr. Kuang-An Chang for their concerns and advices during the PhD program at Texas A&M University. Indeed, I missed Aggie life a lot since I moved out. I am also sincerely grateful to Dr. Moo-Hyun Kim and ANMCCF family who always pray for our family. In LA, we used to recollect all the joyful moments shared with ANMCCF family. Most importantly, none of this work would have been possible without unconditional support from my family. I would like to thank my parents and parents-in-law for giving me love, care and encouragement during my study. I am also deeply appreciate my beloved children, Jaiden Jiho Son and Irene Jiin Son who always welcome me with a big smile at the door. Please forgive your dad, who was often busy with the study. I promise that I will play with you more in the future. And finally, I would like to thank Bohyang Kim for being my wife, my friend, my comfort, and my sanity. I love you so much. iii Table of Contents Dedication ii Acknowledgements iii List of Tables vii List of Figures viii Abstract xii Chapter 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives of Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Chapter 2 Nested and Multi-Physics Modeling of Long Waves 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Shallow Water Equation Model . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Governing Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Boussinesq Equation Model . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1 Governing Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.2 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4 Preliminary Discussion for Mismatches between Models . . . . . . . . . . 20 2.4.1 Mismatch in Physics . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.2 Mismatch in Numerics . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 Coupling Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5.1 Coupling Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5.2 Coupling with Different Grid Size. . . . . . . . . . . . . . . . . . . 29 2.5.3 Special Numerical Interface Treatment . . . . . . . . . . . . . . . . 30 2.6 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.6.1 Gaussian Hump Simulation . . . . . . . . . . . . . . . . . . . . . . 32 2.6.2 Physical and Numerical Setup. . . . . . . . . . . . . . . . . . . . . 33 2.6.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.7 Tsunami Wave Fission Simulation . . . . . . . . . . . . . . . . . . . . . . 42 iv 2.8 2004 Sumatra Tsunami Simulation . . . . . . . . . . . . . . . . . . . . . . 46 2.8.1 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.8.2 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Chapter 3 Interaction of Shallow Water Waves with Weakly Sheared Currents of Arbitrary Profile 58 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Brief Review on Wave-Current Interactions . . . . . . . . . . . . . . . . . 62 3.2.1 Waves over the Current . . . . . . . . . . . . . . . . . . . . . . . . 62 3.2.2 Currents under the Waves . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 Boussinesq Equations for Waves and Currents . . . . . . . . . . . . . . . . 65 3.3.1 Non-dimensionalized Governing Physics and Boundary Conditions 65 3.3.2 Reynolds Stresses Under Combined Wave-Current Flow . . . . . . 68 3.3.3 Subgrid Scale Turbulent Closure Model . . . . . . . . . . . . . . . 69 3.3.4 Depth-Integrated Momentum Equations for Waves and Currents . 70 3.3.5 Depth-Integrated Continuity Equation for Waves and Currents . . 78 3.3.6 Modulation of Dispersion Properties by Currents . . . . . . . . . . 79 3.3.7 Free parameter (b x ,b y ) . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.4.1 Waves over Uniform or Linearly-varying Currents . . . . . . . . . . 84 3.4.2 Waves and Turbulent Currents with Bed Roughness . . . . . . . . 90 3.4.3 Bichromatic Waves and Uniform Currents in a Vayring Depth. . . 95 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 Chapter 4 A Depth-Integrated Model for Free Surface Waves Propa- gatingoverFluidswithWeakVerticalandHorizontalDen- sity Variation 104 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2.1 Governing Physics and Boundary Conditions . . . . . . . . . . . . 109 4.2.2 Derivation of Boussinesq-type Equations for Dispersive Waves over Variable Density Fluid . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.3 Limiting Cases of Derived Model . . . . . . . . . . . . . . . . . . . . . . . 121 4.3.1 Linear, Inviscid Equations . . . . . . . . . . . . . . . . . . . . . . . 121 4.3.2 First-Order Nonlinear, Inviscid Equations . . . . . . . . . . . . . . 122 4.4 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.4.1 Horizontally Varying Fluid Density - Pneumatic Breakwater. . . . 123 4.4.2 Waves Excited by Internal Motion: Linear Surface Waves . . . . . 128 4.4.3 Waves Excited by Internal Motion: Nonlinear Surface Waves . . . 135 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 Chapter 5 Conclusions and Future Works 142 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 v Bibliography 146 Appendix A Numerical Scheme of COMCOT . . . . . . . . . . . . . . . . . . . . . . . . . . 157 Appendix B Second order terms in Boussinesq Equation . . . . . . . . . . . . . . . . . . . . 159 Appendix C Variables in Numerical Scheme of Boussinesq Model . . . . . . . . . . . . . . . 162 Appendix D Derivation of Momentum Equation in Wave-Current Model . . . . . . . . . . . 165 Appendix E Numerical Scheme in Wave-Current Model . . . . . . . . . . . . . . . . . . . . 170 E.1 Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 E.2 Time Marching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 Appendix F Derivation of Momentum Equation of Boussinesq Model for Variable Density Fluid Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 Appendix G Numerical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183 G.1 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 G.2 Spatial Discretization : Finite Volume Method . . . . . . . . . . . . . . . 189 Appendix H Second-order Sub- and Super-harmonic Solution in Two-layer Fluids . . . . . . 193 vi List of Tables 2.1 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.2 Grid setup for 2004 Sumatra Tsunami simulation . . . . . . . . . . . . . . 50 2.3 Fault parameters for 2004 Sumatra earthquake . . . . . . . . . . . . . . . 51 3.1 Wave and Current Conditions in Swan(1990)’s Experiments . . . . . . . . 85 3.2 Wave and current characteristics in Kemp & Simons(1982, 1983) . . . . . 91 3.3 Bichromatic waves and current characteristics in Dong et al.(2009) . . . . 97 4.1 Experiment conditions of pneumatic breakwater(Zhang et al.(2010)) . . . 125 4.2 Parameters for simplified density fields of bubbly area . . . . . . . . . . . 127 4.3 Experiment conditions of internal waves in two-layer(Umeyama(2002)) . . 134 vii List of Figures 1.1 A huge whirlpool generated by 2011 Tohoku Tsunami traps a boat in port of Oarai, Japan (Courtesy of Reuters/Kyodo) . . . . . . . . . . . . . . . . 4 2.1 Finite difference stencils for COMCOT . . . . . . . . . . . . . . . . . . . . 15 2.2 Finite volume stencil for Boussinesq Model . . . . . . . . . . . . . . . . . 21 2.3 Schematic drawing of coupled grid system(“Zone A” is specified for Fig- ure 2.4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 Data exchange schematic between COMCOT and Boussinesq grids (De- tailed view of “Zone A” in Figure 2.3) . . . . . . . . . . . . . . . . . . . . 26 2.5 Flowchart of coupled COMCOT-Boussinesq model calculations . . . . . . 27 2.6 Grid system for different grid sizes between COMCOT and Boussinesq Model (Lower-left corner section only) . . . . . . . . . . . . . . . . . . . . 31 2.7 Physical concept of Gaussian hump simulations; top: front view, bottom: plan view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.8 Definition sketch for L ∗ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.9 Temporal variation of free surface elevation at the center of the wave basin (r =1,ǫ=0.001,Cr =0.01 and e=10) . . . . . . . . . . . . . . . . . . . 38 2.10 Snapshots of water surface at 3 different times(r =1,ǫ=1,Cr =0.01 and e=10). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.11 Calculation of standard deviation when r = 2,ǫ = 0.001,Cr = 0.01 and e=10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.12 Error distribution according to stability index,γ, with fitted curve. . . . . 43 viii 2.13 Definitionsketchofalongwavepropagatingontoashallowshelf(LSW:Linear Shallow Water Equation, NLSW: Nonlinear Shallow Water Equation) . . 45 2.14 Timehistoriesofwatersurfaceelevationsat4differentlocations. Intheleft half of the figure are shown the coupled model results (solid line: Boussi- nesq model (Layer 4, dashed line: Measured data), while the right half of the figure shows the COMCOT only results at the same times (solid line: NonlinearShallowWaterEquationmodel(Layer3),dashedline: Measured data) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.15 Initial surface elevation of 2004 Sumatra tsunami . . . . . . . . . . . . . . 48 2.16 Surface elevation (m) in all layers at time=730min. The lower row shows theoutputfromLayer5,theBoussinesqmodel,ofthefreesurfaceelevation (left) and vorticity (1/s) (right). . . . . . . . . . . . . . . . . . . . . . . . 52 2.17 Vorticity (1/s) evolution inside Oman Salalah harbor at nine different times 53 2.18 Comparisonofvorticity(1/s)evolutionbyBoussinesq-coupledmodel(left) and COMCOT-only (right) . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.1 Definition sketch of long wave propagation over underlying currents of arbitrary profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.2 MeasuredshearcurrentsintheabsenceofwavesinSwan(1990);(a)CASE2F; (b) CASE2A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.3 Comparison of wave records under following, uniform current (Case 1F) betweenexperiment(o)andnumericalsolution(-): surfaceelevation(a)and oscillating velocities at z=-0.1m(b); at z=-0.2m(c); at z=-0.3m(d) . . . . 87 3.4 Comparison of wave records under opposing, uniform current (Case 1A) betweenexperiment(o)andnumericalsolution(-): surfaceelevation(a)and oscillating velocities at z=-0.15m(b); at z=-0.25m(c); at z=-0.35m(d) . . 88 3.5 Comparisonofwaverecordsunderopposing,linearshearcurrent(Case2A) betweenexperiment(o)andnumericalsolution(-): surfaceelevation(a)and oscillating velocities at z=-0.1m(b); at z=-0.2m(c); at z=-0.3m(d) . . . . 89 3.6 Comparisonofmaximumwavevelocityunderfollowing(a)andopposing(b) linear shear current; experiment(o), numerical solution(-) . . . . . . . . . 90 3.7 Mean-velocityprofilesofcombinedwavesandcurrents;(a)WCA1;(b)WCA3; (c)WCA4; (d)WCA5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 ix 3.8 Mean-velocityprofilesofcombinedwavesandcurrents;(a)WDR1;(b)WDR3; (c)WDR4; (d)WDR5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 3.9 Distribution of parameter b according to radiation stress(S xx scaled by current bed stress(τ Cb )); Fitted line is represented by solid and dashed line 96 3.10 Experimental setup of bichromatic waves under an ambient currents by Dong et al.(2009) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.11 Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbound- ing long waves;(without currents) . . . . . . . . . . . . . . . . . . . . . . . 98 3.12 Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbound- ing long waves;(following currents) . . . . . . . . . . . . . . . . . . . . . . 100 3.13 Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbound- ing long waves;(opposing currents) . . . . . . . . . . . . . . . . . . . . . . 101 3.14 Calculatedamplitudespectraofbichromaticwaveswithandwithoutcurrents102 4.1 Sketch of long wave propagation over two-layer fluids with horizontal den- sity variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.2 Experimental setup of pneumatic breakwater(Zhang et al.(2010)) . . . . . 125 4.3 Simplified density field of bubbly fluids; (a)Depth-varying ρ(z) (b)Depth- constant ρ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.4 Comparison of transmission coefficient(K t ) between measurement and cal- culation; T=1.55 sec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.5 Comparison of transmission coefficient(K t ) between measurement and cal- culation; T=1.29 sec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.6 Simplified two-layer system with interfacial wave 1 . . . . . . . . . . . . . 130 4.7 Measured(solid) and approximately calculated(dashed) density profile in Umeyama(2002) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 4.8 Measured(dashed), numerical(solid) and analytical(dot) surface elevation in CASE I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 1 Idea underlying the plot assumes internal-mode interfacial wave(η i ) generates internal-mode surface wave(η s ), For the definition of internal-mode, see P˘ ar˘ au and Dias(2001) x 4.9 Measured(dashed), numerical(solid) and analytical(dot) surface elevation in CASE II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 4.10 Second-order super-harmonic and sub-harmonic amplitude relative to the exact solution(a ± /a ± exact ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 xi Abstract In the coastal area, defined as the region between the shoreline and some offshore limit where the depth can no longer influence the waves, complex behavior of waves is an- ticipated due to various physical effects such as turbulence, wave-structure interaction, wave-current interaction, wave breaking and fluid-density variations. For modeling of nearshore hydrodynamics, many numerical models have been developed so far, but many of such effects are not yet considered appropriately. In this dissertation, depth-integrated numerical models used in long wave simulation are developed for better understanding of complicated hydrodynamics at the nearshore. First, a non-dispersive shallow water equation model and dispersive Boussinesq model are two-way coupled. The fundamental purpose of the coupling effort is to develop the capability to seamlessly model long wave evolution from deep to shallow water with fine scale resolution, without the loss of locally important physics. Second, a set of depth- integrated equations describing combined wave-current flows are derived mathematically and discretized numerically. To account for the effect of turbulent interaction between waves and underlying currents with arbitrary profile, new additional stresses are intro- duced, which represent radiation stress of waves over the ambient current field. Finally, a numerical model for gravity waves propagating over variable density fluids is developed xii by allowing horizontal and vertical variation of fluid density. Throughout the derivation, density change effects appear as correction terms while the internal wave effects on the free surface waves in a two-layer system are accounted for through direct inclusion of the internal wave velocity component. For each of the studied topics, numerical tests are performed to support accuracy and applicability. Consequently, we have developed a comprehensive tool for numerical simulation of complex nearshore hydrodynamics. xiii Chapter 1 Introduction 1.1 Motivation “We will never control the furious Earth, but through our scientific understanding of its nature, we may be able to prevent tragic and costly losses” 1 A series of recent gigantic disasters in the oceanic area, including the 2004 Indian tsunami and 2011 Tohoku tsunami remind us of the importance of understanding nature, especially ocean dynamics. Researchers have made great progress in developing tools for modeling and predicting the meteorological and oceanographic environments covering the immense ocean. Thinking of the lifespan of transoceanic waves, whether generated by wind stress or other source of energy in the deep ocean, travel over the ocean and reach at the shoreline after various types of physical transformations such as shoaling, 1 This quote is from Furious Earth: The Science and Nature of Earthquakes, Volcanoes, and Tsunamis. by Ellen J. Prager, McGraw-Hill, 1st ed., 1999 1 refracting, interacting with other flows or structures, and breaking, many of which occur at the nearshore region. Therefore, it is yet challenging to model a comprehensive evolution of waves from deep ocean to the shallow coastal area without loss of accuracy in physics. Three di- mensional Navier-Stokes equations will be the best choice, unless there is concern about computational resources or numerical issues. In reality, however, it seems too far from practical,sincewearenotcapableofdealingwiththecomplexityandmulti-scalebehavior without any assumptions aiming to simplify the problem. With the purpose of practi- cal application, many numerical tools have been developed for long wave modeling with differentlevelsofapproximation. EstablishedexamplesincludeCOMCOT(CornellMulti- grid Coupled Tsunami Model), MOST(Method of Splitting Tsunami), FUNWAVE(Fully Nonlinear Wave Model), COULWAVE(Cornell University Long and Intermediate Wave Modeling Package) and others. The majority of existing long wave models have their theoretical basis on either non- dispersive shallow water equations or Boussinesq-type equations. Shallow water equation models neglect the dispersive property of waves, permitting depth-uniform velocity pro- files, and a huge efficiency is gained through modeling with these 2HD(horizontal dimen- sion) equations. Yet, it is still argued that the application of a non-dispersive model in the nearshore area seems to be insufficient for accurate representation of certain physics, since nearshore hydrodynamics may involve many complex features. The dispersive property of intermediate waves need to be accounted for, as it can be important in the shallow regions in some circumstances. In the nearshore where the water depth is very shallow and amplitude and wavelength can become high and short, 2 nonlinear and bathymetric interactions occur across a wide range of frequencies. The in- teractionsamongwavesofvariousfrequenciescanlocallygeneratevariousshorter-crested, or dispersive waves components. Thus, an application of the dispersive wave model to the nearshore environment is considered useful and needed for practical purposes. Another complex situation occurs in coastal area when inhomogeneous types of fluid flows such as waves, currents and tides co-exist and interact with each other. Gravity waves arriving at coastal regions interact with the background current flow, which is usually driven by tides, thermocline, salinity variations and river mouth discharges. In consequence,thecombinedeffectsofwavesandcurrentsneedtobeconsideredforaccurate modeling of nearshore physics such as morphodynamic changes, mixing and transport of solutes. Freshwater and seawater are often found in estuaries in coastal regions. Thus the fluid densities are subject to change horizontally and vertically due to thermal and saline variability. This density variation will result in complex physical processes contrary to uniform density fluids. Therefore, density-variation is a primary concern in estuarine hydrodynamics. Wave-breaking, sedimentations, andwaveinteractionwithcoastalstruc- tures are worthy of mention as other interesting issues, as they often occur in nearshore area. The wave-related, physical problems mentioned above tend to make nearshore hydro- dynamics more complex than those in deep water. Sometimes, such dynamics result in unexplainable events, looking like supernatural phenomena and require dedicated consid- erations of the controlling physics. An excellent example can be found during the 2011 3 Figure 1.1: A huge whirlpool generated by 2011 Tohoku Tsunami traps a boat in port of Oarai, Japan (Courtesy of Reuters/Kyodo) Tohoku tsunami, which created a huge whirlpool when it reached a harbour near Oarai city, Japan(see Figure 1.1). This example describes how significantly the local physics can govern the nearshore hydrodynamics and how big their impacts are. Thus, accurate and computationally effi- cient modeling of nearshore hydrodynamics with consideration of physical complexities is extremelyimportant inthat imprecisepredictionwill cause tremendous lossofproperties and sometimes threaten our lives. Thinking of the rapidly-expanding populations toward shorelines, such efforts seem to be urgently in need. 4 Hydrodynamic complexities in the nearshore are recognized as important and so are popular topics in coastal and ocean engineering. However, comprehensive modeling of complexbehaviorofwavesisnotyetachievedandmostofcomplexitiesremainunrealized in present numerical models. As Boussinesq models include nonlinearity and frequency dispersionofwavesattheshallowregion,otherphysicalfactorsembeddedinnearshorehy- drodynamics can be elucidated physically and mathematically in the Boussinesq context. Thus it may be valuable to study various hydrodynamic complexities at the nearshore in this framework. 1.2 Objectives of Study The main objective of this research is to introduce a set of numerical models used in long wave simulation without the loss of locally important physics in the nearshore. To this end, we extended the conventional Boussinesq-type investigations in three ways. First, a non-dispersive shallow water model and dispersive Boussinesq model will be two-way coupled to develop a seamless model for long wave evolution from deep to shallow water with fine scale resolution, without the loss of locally important physics. In such an effort, it needs to be postulated that a more physically complete attempt at long wave mod- eling can be achieved through the integration of a shallow water equation model with a Boussinesq-type model. The former is computationally in charge of propagation of waves in the deep ocean, while the latter can be concentrated on a specific area of inter- est, typically nearshore where waves are prone to high nonlinearity and local frequency dispersion. 5 Second, a set of depth-integrated equations describing combined wave-current flows willbederivedmathematicallyanddiscretizednumerically. Inthenearshorearea,gravity wavesheavilyinteractnotonlywiththebottomgeometrybutalsowiththecurrentswhich are usually driven by tides, thermocline, salinity variations and river mouth discharges. To account for the effect of turbulent interaction between waves and underlying currents, additional stresses will need to be introduced to represent the intensity of turbulent interaction between waves and currents. Finally, a numerical model for gravity waves propagating over variable density fluids will be developed by allowing horizontal and vertical variation of fluid density. Through- out the derivation, density change effects as well as internal wave effects on the free surface waves will be accounted for and will appear as correction terms to the conven- tional Boussinesq equations of uniform density. For each of the studied topics, numerical tests will be performed to support its accuracy and applicability. 1.3 Organization In Chapter 2, non-dispersive shallow water and dispersive, Boussinesq-type numerical models used in long wave modeling, as well as an approach to two-way couple these models together, are introduced. The two model components are briefly introduced, and the physical mismatch between the two models is examined analytically. Then, a general benchmark test has been undertaken to provide a parameter range for expected accuracy andstabilityofthecoupledmodel. Finally,themodelisappliedtothe2004IndianOcean tsunami. 6 In Chapter 3, a set of depth-integrated equations describing combined wave-current flows is derived and validated. To account for the effect of turbulence induced by nonlin- ear interaction between waves and currents, additional stresses are introduced. Using a parameterb, an additional stress is defined to represent the intensity of turbulent inter- action between waves and currents. An appropriate estimation onb is provided through KempandSimons(1982,1983)’sexperiments. Accuracyofthemodelisexaminedthrough three experimental data sets, which resemble various types of hydraulic situations in the nearshore. In Chapter 4, by allowing horizontal and vertical changes of fluid-density, depth- integrated model equations for long surface waves over variable density fluid are derived mathematically and discretized numerically. Proposed model is applied to surface wave propagations over either horizontally or vertically varying density fluids for the verifica- tion. In Chapter 5, the conclusions of the dissertation are summarized. Also included are suggestions for the future works. 7 Chapter 2 Nested and Multi-Physics Modeling of Long Waves In this chapter, non-dispersive shallow water and dispersive, Boussinesq-type numerical modelsusedinlongwavemodeling,aswellasanapproachtotwo-waycouplethesemodels together are introduced. The fundamental purpose of the coupling effort is to develop thecapabilitytoseamlesslymodellongwave(e.g., tsunamis)evolutionfromgenerationto inundation with fine scale resolution, without the loss of locally important physics. The twomodelcomponentsarebrieflyintroduced,andthephysicalmismatchbetweenthetwo modelsisexaminedanalytically. Ascouplingofnumericallyandphysicallyheterogeneous models may result in undesirable errors, a general benchmark test has been undertaken toprovideaparameterrangeforexpectedaccuracyandstability. Longwavepropagation ontoashallowshelfissimulatedtovalidatethecoupledmodel,examiningtheimportance ofdispersiveandnonlineareffectsinthenearshorearea,aswellastheutilityofthecoupled modeling system. Finally, the model is applied to the 2004 Indian Ocean tsunami. In this test, the local dynamics experienced in the Port of Salalah in Oman, as documented by Okal et al. (2006), are recreated. 8 2.1 Introduction As a long gravity wave propagates over a non-uniform ocean bottom, shoaling and re- fraction can have a transforming effect. The response of waves to bathymetric changes results in deformation of the amplitude and wavelength, permitting waves to conserve mass and momentum. Many efforts have been made to construct a relationship between wave height and water depth, and using various levels of approximation (e.g. linear waves) it is possible. For the approximation of long and intermediate length (or shallow and intermediate depth) waves, two physical characteristics of waves, nonlinearity and frequency dispersion, are generally employed. Undera“true”longwave,frequencydispersionisnegligible. Thisassumptionyieldsa hydrostaticpressurefieldandahorizontalvelocitythatisuniformoverdepth. Atsunami is often considered a long wave. Frequency dispersion in a tsunami can be ignored when thetsunamiwavelength,typicallyontheorderof100kminthedeepocean,isconsiderably larger than water depth. Therefore, the usual approach to describe tsunami evolution is to take either the linear or nonlinear shallow-water models as the governing equations. A number of computational models based on this approximation exist, and some are introducedhere. MOST(MethodofSplittingTsunami)developedbyTitovandSynolakis (1998) is capable of predicting wave height or inundation using a technique where two- dimensionalequationsaresplitintoapairofone-dimensionalequations. Liuetal.(1998), on the other hand, presented COMCOT (Cornell Multi-grid Coupled Tsunami Model) adopting the staggered leap-frog integration with an upwind scheme for the nonlinear 9 convective terms. COMCOT can also model tsunami propagation and some nearshore- dynamics such as run-up. GTM (Global Tsunami Model) was designated for assessment of tsunami hazard, inundation, mapping and prediction of the tsunami arrival time by Kowalik et al.(2005). More recently, aided by adaptive finite volume methods for wave propagation, TsunamiClaw (Conservation Laws) has been created as a work of George andLeVeque(2006). Lastly,TsunAWIusesthefiniteelementmethodwiththeadvantage of flexibility in grid generation, and was found to be comparable to the multi-grid, nested modeling approach (Harig et al., 2007). Even though all of these models employ different numerical techniques, all solve the linear and/or nonlinear shallow water equations. Depending on the wavelength of the tsunami, however, frequency dispersion effects can be significant. Specifically, nei- ther hydro-static pressure nor depth-constant horizontal velocity can be presumed. For transoceanic propagation of a tsunami as well as landslide-generated tsunami, the dis- persive effects, estimated through the ratio of water depth to wavelength, should be included to yield more accurate results (e.g. Yoon(2002), Lynett et al.(2003), Grilli et al.(2007)). For this reason, some efforts to add the frequency dispersion effect into non- dispersive models through numerical truncation error have been made(e.g. Yoon(2002), Burwell(2007)). Despite such attempts to mimic physical dispersion, it is still an attrac- tive challenge to model tsunami with the Boussinesq or Navier-Stokes equations with the aim to, hopefully, obtain more realistic wave predictions. Corresponding examples include COULWAVE(Cornell University Long and Intermediate Wave Modeling Pack- age) by Lynett et al.(2003), GEOWAVE, which is equipped with the FUNWAVE(Fully 10 Nonlinear Wave Model) engine, by Grilli et al.(2007), and recent work by Saito and Furumura(2009). The above mentioned computational models aim, of course, for accurate prediction of nearshore physics such as wave shoaling, wave-diffraction and refraction, run-up and nonlinearinteractionsofwaves,eachwiththeirspecializedadvantages,e.g.,smallcompu- tational time of COMCOT or rigorous representation of physics in the Boussinesq model. In coastal regions, where the water depth is very shallow and amplitude and wavelength can become high and short, nonlinear and bathymetric interactions across a wide range of frequencies occur. These interactions can locally generate various shorter-crested, or dispersive waves components. A well known example is the transformation of a tsunami front into an undular bore. Thus, the nearshore is expected to be nonlinear and (pos- sibly) dispersive, and Boussinesq model is appropriate, as addressed in some literature (e.g. Lynett(2006)). Inrelatedefforts,Kimetal.(2009)havepresentedadepth-integrated model for weakly dispersive, turbulent and rotational fluid flows. This approach permits the explicit inclusion of viscous effects in shallow water, coupled with the nonlinear and weakly-dispersive physics of the Boussinesq model. With accuracy, this model can simu- latenonlinearandweaklydispersivenearshoredynamics,aswellaslargeeddiesgenerated by long waves and currents. Here it is postulated that a more physically complete attempt at tsunami model- ing can be achieved through the integration of a shallow water equation model with a Boussinesq model. COMCOT is computationally ”in charge” of generation and propaga- tion of tsunamis in the deep, open ocean, which in general will be the huge majority of a simulation domain. On the other hand, the Boussinesq effort can be concentrated on a 11 specific area of interest, typically nearshore where waves are prone to high nonlinearity, turbulence, and local frequency dispersion. Thischapteraimstointroduceasetofnumericalmodelsusedintsunamimodeling,as wellamethodtocouplethemtogether. Thefundamentalpurposeofthecouplingeffortis to develop the capability to seamlessly model tsunami evolution from generation to inun- dation with fine scale resolution, without the loss of locally relevant physics. In addition to this, local turbulent structures, such as eddies and gyres, generated by tsunamis in the nearshore area or around coastal structures can be studied with the coupled model. The outline of this chapter is as follows. In the first two sections, the model com- ponents, the shallow water wave equation model (COMCOT) and the Boussinesq-type model, will be briefly introduced with their numerical scheme. Physical and numerical “mismatches” between the two models will be discussed analytically, which is followed by the coupling method given in detail. The next section is devoted to the validation of the coupled model through a benchmark test with wide-varying conditions and re- sulting guidance for general use. A typical problem of long wave propagation into the coast is given in following section. Finally, the presented model is applied to the 2004 Sumatratsunamitoinvestigatenearshoredynamics,withaparticularfocusontheeddies generated inside a harbor basin. 2.2 Shallow Water Equation Model As introduced above, Liu et al.(1998) presented a nested multi-grid model which has the option of using either the linear or the nonlinear shallow water equations (NLSW) with 12 two different types of coordinate systems, namely Cartesian or spherical. This general framework includes the effects of bottom friction as well as a special treatment for the moving shoreline. The model named COMCOT v.1.6 has been adapted here to simulate tsunami propagation across oceanic basins. 2.2.1 Governing Physics The nonlinear shallow equations including bottom frictional effects in conservative form are: ∂ζ ∂t + ∂M ∂x + ∂N ∂y =0 (2.1) ∂M ∂t + ∂ ∂x M 2 H + ∂ ∂y MN H +gH ∂ζ ∂x +τ x =0 (2.2) ∂N ∂t + ∂ ∂x MN H + ∂ ∂y N 2 H +gH ∂ζ ∂y +τ y =0 (2.3) in which ζ is the surface elevation, h is the still water depth, H = h +ζ is the total water depth, and M, N are the volume fluxes in the x and y directions, defined respec- tively as Hu and Hv. The bottom friction terms τ x ,τ y in the momentum equations are approximated in COMCOT via Manning’s formulation τ x = gm 2 H 7/3 M M 2 +N 2 1/2 (2.4) τ y = gm 2 H 7/3 N M 2 +N 2 1/2 (2.5) 13 wheremistheManning’srelativeroughnesscoefficient. Notethattheaboveequationset is the nonlinear solver in COMCOT; the linear solver of course neglects the convection terms in the momentum equations, and does include the Coriolis force when solving in spherical coordinates. 2.2.2 Numerical Scheme ThenumericalschemeemployedbyCOMCOTistheexplicitleap-frogdifferencemethod. Nonlinear terms in the model are approximated with upwind finite differences and linear termsbytwo-pointcenteredfinitedifferences. Thisnumericalschemeisstableandrobust butisalow-orderaccuratemethod,meaningthatitissusceptibletonumericaldispersion and dissipation errors. The finite difference forms for the continuity and momentum equations are described in Appendix A. The finite difference stencil of this scheme is depicted in Figure 2.1, suggesting two neighboring points on each side of a calculation point are necessary for each location calculation of derivatives. For the present study, COMCOT has been parallelized for use on shared-memory computers, such as multi-processors and/or multi-core computers. OpenMP was used for the parallelization, which is the standard method for shared-memory parallelization. The parallel model has been tested up to 8 processors, and shows a near linear speed-up (using 8 processors reduces CPU time by a factor of 1/8). Togeneratethetsunamifromanunderseaearthquake,COMCOTusesthefaultmodel ofOkada(1985). Themainassumptionsofthismodelarearectangularfaultplanewithin an elastic deformation. The fault model predicts the deformation of the seafloor, which corresponds directly to the initial deformation of the ocean water free surface. There are 14 Conservation of Mass 1/2 , n i j ! 1/2, n i j M 1/2, n i j M , 1/2 n i j N , 1/2 n i j N i i+1 i+1/2 i-1 i-1/2 i-3/2 i+3/2 j j-1 j+1 j-3/2 j-1/2 j+1/2 j+3/2 i i+1 i+1/2 i-1 i-1/2 i+3/2 j j-1 j+1 j-3/2 j-1/2 j+1/2 j+3/2 Conservation of Momentum(x-direction) 1 1/2, n i j M 1/2 , n i j ! 1/2 , 1 n i j ! " 1/2 1, 1 n i j ! ! " 1/2 1, n i j ! ! 1/2 1, 1 n i j ! ! ! 1/2 , 1 n i j ! ! 1/2 1, n i j ! " i+2 1/2 2, n i j ! ! 3/2, n i j M 1/2, n i j M , 1/2 n i j N , 1/2 n i j N 1, 1/2 n i j N ! 1, 1/2 n i j N 1 , 1/2 n i j N 1/2 , n i j ! 1/2 1, n i j ! " 1/2 1, 1 n i j ! " ! 1/2 1, n i j ! ! 1/2 1, 1 n i j ! ! ! 1/2 , 1 n i j ! ! i i+1 i+1/2 i-1 i-1/2 i-3/2 i+3/2 j j-1 j+1 j+2 j-1/2 j+1/2 j+3/2 Conservation of Momentum(y-direction) , 1/2 n i j N , 3/2 n i j N 1/2 , 1 n i j ! " 1/2 , 2 n i j ! ! 1/2, n i j M 1/2, n i j M 1/2, 1 n i j M ! 1/2, 1 n i j M Figure 2.1: Finite difference stencils for COMCOT 15 a number of parameters which govern the fault model: • Latitude and longitude of epicenter • Focal depth • Length and width of fault plane • Dislocation • Strike angle • Slip angle • Dip angle Once the earthquake has been described with the above parameter set, COMCOT is able to propagate the initial disturbance across oceans. For propagation across deep ocean waters, COMCOT gives the option of using the linear version of the shallow water equations. This version is solved considerably faster, in the computational sense, than the nonlinear version, and can be used with confidence as long as the tsunami wave height is a very small fraction of the depth, practically less than 1/25∼1/50 of the local depth. Whenthisthresholdisexceeded,thenonlinearversionofCOMCOTisrequiredfor accurate results. Generally, if runup or nearshore wave heights are needed, the nonlinear version of the model should be used. For runup calculations, COMCOT utilizes a simple but accurate moving shoreline algorithm. The continuously sloping beach profile is approximated as a stair-stepped profile. When the water level exceeds the elevation of the “stair” above, the water floods 16 that “step” and the shoreline moves landward (inundation). This approach has shown to re-produce analytical solutions reasonably and field data as well as any other published, shallow-water equation model. 2.3 Boussinesq Equation Model Recently,Kimetal.(2009)havepresentedadepth-integratedmodelforweaklydispersive, turbulentandrotationalfluidflows. Itisderivedfromthespatially-filteredNavier-Stokes equationsinordertoconsiderviscouseffectsofaturbulent fluid. Accordingly, thismodel includes approximated bottom-induced turbulence and thereby the associated vertical and horizontal rotational effects can be captured. In the present study, we have adopted the Boussinesq model of Kim et al.(2009) to simulate the nearshore hydrodynamics and turbulence effects such as large eddies and wakes generated in the nearshore. 2.3.1 Governing Physics TheBoussinesq-typeequationsincludingturbulentviscosityandtheassociatedhorizontal and vertical vorticity terms are given in conservative form below: ∂H ∂t + ∂HU α ∂x + ∂HV α ∂y +D c =0 (2.6) ∂HU α ∂t + ∂HU 2 α ∂x + ∂HU α V α ∂y +gH ∂ζ ∂x +HD x m +U α D c =0 (2.7) 17 ∂HV α ∂t + ∂HU α V α ∂x + ∂HV 2 α ∂y +gH ∂ζ ∂y +HD y m +V α D c =0 (2.8) where U α and V α are the x and y component velocities at z α = −0.531h and D x m ,D y m are 2nd order correction terms of the depth-integrated momentum equations as defined in Kim et al.(2009). Likewise D c includes 2nd order correction terms in the continuity equation. It is noted that the dispersive, viscous, and vorticity corrections are included as these 2nd order terms. All 2nd order terms can be found in Appendix B. 2.3.2 Numerical Scheme To numerically solve the governing equations in conservative form, a highly accurate and stable model is developed. The numerical method uses a fourth-order MUSCL-TVD (Monotone Upstream-centered Schemes for Conservation Laws- Total Variation Dimin- ishing) scheme to solve the leading order (shallow water) terms, while for the dispersive terms, a cell averaged finite volume method is implemented. For the time integration, a third order Adams-Bashforth predictor and the fourth-order Adams-Moulton correc- tor scheme has been used to keep numerical truncation errors small. It is noted that Boussinesq-type models such as the one solved here, which include up to third-order spa- tial derivatives, require a high-order solution scheme to keep the derivatives associated with the numerical truncation error at least an order below those contained in the model equations. 18 The explicit predictor step is ζ n+1 =ζ n + Δt 12 23E n −16E n−1 +5E n−2 (2.9) P n+1 = P n + Δt 12 23F n −16F n−1 +5F n−2 +2F n 1 −3F n−1 1 +F n−2 1 +F p v (2.10) Q n+1 = Q n + Δt 12 23G n −16G n−1 +5G n−2 +2G n 1 −3G n−1 1 +G n−2 1 +G p v (2.11) The implicit corrector step is ζ n+1 =ζ n + Δt 24 9E n+1 +19E n −5E n−1 +E n−2 (2.12) P n+1 = P n + Δt 24 9F n+1 +19F n −5F n−1 +F n−2 +F n+1 1 −F n 1 +F c v (2.13) 19 Q n+1 = Q n + Δt 24 9G n+1 +19G n −5G n−1 +G n−2 +G n+1 1 −G n 1 +G c v (2.14) where P, Q are defined as P =HU α + H 2 z 2 α −ζ 2 U αxx + H(z α −ζ)(hU α ) xx − Hζ x {ζU αx +(hU α ) x } (2.15) Q=HV α + H 2 z 2 α −ζ 2 V αyy + H(z α −ζ)(hV α ) yy − Hζ y n ζV αy +(hV α ) y o (2.16) where the superscript n denotes time level and the subscripts x and y imply derivatives in the x and y direction, respectively. E, F, G, F 1 , G 1 , F p v , G p v , F c v , G c v in the above equations include a number of spatially discretized terms; all can be found in Appendix C. The finite volume stencil for this scheme is displayed in Figure 2.2, which shows that 4 neighboring points are required for each local calculation. 2.4 PreliminaryDiscussionforMismatchesbetweenModels It is necessary to compare models in terms of both physical limitations and numerical properties, as this will provide the basic guidance for coupling. There exist two major 20 i j Figure 2.2: Finite volume stencil for Boussinesq Model groups of errors, one from approximated equations and the other from the numerical scheme, the so-called truncation error. Prior to our coupling work, we shall thus consider the physical and numerical differences between COMCOT and the Boussinesq model, since the models, of course, have different governing equations as well as very different numericalsolutionschemes,whichwillcauseaphysicalandnumericalaccuracymismatch along the coupling interface. 2.4.1 Mismatch in Physics COMCOT, based on the shallow water equations, approximates the horizontal velocities and pressure gradient to be constant with depth, so one can ignore vertical variation of physics. On the other hand, Boussinesq-type equations allow (weak) vertical change of horizontal velocity, expressed as quadratic function of z. This allowance of vertical variation in the flow permits the model to include the effects of frequency dispersion. In 21 addition to this, Kim et al.(2009) included viscous and rotational effects, which originate from bottom-induced stress. Both frequency dispersion and viscous effects exist as cor- rections to the leading order terms - the inviscid nonlinear shallow water equations. This can be seen straightforwardly by eliminating 2nd order terms in Equations (2.6), (2.7) and (2.8) and then comparing with Equations (2.1), (2.2) and (2.3). Also note that the bottom friction terms in both models are included, but in a different manner; COMCOT has ad-hoc added friction terms in the momentum equations, while the Boussinesq model has both the bottom stress and a number of additional terms resulting from an explicit inclusion of bottom stress in the derivation. Consequently, in order to avoid errors from these physical differences, dispersive and viscous effects should be sufficiently small in order for governing physics to be continuous across the model interface. In other words, physics-driven model errors can be mitigated when the local relative depth(h/L) and bottom friction along interface are small. Of course in a general, nonlinear simulation, one can not for example guarantee the minimum simulated value of L a priori, but such relations can be used as guidance for constructing a simulation. 2.4.2 Mismatch in Numerics Different numerical solution schemes will produce different output even if solving the same algebraic equations. Discretizing equations using any sort of numerical method includes an error from truncation, and such errors of course depend on the numerical scheme itself. As explained in the previous section, COMCOT and the Boussinesq model employ different types of schemes; while both are fixed grid solvers, the spatial stencils of the two schemes are very different, as are the time integration methods. 22 Again, the Boussinesq solver is a fourth-order accurate scheme and COMCOT a second-order scheme. Since all numerical truncation errors from the Boussinesq solver are at least two orders of differentiation higher than those from the COMCOT solver, the primary numerical error will originate from the nonlinear COMCOT model. Specifically, the leading order truncation error arising from the upwind differencing can be given as: Er =0.5(1−Cr)uΔx ∂ 2 u ∂x 2 (2.17) where Er is the numerical truncation error of the upwind difference and Cr is the local courant number, given as (Δt √ gh)/Δx. In general, it is not possible to ensure Cr ≈ 1 in COMCOT for an arbitrary grid with variable bathymetry, and furthermore the maximum allowable Cr in the COMCOT scheme according to stability analysis is 0.7 (in the Boussinesq it is 0.5 according to Kim et al. (2009) ). It is reasonable to assume that Er ≈ O(uΔx ∂ 2 u ∂x 2 ) for a generic geophysical simulation. Thus, the only solution to ensureaprecisenumericalmatchacrossthecouplinginterfaceistohaveatruelongwave, with negligible depth-averaged velocity curvature in the horizontal plane, at the interface location. This conclusion is not surprising, as the shallow water based COMCOT model has this as a general requirement for accuracy in any and all applications. By requiring that, at the interface, the physics represented by COMCOT are valid for that location, the mismatch in numerics essentially vanishes. If one was able to isolate the dispersive and viscous effects in the Boussinesq domain, and such effects were small at the coupling interface, the model matching would be best. 23 From above discussion, it can be reasonably expected that the primary errors ex- pectedatthecouplinginterfacearedrivenbyphysicaldifferencesinthemodelequations, and specifically use of the COMCOT model for a hydrodynamic situation where it is, strictly speaking, beyond its physical validity limits. Note that possibly large numerical differences might also exist; however as discussed above these arise when the wave is not practically long, and the velocity curvature is not negligible. This expectation implies that some sort of special treatment to deal with the physics mismatch will be required at the interface. This will be described in more detail in the next session. The general approach will be to turn “off” of the high order terms in the Boussinesq model at the interface and slowly turn them back “on” as one moves inside the Boussinesq domain. The remaining model mismatch errors, controlled by slightly different viscous closures and very different numerical schemes, will be mitigated by spatial filtering. Stability and accuracy of this special treatment will be discussed as well. 2.5 Coupling Approach The coupling method in which the shallow water model and Boussinesq model are inte- grated is presented here. The constituents are two-way coupled. Boundary conditions on the interfacing side of each model are provided by its counterpart model through data exchange and overlapping grid points. 2.5.1 Coupling Method Toaccommodatedataexchangebetweenthetwomodels,thecomputationalgridsofboth models should be overlapping, as shown in Figure 2.3. Since each model has derivatives 24 Domain calculated by COMCOT Boundary condition from Boussinsq model Boundary condition from COMCOT Domain calculated by Boussinesq COMCOT domain Boussinesq domain Imaginary domain for coupling Zone A Figure2.3: Schematicdrawingofcoupledgridsystem(“ZoneA”isspecifiedforFigure2.4) of different order in the corresponding governing equations, they each need a different number of overlapping points. These overlapping points act as exterior boundary con- ditions on the spatial edges of the computational domain; they are effectively imaginary grid points with data values taken from the neighboring model. As seen in Figure 2.4, COMCOT needs two points as a boundary condition whereas the Boussinesq model re- quires four neighboring points due to the 4th order MUSCL-TVD scheme and 3rd order spatial derivatives. Also, special attention must be paid to the calculation of velocity as each model defines velocity (or flux) at a different location relative to a cell (grid points are defined at the center of a cell). In the Boussinesq model, surface elevation as well as 25 1 2 3 4 5 6 N-4 N-3 N-2 N-1 N N+1 1. Calculation of COMCOT domain 2. Pass information from COMCOT to Boussinesq 3. Calculation of Boussinesq domain 4. Get back information from Boussinesq Boussinesq grid COMCOT grid Imaginary grid Surface elevation and velocity at Boussinesq grid Surface elevation at COMCOT grid “Time-averaged” flux at COMCOT grid Transfer of information Figure 2.4: Data exchange schematic between COMCOT and Boussinesq grids (Detailed view of “Zone A” in Figure 2.3) velocity components are defined at center of the cell based. On the other hand, flux has been placed at the interface of each cell in COMCOT model due to the staggered grid. With this concept of the interface treatment, we propose the calculation algorithm as shown in Figure 2.5. The algorithm consists of two main parts. The first part is the COMCOTmodelcalculationontheleftsideofFigure2.5andtheotheristheBoussinesq modelcalculationontherightside. Theyexchangedataeverytimestepthroughtwo-way coupling, as indicated by the boxes in the middle of Figure 2.5. The coupling algorithm can be explained in more detail by proceeding step by step. Let it be assumed that information up to time level t = t n is known. Note that the numerical scheme of COMCOT is staggered not only in space but also in time. Therefore 26 Read data Solve continuity equation Boundary condition (Using information at the imaginary cells) Solve conservation of momentum equation Determine initial conditions (ζ, U α ,V α ) 1/2,3/2,5/2 for Boussinesq model Update variables in COMCOT ζ n-1/2 = ζ n+1/2 (M, N) n =(M, N) n+1 Pass (ζ, U α ,V α ) n+1/2 at imaginary cells for Boussinesq Predictor for (ζ, U α ,V α ) n+1/2 Corrector for (ζ, U α ,V α ) n+1/2 Update variables in Boussinesq (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-7/2 = (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-5/2 (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-5/2 = (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-3/2 (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-3/2 = (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-1/2 (ζ, U α ,V α ,E,E1,F,F1,G,G1) n-1/2 = (ζ, U α ,V α ,E,E1,F,F1,G,G1) n+1/2 Start End Repeat until error<10 -4 N=3,….., N 2-way coupling COMCOT Boussinesq model Pass (ζ n+1/2 ,U α n ,V α n ) at imaginary cells for COMCOT Figure 2.5: Flowchart of coupled COMCOT-Boussinesq model calculations 27 in COMCOT, flux is known at time level t = t n and the surface elevation at time level t = t n−1/2 . For the Boussinesq model, the surface elevation and velocities at time level t=t n−1/2 are known. All dependent variables in the Boussinesq model are calculated at the same time level, t = t n+1/2 , whereas COMCOT simulates the surface elevation and flux terms at the time level t =t n+1/2 and t =t n+1 , respectively, following the leap frog scheme. The coupled solution scheme is outlined below. 1. Calculate COMCOT free surface at time level t = t n+1/2 by solving the NLSW continuity equation with flux information at time levelt=t n and surface elevation at time level t=t n−1/2 . 2. Calculate all Predictor step values in the Boussinesq model, yielding initial predic- tions at time level t=t n+1/2 for surface elevation as well as velocity. 3. Transfer Predictor Boussinesq surface elevation values and fluxes along the inter- face at time level t =t n+1/2 and t =t n , respectively, into COMCOT as boundary conditions. Note that Boussinesq flux term should be interpolated not only spa- tially but also in time. Linear one-dimensional interpolation is used, numerically consistent with the upwind differencing in COMCOT. If the grid size in each model is not constant, utilize a two-dimensional (bi-linear) interpolation technique to give appropriate boundary condition at the interface. 4. Calculate the COMCOT flux at time level t = t n+1 by solving the shallow water momentum equation with surface elevation information at time level t = t n+1/2 along with flux at time level t =t n . In this step, the information transferred from Boussinesq model has been used. 28 5. Extract the boundary condition at time level t = t n+1/2 to be transferred from COMCOT into the Boussinesq model. Again, note that the Boussinesq model does not employ a staggered-grid in time. Having already calculated the COMCOT surface elevation at time level t = t n+1/2 , our interest will focus on getting the COMCOT flux at time levelt=t n+1/2 . This can be done simply by taking average two values both at time level t=t n and t=t n+1 . 6. Transfer the COMCOT surface elevation and flux along the interface at time level t = t n+1/2 into the Boussinesq as boundary conditions required for the implicit corrector step. 7. Calculate the surface elevation and velocities at t = t n+1/2 from the Boussinesq Corrector using COMCOT boundary conditions at time level t=t n+1/2 . 8. Optionally, a filtering technique can be applied in order to remove spurious two- grid wave components with high frequency. See the “Special Numerical Interface Treatment” section below for details. 9. Return to step 1 for the next time step. 2.5.2 Coupling with Different Grid Size Coupling with different grid sizes can also be accommodated. For the estimation of information at an interface, an interpolation technique has been used. If the relative 29 position of a desired point in a grid is known,f a at the desired point(x a ,y a ) is calculated by using bilinear interpolation as follows. f a =(1−t)(1−u)f BL +t(1−u)f BR +tuf TR +(1−t)uf TL (2.18) in which t≡ x a −x L x R −x L , u≡ y a −y B y T −y B (2.19) where the subscript L,R,B and T in the above definition respectively means left, right, bottom and top. Figure 2.6 shows the grid system for coupling with different grid sizes. 2.5.3 Special Numerical Interface Treatment Each of the two numerical models coupled here has its own governing equations (ap- proximated physics) and numerical scheme, and this can result in both a physical and numericalmismatchalongtheinterface,asdiscussedabove. Fromaqualitativeanalysisof theequationsandnumericaldifferences,itisextrapolatedthatlargenumericaldifferences arise only when significant equation (physics) differences exist. This expectation was ob- served during early testing of the interface. To reduce this primary error dependency on the mismatch of model physics, the high-order dispersive terms in the Boussinesq model are neglected at the interface. These terms are linearly ramped back into the equations over a length of 20 grid points moving away from the COMCOT interface. While forcing the two numerical schemes to solve similar governing equations at the interfaceeliminatesalargefractionofthe interfaceerror, thedifferent viscoustreatments and numerical solvers will still incur some error in the simulation. This error commonly 30 1 2 3 4 5 6 7 8 Surface elevation of COMCOT (at the centre of grid) Flux of COMCOT (at the boundary of grid) Surface elevation and flux of Boussinesq Model (at the centre of grid) 8 7 6 5 4 3 2 1 Imaginary Grid for Coupling Imaginary Grid for Coupling Figure2.6: GridsystemfordifferentgridsizesbetweenCOMCOTandBoussinesqModel (Lower-left corner section only) takestheformofspurioushigh-frequency,two-gridwavecomponents;inessencenumerical noise reflected off “improper” boundary conditions. To remove this spurious two-grid wave component, the nine-point spatial filter is employed as suggested in Shapiro (1970). ζ i = 1 256 {186ζ i +56(ζ i−1 +ζ i+1 )−28(ζ i−2 +ζ i+2 ) +8(ζ i−3 +ζ i+3 )−(ζ i−4 +ζ i+4 )} (2.20) 31 Note that application of this filter is common in high-order Boussinesq applications, which tend to be prone to high frequency instability (e.g. Gobbi and Kirby, 1999), and is not used in the COMCOT domain. The filter is typically applied throughout the entire Boussinesq domain once every 100 time steps. 2.6 Validation As discussed above, coupling of two heterogeneous models is subject to the generation of undesirable errors. These errors are a function of wave nonlinearity and dispersion, and are difficult to quantify directly for model operation. To define validity for practical application, a general benchmark test is proposed with various initial, geometric, and nu- merical conditions. For this simulation experiment, COMCOT’s nonlinear shallow water equation model is coupled with the fully nonlinear version of the Boussinesq model. This approach is taken in order to examine the typical applicability space, since this combi- nation can be regarded as both the most physically well-matched coupling (compared to using the linear COMCOT), and likely also the most common matching setup. The ad-hoc modifications presented in the “Special Numerical Interface Treatment” are not used here, to provide a more conservative result. Finally, output from these “validation” simulations are evaluated, with a strong focus on stability and accuracy. 2.6.1 Gaussian Hump Simulation The Gaussian hump initial condition is very useful for this test as the resulting water surface disturbance radiates in all directions, forcing cross-derivatives in model to be 32 non-zero; there is no dominant propagation direction. The Gaussian hump used here has the initial free surface condition defined as, ζ(x,y,t)=H 0 exp − 1 α 2 (x−x 0 ) 2 +(y−y 0 ) 2 (2.21) whereH 0 is the initial height of the hump at its center (x 0 ,y 0 ) andα is the characteristic horizontal lengthscale of the Gaussian hump. All initial velocities are set to zero. The analytical solution for this case can be derived with the assumption of a small amplitude wave, or equivalently a linearization of the governing equations. A solution using Fourier decomposition can be found in Wei et al.(1995). 2.6.2 Physical and Numerical Setup The test cases include various physical conditions with different initial H 0 /h and α/h, in order to consider the effect of nonlinearity and frequency dispersion on the results. Figure 2.7 shows the physical layout of the basin and computational grids. The length of the basin has been fixed at 100m over all the simulations while the various other parameters are changed. Nonlinearity is expressed in the typical format asǫ=H 0 /h and dispersion as µ =h/L ∗ . Here L ∗ , is a characteristic length scale of the initial condition, given as the length between wave points at 5%H 0 , and graphically defined in Figure 2.8. The tested range of this parameter is 0.0002 to 0.0193, a range common for tsunamis where significant dispersion is not expected to be important. µ is not studied in this parametric exercise; these results are only valid for incident long waves where COMCOT 33 L y =100m L x =100m Boussinesq domain COMCOT domain 50m 50m h H wall wall Figure 2.7: Physical concept of Gaussian hump simulations; top: front view, bottom: plan view 34 H 0.95H L* Figure 2.8: Definition sketch for L ∗ can be expected to yield accurate predictions. Nonlinearity ranges from ǫ = 0.001 to ǫ=1.0. Along with different combinations of physical conditions, numerical parameters are varied as well. These numerical parameters, represented by the grid size ratio be- tween COMCOT and Boussinesq model(r = Δx C /Δx B ) and the CFL condition(Cr = √ ghΔt/Δx) are set to cover a wide range of possible configurations. Additionally, the ratio of the Boussinesq grid size to the water depth(e = Δx B /h) is tracked to examine some observed instability possibly due to a relatively small Boussinesq grid size. Consequently, throughout the simulations, four dimensionless parameters (r,ǫ,Cr, and e) are controlled so as to characterize factors affecting the numerical results. Those parameters and their ranges are listed in Table 2.1, producing 320 unique parameter combinations. For any four parameter set, a simulation is completely described in terms 35 Table 2.1: Simulation Setup r ǫ Cr e 1 0.001 0.01 1 2 0.01 0.05 2 4 0.1 0.1 5 6 1 0.5 10 10 - - - of its physical and numerical configuration. The order by which all the parameters can be calculated is as follows: 1. Calculate Δx B using r (Δx C =1.0m, nx C =100 throughout all the simulations) 2. Next, the number of grid points in the Boussinesq domain can be determined as nx B =r×(50−1)+1+4+4. Note that for imaginary, overlapping grids, 4 grids are added on each edge 3. Calculate h using e 4. H 0 can be obtained by ǫ(=H 0 /h) 5. The initial Gaussian surface is generated using α=15 m 6. Finally, Δt is determined using Cr(= √ ghΔt/Δx B ) whereΔx,nxdenotegridsizeandgridnumber,respectively,Δtthetimeincrement,hthe water depth, and H 0 the height of hump, respectively. The non-dimensional simulation time, t ′ (=t √ gh/L ∗ ) = 30 has been used for all cases; all individual cases have approxi- matelythesamenumberofwatersurfacefluctuations(characteristicperiods)duringeach simulation. 36 Among the 320 runs, the most computationally expensive case will be the one with r = 10,Cr = 0.01, and ǫ = 1.0 which corresponds to Δx B = 0.1m, Δt = 0.01sec and nx B = 499. On the other hand, the set having r = 1,Cr = 0.5 and e = 10 which gives Δx B = 1.0m, Δt = 0.505sec and nx B = 58 will be the most rapid. Runs required a few minutes to a few hours to complete on a desktop computer; the longest runs where those which developed numerical instabilities, causing excessive iteration of the Boussinesq corrector step. 2.6.3 Simulation Results Comparison with an analytical solution is, of course, the ideal method to evaluate the accuracy of a numerical model. However, the aforementioned analytical solution is be obtained using linearization; most of the (ǫ) values tested would violate the required assumption for linearization. Hence, the analytical solution will only be used to show that the coupled numerical model is producing accurate results for the smallest (ǫ) cases. Figure 2.9 shows the time series of water surface elevation at the center of the Boussiensq domainwhenr =1,ǫ=0.001,Cr =0.01ande=10. Thecomparisonbetweennumerical and analytical data is excellent for this small amplitude case and the coupled model works quite well. Snapshots of water surface elevation when r = 1,ǫ = 1,Cr = 0.01 and e = 10 are also given in Figure 2.10. These snapshots are showing both the Boussinesq and COMCOT surfaces, and there are no evident numerical errors along the interface, or anywhere else in the domain. Now, interest turns to evaluating output from all 320 simulations in some characteris- tic way. Here, the time series of surface elevation at the center of the domain will be used 37 0 200 400 600 800 1000 1200 1400 1600 −5 0 5 10 x 10 −5 time(sec) ζ(m) Numerical Analytical Figure 2.9: Temporal variation of free surface elevation at the center of the wave basin (r =1,ǫ=0.001,Cr =0.01 and e=10) 38 20 40 60 80 100 20 40 60 80 100 −0.1 −0.05 0 0.05 0.1 t=0s x(m) y(m) z(m) 20 40 60 80 100 20 40 60 80 100 −0.1 −0.05 0 0.05 0.1 t=70s x(m) y(m) z(m) 20 40 60 80 100 20 40 60 80 100 −0.1 −0.05 0 0.05 0.1 t=300s x(m) y(m) z(m) Figure 2.10: Snapshots of water surface at 3 different times(r = 1,ǫ = 1,Cr = 0.01 and e=10) 39 as the comparison basis. Considering that use of the linearized analytical solution is not helpful for assessment of highly nonlinear cases, some other available data for accuracy evaluation must be found. Here, a “base” case for eachǫ andµ combination is proposed. Each base case uses r = 1 and Cr = 0.01 , which is expected to produce most accurate and stable result for any physical parameter combination. Therefore, each simulation is able to be evaluated relative to its base case, assumed to be the “correct” solution. Note, however, that with this approach, it is not possible to state with absolute confidence that the “base” case represents an accurate solution, only a stable, converged one. It is likely thatstabilityandaccuracyoccurintandem, butthisisnotguaranteed. Simulationsthat closely resemble the “base” case can be considered stable simulations, while simulations that do not resemble the ”base” case are most likely unstable and inaccurate. Each nu- merical result is rated with respect to the root mean square (RMS) difference from its base case averaged through the last two dimensionless time units, from t ′ = 28 to 30 . An example depicting this procedure is shown in Figure 2.11. Note that if a particular simulation crashes due to instability, this RMS difference is set to 1.0. Despite the straightforward assessment of the simulation data, it is not a simple matter to demonstrate the relation between the error and the 4 parameters. As indicated inTable2.1andmentionedearlier,the4parametersareinterrelated,yetaffecttheresults in different ways. For an explicit measure of accuracy and stability using the 4 correlated parameters, adimensionless parameter, the stability indexγ, isintroduced. The stability index is a product of r,ǫ,Cr, and e and is expressed as γ =r a ǫ b Cr c e d (2.22) 40 0 500 1000 1500 2000 −0.5 0 0.5 1 time(sec) ζ / H (−) Normalized surface elevation raw base 0 500 1000 1500 2000 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 (RMS Error=0.0092044 Correlation=0.99979) time(sec) Difference Figure 2.11: Calculation of standard deviation when r = 2,ǫ = 0.001,Cr = 0.01 and e=10 41 where the exponents a,b,c and d are to-be-determined by nonlinear regression analysis on the processed data set. Through the stability index γ, the accuracy and stability properties of a coupled model simulation can be characterized, approximately, before said simulation is run. The higher the γ value, the more poorly behaved the simulation should be. The exponentsa,b,c andd are found through a best fit with the RMS error of each simulation; the calibrated stability index should yield the RMS error expected late in a simulation. Through the nonlinear regression analysis, a,b,c and d are found as 0.0001, 0.0328, 0.0621 and -0.0040, respectively. This implies that the Courant number (Cr) is the most dominant factor in a simulation, whereas the stability of the coupled model is insensitive toboththegridratiofactor(r)andthedepth-scaledgridlengthintheBoussinesqmodel (e). Figure 2.12 shows the error distribution according to γ. It is not surprising that the distribution of RMS error in Figure 2.12 resembles a typical cumulative probability curve. There is a transition between stable simulations (RMS error 0) and unstable, inaccurate simulations (RMS error 1), and this error “accumulates” with increasing γ. From this result, it can roughly but conservatively be stated that the stability indexγ is recommended to be less than 0.9 to provide a high likelyhood of a stable simulation. 2.7 Tsunami Wave Fission Simulation Asmentioned inthe Introduction, the Boussinesq model has the ability toyieldareason- ably complete representation of coastal hydrodynamics. In this section, which is focused on the demonstration of these properties through an efficient use of the coupled model, 42 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 γ =r 0.0001 ǫ 0.0328 Cr 0.0621 e −0.0040 RMS error data fitted sigmoidal curve Figure 2.12: Error distribution according to stability index,γ, with fitted curve 43 model verification will be performed using laboratory data. Matsuyama et al.(2007) conducted an experiment in a 205m long channel using a large and undistorted scale to investigate tsunami shoaling on the continental shelf. The incident waveform uses a sinusoidal-shape wave with single cycle defined as ζ = Asin 2π T t , 0≤t≤T = 0, t>T (2.23) where A is wave amplitude, T is period, and t is time. A single experimental case with A=0.03m and T=20sec, which exhibits significant tsunami shoaling on the continental shelf, is utilized for this study. The experimental set-up is depicted in Figure 2.13, where the bathymetry includes a depth-varying shelf connected by two mild-slopes. Long waves from the deep water depth become steeper and possibly short-crested when propagating ontoashallowshelf; nonlinearityanddispersioneffectsmayneedtobetakenintoaccount for an accurate representation of the long wave transformation. For the offshore propagation region (Layers 1 and 2), the linear shallow water version ofCOMCOTisappliedusingarelativelycoarsegridsizeof1.5mand0.3m,respectively. Along the nearshore area, the nonlinear shallow water equations (Layer 3) are coupled with the Boussinesq model (Layer 4), both using the same relatively fine grid of 0.075 m. Notethatasecondsimulation,withoutusingtheBoussinesqmodel(i.e. Layer4removed), wasalsoperformed. Thiswillallowforadirectcomparisonbetweennearshorepredictions ofCOMCOTandtheBoussinesqmodel,withbothusingpreciselythesameincidentwave 44 Bathymetry (side view) h=4m 1 10 1 15 H=0.06m T=20 sec Shelf region 90m x=0m coastline 0.45m 1 200 10.81m 35.5m Grid setting* (plan view) Applied Model * Grid size ratio is not actually scaled LSW(Layer 1, x=1.5m, t=0.02sec) LSW(Layer 2, x=0.3m, t=0.01sec) NLSW(Layer 3, x=0.075m, t=0.005sec) Boussinesq (Layer 4, x=0.075m, t=0.005sec) x=90m x Figure 2.13: Definition sketch of a long wave propagating onto a shallow shelf (LSW:Linear Shallow Water Equation, NLSW: Nonlinear Shallow Water Equation) 45 condition and numerical grid sizes. Figure 2.14 presents time series comparisons of water surfaceelevationbetweenthemodelresultsandmeasurementsatdifferentlocations. The coupled model comparisons are given in the left half of the figure, while the COMCOT- onlyresultsareshownintherighthalf. Clearlydifferentbehaviorofthewavefrontinthe shallow shelf region is predicted by the coupled and COMCOT-only models. The front of the long wave becomes short-crested and generates (or strictly speaking, disintegrates into)severalsolitonsofdifferentsize(e.g.,MadsenandMei(1969));thisprocessisreferred to as tsunami wave fission in the literature. This transformation is the classic undular bore formulation which is dispersive in nature, and thus not predictable by the shallow water wave equations solved by COMCOT. The coupled model predicts a maximum sea surface elevation at the front of the tsunami which is 2.0 times larger than COMCOT alone, yielding a good agreement with measured data in both amplitude and speed. This type of difference is highly local in nature, and provides a reasonable picture of the magnitude and scale of dispersion-driven physics during nearshore tsunami evolution. 2.8 2004 Sumatra Tsunami Simulation Asapracticaltest,thecoupledmodelisappliedtothehistoricaltsunamieventofDecem- ber 2004 in the Indian ocean, through which the comprehensive lifespan of a tsunami, from its generation, propagation, shoaling, and run-up, might be investigated in true scale. Our specific geographic focus is Port Salalah, along the southeastern Omani coast- line. Asnotedinthemodeldescription,theBoussinesqmodelenhancedwiththeviscosity and vorticity terms (i.e. Kim et al, 2009) is capable of simulating turbulence effects such 46 10 20 30 40 50 60 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=80m Numerical(Boussinesq) Measured 10 20 30 40 50 60 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=80m Numerical(COMCOT) Measured 40 50 60 70 80 90 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=40m 40 50 60 70 80 90 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=40m 40 50 60 70 80 90 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=30m 40 50 60 70 80 90 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=30m 50 60 70 80 90 100 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=20m 50 60 70 80 90 100 −0.05 0 0.05 0.1 0.15 time(sec) ζ(m) x=20m Figure 2.14: Time histories of water surface elevations at 4 different locations. In the left halfofthefigureareshownthecoupledmodelresults(solidline: Boussinesqmodel(Layer 4, dashed line: Measured data), while the right half of the figure shows the COMCOT only results at the same times (solid line: Nonlinear Shallow Water Equation model (Layer 3), dashed line: Measured data) 47 Figure 2.15: Initial surface elevation of 2004 Sumatra tsunami as large eddies and wakes generated in the nearshore or harbors. Hence, another point of interest, aside from the dispersion differences noted in the previous comparison, is pre- dicting such complex turbulent physics, and this is one of the major advantages arising from integrating the two models. 2.8.1 Simulation Setup The bathymetric data and grid system has been organized to simulate 2004 Sumatra tsunami using the multigrid system in COMCOT. Figure 2.15 shows the initial sea sur- face elevation induced by the 2004 Sumatra earthquake within the entire computational domain, covering from (45 ◦ E, −10 ◦ N) to (105 ◦ E, 30 ◦ N). Open ocean bathymetry and 48 topography is taken from the GEBCO database, while shallow bathymetry off the coast of Oman is taken from digitized nautical charts. The parent domain, numbered Layer 1, has 4 subdomains given as Layers 2 through 5; higher numbered grids are nested within lower numbered grids. Layer 3, for example, is nested within Layer 2, and has a finer grid size and smaller time step. All the parameters necessary for simulation are listed in Table 2.2. From Layer 1 to 4, COMCOT is applied. Near Port Salalah in Oman (Layer 5), the main area of interest, the Boussinesq model with the highest grid resolution is appliedinordertocapturelocalandturbulentdynamics, suchasvortices, insidethehar- bor. It is worthwhile to note that the grid resolution is decreased by a factor of 400 from Layer 1 to Layer 5. Additionally, a readily calculated maximum stability index value, γ, is 0.8, assuming that the nonlinearity (ǫ) and water depth are the conservative values of 1 and 23m, respectively. While this stability index is near the limit of our recommended range, the use of the conservative values implies expected stability for this numerical configuration. For the generation of an initial surface condition for the tsunami, the three-subfault- source condition of Wang and Liu(2006) has been applied; the parameters of which are listed in Table 2.3. The runtime of the simulation is set to 830 minutes of physical time. Additionally, for the purpose of comparison between the results with and without coupling, another simulation using only COMCOT has been implemented with the same configuration except the grid size of Layer 4. In this COMCOT-only simulation, the grid resolution of layer 4 is set to 9.3 m, which is the resolution used by the Boussinesq layer in the coupled model. 49 Table 2.2: Grid setup for 2004 Sumatra Tsunami simulation Layer No. x range(longitude, ◦ E) y range(latitude, ◦ N) nx×ny dx dt(sec) Model Layer 1 45 ◦ 00 ′ 00 ′′ ∼105 ◦ 00 ′ 00 ′′ −10 ◦ 00 ′ 00 ′′ ∼30 ◦ 00 ′ 00 ′′ 1801×1201 2 ′ 1 LSW(S) a Layer 2 45 ◦ 17 ′ 12 ′′ ∼61 ◦ 38 ′ 48 ′′ 12 ◦ 17 ′ 12 ′′ ∼28 ◦ 18 ′ 48 ′′ 2455×2405 24 ′′ 0.5 LSW(S) Layer 3 52 ◦ 56 ′ 38 ′′ ∼56 ◦ 28 ′ 58 ′′ 16 ◦ 08 ′ 38 ′′ ∼18 ◦ 04 ′ 58 ′′ 2655×1455 4.8 ′′ 0.25 LSW(S) Layer 4 53 ◦ 58 ′ 00 ′′ ∼54 ◦ 04 ′ 00 ′′ 16 ◦ 54 ′ 00 ′′ ∼17 ◦ 00 ′ 00 ′′ 600×600 18.5m 0.125 NLSW b Layer 5 53 ◦ 59 ′ 46 ′′ ∼54 ◦ 01 ′ 18 ′′ 16 ◦ 55 ′ 46 ′′ ∼16 ◦ 57 ′ 19 ′′ 309×309 9.3m 0.125 BOUSS c a Linear Shallow Water Model in Spherical Coordinates b Non-Linear Shallow Water Model in Cartesian Coordinates c Boussinesq Model in Cartesian Coordinates 50 Table 2.3: Fault parameters for 2004 Sumatra earthquake Parameter Fault 1 Fault 2 Fault 3 Latitude of epicenter( ◦ N) 7.6 4.15 11.85 Longitude of epicenter( ◦ E) 93 94.55 92.3 Focal depth(km) 5 5 5 Length of fault plane(km) 670 200 300 Width of fault plane(km) 150 150 150 Dislocation(m) 15 15 15 Strike angle( ◦ ) 345 300 365 Slip angle( ◦ ) 90 90 90 Dip angle( ◦ ) 13 13 13 2.8.2 Results and Discussion Fromthesimulationresult,thefirsttsunamiwavesarriveatSalalahportinOman(Layer 5) approximately 420 minutes after the earthquake. This is comparable to the initial ar- rival time at the Port of 433 minutes post-earthquake, as reported by Okal et al. (2006). Afterward, successive attacks by a long train of tsunami waves caused a significant dis- turbance inside and immediately near the harbor. These disturbances are given in Okal et al.(2006) who discuss various ship incidents during the tsunami attack. They reported that a freighter docked at the berth had broken its mooring lines and drifted in- and outside the harbor, caught in a complex system of eddies and currents. The coupled model system appears to be able to represent these chaotic dynamics. Figure 2.16 shows one result of sea surface in each layer at 730minutes after tsunami generation, with an extra plot for the vorticity inside the harbor taken from the Boussinesq layer. Along the breakwater and coastline, the tsunami generates eddies of various sizes, and the flow is chaotic. Vorticity evolution as the tsunami propagates into the harbor has been reason- ably captured by the Boussinesq model and is depicted in Figure 2.17. 51 Figure 2.16: Surface elevation (m) in all layers at time=730min. The lower row shows the output from Layer 5, the Boussinesq model, of the free surface elevation (left) and vorticity (1/s) (right). 52 Figure 2.17: Vorticity (1/s) evolution inside Oman Salalah harbor at nine different times 53 A comparison between the vorticity results of COMCOT and the Boussinesq model is shown in Figure 2.18. Eddies are very weakly generated in COMCOT relative to the Boussinesq model. While bottom drag, which generates the boundary shear layers that curl up into the large eddies, are modeled differently in the two models, this is likely not the reason for the large difference; for a given velocity the bottom friction from the two formulations will be similar at this geophysical scale. The most likely cause of this largedifferenceisthenumericaltruncationerroroftheupwinddifferencinginCOMCOT, given as variable Er in section 4.2. This error can be expressed as Er =0.5(1−Cr)uΔx ∂ 2 u ∂x 2 =ν num ∂ 2 u ∂x 2 (2.24) and so can clearly be viewed as a diffusion term. For aCr≈0.5 and flow speeds ranging from 1-5 m/s, the numerical eddy viscosity, ν num , varies from≈ 2−10 m 2 /s. This is a verylargediffusioncoefficient,andtakenwiththeexpectationthatthevelocitycurvature is large inside boundary shear layers and eddies, it is evident that the numerical diffusion in the COMCOT model is driving the large differences in the vorticity patterns. As the two models predict different eddy patterns, the velocity predictions inside the harbor will be equally varied. Specifically, with the large and interacting eddies predicted by the coupled model, the simulated velocities are much larger. This kinematic aspect is of great importance to harbors during tsunamis, as it is the currents that lead to drag forces great enough to snap mooring lines, and transport large freighters as randomly meandering “ghost” ships. If one is interested in simulating this rotational features, numerical errors from low-order upwinding should be avoided. 54 Figure 2.18: Comparison of vorticity (1/s) evolution by Boussinesq-coupled model (left) and COMCOT-only (right) 55 2.9 Summary For the purpose of seamlessly modeling tsunami evolution from generation to inundation withfinescaleresolution,withoutthelossofimportantphysics,atwo-waycoupledmodel for tsunami simulation has been developed. The two components are the shallow-water solverCOMCOTandadispersive,turbulent,androtationalBoussinesqmodel. Ageneral framework in which the coupled model is implemented is as follows: Since COMCOT is well designed for generation and propagation of a tsunami in the open ocean, it will be responsible for the computation of oceanic evolution. On the other hand, the final stage of tsunami life, including nearshore dynamics such as inundation, nonlinear wave inter- actions, steep bore fronts, and turbulent activity, will be described by a fully nonlinear Boussinesq-type model. The Boussinesq model that can describe nearshore evolution of a tsunami with high physical detail is designed to be located flexibly within COMCOT as a nested layer. As coupling of two heterogeneous models may result in undesirable errors, a general benchmark test has been completed with various conditions provided for validity of the coupled model application. With regard to stability and accuracy, the simulation output is evaluated and general guidance for the coupled models application space has been presented; the so-called stability indexγ should be less than 0.9. As a further validation of the coupled model in the nearshore region, long wave propagation onto a shallow shelf has been examined and compared with laboratory data. Distinct dispersive effects at the leading wave front have been observed through use of the Boussinesq model, which demonstrates that near coastal areas, dispersive effects may be locally important. 56 Finally, a recent tsunami event, the 2004 Sumatra tsunami, has been simulated with a far field focus on the Port of Salalah. The coupled model has successfully simulated various sizes of eddies generated by the tsunami through turbulence activity. The results are further supported by observations addressed in Okal et al.(2006). It is found that one needs to be very careful when using numerical solution schemes with leading order diffusion errors to predict such rotational features, as this numerical error can rapidly remove intense shear layers and strong eddies from the current field. 57 Chapter 3 Interaction of Shallow Water Waves with Weakly Sheared Currents of Arbitrary Profile A set of depth-integrated equations describing combined wave-current flows is derived in this chapter. To account for the turbulent effects by nonlinear interaction between waves and underlying currents with arbitrary horizontal vorticity, additional stresses are introduced. A parameter b is included in the additional stresses to represent the radiation stress of waves over the ambient current field. Doppler shift effect is theoret- ically proved to be retained appropriately in the equation set. To solve the equations, a fourth-order MUSCL-TVD scheme incorporated with approximated Riemann solver, is adopted for leading order terms while a cell-averaged finite volume method is utilized for higher-order terms. The model results are validated through comparisons with three experimental data sets. The first simulation is about the propagation of a long wave over a depth-uniform or linearly-sheared current and the comparison with the measurements shows good agreements. Then, based on Kemp and Simons(1982, 1983)’s experiments, a reasonable value for b is estimated. Finally, the model is applied to a more complex configuration where bichromatic waves are interacting with a spatially varying current. 58 The simulated results clearly indicate that the model is capable of predicting the wave and current flows more accurately. 3.1 Introduction The coastal area involves hydrodynamic situations in which different flows such as waves, currents and tides are coexisting and interacting with each other. In the shallow coastal zones, inparticular, agravitywaveissubjecttointeractwiththecurrentwhichisusually driven by tides, thermocline, salinity variations and river mouth discharges. Morphody- namic changes, mixing and transport of solutes in coastal region will be governed by the combined wave-current flows. Therefore, the interaction between wave and current is of prime importance to the physical processes in the coastal zone. Waves and currents have different hydrodynamic properties which will allow the dis- tinctions between them. For instance, waves transfer energy without conveying mass while the currents transfer both energy and mass. Also, they have effects of different timescales on the coastal environments. Sediment transport at nearshore area is driven mainly by currents, and therefore evolves on the long-term scales (commonly on the or- der of days). On the other hand, waves contribute to the sediment transportation on shorter-time scale by the radiation stress. Therefore, both long-time scale and short-time scale effects have to be included while studying the physical processes in wave-current co-existing situations. 59 In wave and current environments, the nonlinear interaction plays a significant role andwecannotsimplysuperposethetwocomponents. Furthermore,thenonlinearinterac- tionshavenon-negligibleimpactsonhydrodynamicsofwave-currentsystem, especiallyin the turbulent boundary layer with high roughness(e.g., Grant and Madsen(1979), Davies et al.(1988)). Therefore, turbulent processes as well as velocity profiles are of particular interest in the wave-current situation. Turbulent interaction of wave and current has been investigated quite extensively in the past decade, theoretically and experimentally. Laboratory experiments includes Thomas(1981, 1990),Kemp and Simons(1982, 1983), Klopman(1994), Swan et al.(2001). Most recently, Fernando et al.(2011) presented experimental results on wave-current in- teraction at an angle along with comprehensive reviews on previous works. One inter- esting finding in these experiments is that when the waves and currents are co-flowing (or counter-flowing), mean velocity near the free surface tends to curl back (or forward), giving negative (or positive) gradient of velocity (∂u/∂z). Smith(2006) mentioned that the backward flows are introduced compensate for the Stokes drift effect by waves. A possible explanation for this mixing-like process induced by waves is that the presence of waves introduces additional shear stress over the underlying mean flow, yielding a modi- fiedReynoldsstressesinthecombinedwave-currentfield(e.g.,You(1996),Groenewegand Klopman(1998), Huang and Mei(2003), Umeyama(2005), Yang et al.(2006), Lin(2008)). Umeyama(2005)providedafeasibleinsightregardingthesepreviously-mentionedphysics throughlaboratoryinvestigations. Theyconductedseriesofexperimentsonwave-current turbulent intensities and Reynolds stress in combined wave-current flows. Experiments 60 showed the modification of Reynolds stress (i.e., u ′ v ′ ) by the action of waves on mean flows. On the other hand, numerical models for wave-current interaction have been devel- oped for deep or for finite depth waves(e.g., Swan et al.(2001), Swan and James(2001), Nwogu(2009)) and for long waves(e.g., Benjamin(1962), Freeman and Johnson(1970), Shen(2001)). Grant and Madsen(1979) suggested different eddy viscosities in- and out- side bottom boundary layers and some works stem from this(e.g., Christofferson and Jonsson(1985), Davies et al.(1988)). Majority of the turbulent wave-current model extensively examined near-bed physics (e.g., Kim et al.(2001)) rather than covering entire water column. So far, some numerical models have shown the capability to recreate aforementioned turbulence mixing process throughoutthedepth, inducedbywavesoncurrents(e.g., Dingemansetal.(1996), Groen- weweg and Klopman(1998), Olabarrieta et al.(2010)). In addition, some equations have been proposed to describe the velocity profile of the mean flow either empirically (e.g., You(1996)) or analytically(e.g., Huang and Mei(2003), Yang et al.(2006)). Interestingly, most of the previously-mentioned equations have added the higher or- der correction terms to the leading order solution. Despite of the elaborative efforts to understand these physics, little is developed to account numerically for turbulent wave- currenteffectinshallowwaterflows. ModelssuchasFUNWAVE(e.g., Kirbyetal.(2003)) are capable of predicting wave-current interactions, however only few of these are pos- sibly able to recreate such a non-explainable phenomena. Therefore, including higher order correction effects will enhance the capability of the hydrodynamic model to study wave-current flow field. 61 In this chapter, we derived a set of depth-integrated equations describing combined wave-current flows. The effects of turbulence introduced not only by bottom friction but also by nonlinear interaction between the waves and the underlying currents with arbitrary horizontal vorticity, are included. This chapter is organized as follows. Firstly, wave-current interactions will be briefly exploredinfollowingsection. Nextsectionwillbedevotedtothedescriptionofderivation ofthegoverningequations. Then,validityofthecurrentmodelwillbeperformedthrough three numerical simulations. Final conclusion is made at the end of this chapter. 3.2 Brief Review on Wave-Current Interactions Waves and currents are known to influence each other when existing together. For ex- ample, waves propagating on an opposing current undergoes some transformations as they become shorter and steeper. In the case where the opposite current speed exceeds the wave group speed, the waves are more likely to break and thereby provide another source of turbulence, compared to a wave on calm water(e.g., Yao and Wu(2005)). These changes will hydrodynamically affect the current as it will be briefly illustrated below. 3.2.1 Waves over the Current Wavestravelingoverunderlyingcurrentsexperienceamodulationintheirkinematicsand dynamics such as a change in wave number, frequency, and wave height. Waves become steeper and higher on following currents whereas the opposite for opposing currents. Moreover, underlying currents effectively control the wave frequency in such a way that the wave period will be longer over the following currents and shorter over the opposing 62 one. Doppler shift is a quite common concept to explain such modulations in dispersion relationship. For uniform background currents, this effect can be mathematically proved and expressed as, σ 2 =(ω−ku C ) 2 =gktanhkh (3.1) where σ is intrinsic(or relative) angular frequency, ω the apparent(or absolute) angular frequency,k the wave length,u C the current speed andh water depth. For linearly shear current, it was found to be, σ 2 = ω−ku Cs + W 0 2 tanhkh 2 =gktanhkh (3.2) where u Cs is current velocity at free surface andW 0 is current’s constant vorticity. Anothereffectofthecurrentonthewaveistherefractionduetothecurrentvariation in space, which is much similar to that by the bathymetric changes(Lin (2008)). 3.2.2 Currents under the Waves Current fields also tend to be deformed by wave actions. The wave riding on the current has many potential factors which may affect the mean flow field(e.g., a radiation stress or bottom friction enhancement). The mechanism by which the wave change the currents is yet unclear; however, it is believed that an additional shear stress exerts on the mean flow owing to nonlinear wave-current interactions, feeding horizontal vorticity to the interior flow. This stress is considered to be generated by the radiation stress of waves and mean flows(Lin(2008)). Many experimental results showed consistently that this 63 stress is capable to tilt forward or backward the current velocity depending on the wave- current directions, and is maximized near the free surface with decreasing into the depth. The vorticity of the current is also known to affect this mechanism(Swan et al.(2001)). However,near-bedvelocityisnotsignificantlyaffectedbysuchstress,soitcanbemodeled using the general ‘log-law’ profile(Kemp and Simons (1982), Fernando et al. (2011)). Although the above effect is relatively concentrated in the upper layer of the fluid, this should be included for the complete description of velocity profile for accurate modeling. Some equations have been proposed to identify the resultant velocity profile of mean flow under wave action. You(1996), for instance, suggested semi-empirical equation based on experimental record as, u(z) u ∗ = 1 κ log z+h δ +C h κu ∗ |u ∗ | log −z h (3.3) where u(z) is mean flow velocity profile, u ∗ is shear velocity, z is vertical axis directing upward from free surface,κ is von Karman constant,h is water depth andδ is roughness height. Dimensional parameterC can be obtained by empirical formula from You(1996), but with some lack of consideration on waves properties. The second term of Equa- tion 3.3 implies higher-order correction component to the first term(i.e., log-law profile which has its origin in Prandtl’s mixing length hypothesis) under the wave-free situation. Similar formulations can be found in Umeyama(2005) and Yang et al.(2006) with minor differences. 64 Typically,bedroughnessissupposedtobeenhancedbywave-currentco-existenceand this idea can be realized theoretically by considering the bottom friction for combined wave-current flows(e.g., see Grant and Madsen(1986)). 3.3 Boussinesq Equations for Waves and Currents Asetofdepth-integratedequationsforlongwavesunderthecurrentfieldsisderivedinthis section. The perturbation approach to manipulate a primitive equation into a derivative one is used in the present study. This technique is adopted to develop Boussinesq-type equations including the effect of wave-current interactions. As are ‘standard’ Boussi- nesq equations, nonlinearity and dispersive effects of long waves would be the basis in the perturbation procedure. The viscous terms can be added as perturbed terms into the inviscid Boussinesq equations to explain bottom-induced turbulence effects(Kim et al.(2009)). The primary procedure to introduce a turbulence induced by wave-current interaction, follows Kim el al.(2009). 3.3.1 Non-dimensionalizedGoverningPhysicsandBoundaryConditions PhysicalvariablesaredefinedasshowninFigure3.1todescribethepropagationofwaves over depth-varying currents. The variables are normalized by characteristic variables introduced below. Typical length scale ℓ 0 and h 0 are used for horizontal and vertical 65 coordinates, respectively. Using this set of scaling parameters, nondimensional variables are obtained as, (x,y)= (x ∗ ,y ∗ ) ℓ 0 , z = z ∗ h 0 , ζ = ζ ∗ h 0 , h= h ∗ h 0 , t= t ∗ √ gh 0 ℓ 0 , (u,v)= (u ∗ ,v ∗ ) √ gh 0 , w = w ∗ µ √ gh 0 , p= p ∗ ρgh 0 , µ = h 0 ℓ 0 , ν h t = ν h ∗ t αh 0 √ gh 0 , ν v t = ν v ∗ t βh 0 √ gh 0 (3.4) where () ∗ denotes dimensional variable, (u,v) and w represents velocity components in (x,y) and z directions, respectively, ζ is surface elevation, t is time, p is pressure, g is gravitational acceleration andρ is density. ν h t andν v t represent turbulent eddy viscosities in horizontal and vertical directions, respectively. Parameter µ is chosen to scale the dispersion property of long waves. It is noted that no parameter for nonlinearity is gained above, so that leads to fully nonlinear equations. In addition, small parameters α and β are used for the viscous terms, whose full descriptions can be found in Kim et al.(2009). While the procedures throughout the entire derivation fairly resemble Kim et al.(2009), velocityu represents the instantaneous horizontal velocity including both wave and current components. Therefore the velocityu means combined wave-current velocity unless otherwise stated. Within these parametric frameworks, continuity and Navier-Stokes equation can be casted into normalized versions as follows. ∇·u+w z =0 (3.5) 66 uC * h * x * y * z * z * z * = 0 z * = -h * Figure 3.1: Definition sketch of long wave propagation over underlying currents of arbi- trary profile u t +u·∇u+wu z +∇p=αµ ∇· ν h t ∇u + β µ (ν v t u z ) z (3.6) µ 2 w t +µ 2 u·∇w+µ 2 ww z +p z +1=αµ 3 ∇· ν h t ∇w +βµ (ν v t w z ) z (3.7) where ∇ = (∂/∂x,∂/∂y) is differential operator in horizontal plane and, u = (u,v) is horizontal velocity vector of combined wave-current flows. Subscript z and t mean derivative operators in vertical coordinate and time, respectively. Conditions applied at free surface and at bottom boundaries are expressed in dimen- sionless form, as well. w =ζ t +u·∇ζ at z =ζ (3.8) 67 w+u·∇h=0 at z =−h (3.9) 3.3.2 Reynolds Stresses Under Combined Wave-Current Flow A definition of Reynolds stress under the wave-current condition is made in this section, priortoworkingonderivationsoftheequationset. AsbrieflyaddressedinChapter3.2,an underlying current influenced by the wave has a different velocity profile from that under wave-free condition. This difference is mainly owing to the modified Reynolds stress, which is suggested in Umeyama(2005) and Yang et al.(2006). These turbulence effects by nonlinear wave-current interaction can be modeled through the modified Reynolds stress. Based on Umeyama(2005) and Yang et al.(2006), the turbulent shear stress concerning both bed roughness and nonlinear wave-current interaction is described as, τ =τ b ζ−z ζ +h +τ bb z+h ζ +h (3.10) where τ = (τ x ,τ y ) is the shear stress, and τ b = (τ x b ,τ y b ) is bottom stress. τ bb = (b x τ x b ,b y τ y b ) is additionally defined for convenient expression. Dimensionless parameter b=(b x ,b y ) which is identical to that of Yang et al.(2006), is used to represent the linear distribution of turbulent shear stress by the wave-current interaction. The first term on the right hand side represents simply the stress due to the bottom friction, whereas 68 the second represents the stress induced by wave-current interaction. Therefore, ignor- ing the second component will result in abandonment of turbulent mixing through the wave-current interaction. Quadratic drag law is assumed to estimate the bottom stress as, τ b =ρf cw u|u| (3.11) where u is ‘combined’ wave-current velocity, ρ is a density of fluid. Friction factor f cw can be gained by thef cw =f/4, wheref is estimated from a formula by Haaland(1983). 3.3.3 Subgrid Scale Turbulent Closure Model In shallow water flows, the vertical eddy viscosity is assumed to be depth-constant (El- der(1959)) and this allows us to simplify the three dimensional turbulence in nature into two dimensional problem in horizontal plane. Hinterberger et al.(2007) developed an effi- cienttwodimensionalmodelinshallowwaterflowsandfoundthattwodimensionaldepth- averaged LES(Large Eddy Simulation) model produced sufficient accuracy for practical purposes compared to the three dimensional LES. In the present work, Smagorinsky’s model is utilized for horizontal subgrid eddy viscosity(ν h t ), while Elder(1959)’s model adopted for vertical one(ν v t ) in order to essentially capture the vortical features of wave- current flows without the loss of the flow details(see Kim et al.(2009)). 69 3.3.4 Depth-Integrated Momentum Equations for Waves and Currents The assumption is made that the turbulent viscosity effects are as weak, and of a similar order as dispersion. That is, O(µ 2 )=O(µβ )≪1 (3.12) Either of µ 2 and µβ can be utilized as a small parameter in the derivation, but the former has been chosen(Kim et al.(2009)). Each physical variable is thereby expanded as a power series of small parameter µ 2 as follows. f = ∞ X n=1 µ 2n f n , f =(p,u,v,w) with all f n =O(1) (3.13) Limiting Equation 3.7 to the lowest order yields, (p 0 ) z +1=0 (3.14) which refers to hydrostatic status. Since ∇p 0 is independent of vertical coordinates, we can notice allz-dependency will be eliminated from leading order terms in Equation 3.6, includingu(i.e.,u(x,y,z,t)=u 0 (x,y,t)+O(µ 2 )). Depth integration of continuity equation, Equation 3.5, with applying the boundary conditions, Equation 3.8 and 3.9, will produce, 70 w 0 =−z(∇·u 0 )−∇·(hu 0 )=−zS−T (3.15) Vertical velocity expressed in terms of horizontal velocity,u 0 , allows us to eliminate all the leading ordered vertical velocity in Equation 3.6 and 3.7, which means horizontal vorticity(Ω) is at least second order(i.e. Ω 0 =0). This can be seen as follows, as well. Ω ∗ =u ∗ z −∇w ∗ = µ 2 r g h 0 (u 1z −∇w 0 )+O(µ 4 ) = µ 2 r g h 0 Ω 1 +O(µ 4 ) (3.16) To construct the horizontal velocity structure, we start from horizontal vorticity. Depth integration ofΩ 1 from−h to z gives, u 1 = u 1 | z=−h − 1 2 z 2 ∇S+z∇T − 1 2 h 2 ∇S−h∇T + Z z −h Ω 1 dz ′ +O(µ 2 ) (3.17) Now, we can get the full description of horizontal velocity(u) by substituting Equa- tion 3.17 into Equation 3.13 as follows, 71 u=u 0 + µ 2 u 1 | z=−h − 1 2 z 2 ∇S−z∇T + 1 2 h 2 ∇S−h∇T (3.18) + Z z −h Ω 1 dz ′ +O(µ 4 ) Forarepresentativehorizontalvelocity,Nwogu(1993)’sdepth,z α =−0.531hischosen and this yields u α =u 0 + µ 2 u 1 | z=−h − 1 2 z 2 α ∇S−z α ∇T + 1 2 h 2 ∇S−h∇T (3.19) + Z zα −h Ω 1 dz ′ +O(µ 4 ) Then, subtracting Equation 3.19 from Equation 3.18 will produce the expression ofu in terms ofu α as, u=u α + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T + Z z zα Ω 1 dz ′ +O(µ 4 ) (3.20) Sinceuisthecombinedwave-currentvelocity,wecandecomposeitintoseparateones. The underlying current velocity(u C ) of arbitrary profile is assumed be constant in time, in the present study. Therefore, the horizontal vorticity of current will be of arbitrary shape and assumed to be ofO(µ 2 ), i.e., weakly sheared according to Equation 3.16. The horizontal velocity of combined wave-current field is expressed as, 72 u = u Wα +u Cα (3.21) + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T + Z z zα Ω 1 dz ′ + O(µ 4 ) whereu Wα andu Cα are the horizontal velocity vector at z =z α by waves and currents respectively, satisfyingu α =u Wα +u Cα . Note that we usedu 0 =u α +O(µ 2 ) here. For the simplicity of calculation, we will useu α instead of the decomposed ones. Velocity of O(µ 2 ) can also be interpreted as rotational and irrotational parts. From Equation 3.16, R z zα Ω 1 dz ′ is the rotational velocity which is generated either through the turbulentshearstress(τ)orthroughtheintrinsicallyshearedcurrent(u C )itself. Therefore it follows that τ +ρν v t (u C ) z =ρν v t Z z zα Ω 1 dz ′ z =ρν v t Ω 1 (3.22) which also leads to Ω 1 = τ b ρν v t ζ−z ζ +h + τ bb ρν v t z+h ζ +h +(u C ) z (3.23) usingEquation3.10. TherotationalvelocitycomponentofO(µ 2 )willbeobtainedthrough depth-integration ofΩ 1 . 73 Z z zα Ω 1 dz ′ = τ b ρν v t 1 ζ +h 1 2 z 2 α −z 2 +ζ(z−z α ) + τ bb ρν v t 1 ζ +h − 1 2 z 2 α −z 2 +h(z−z α ) + (u C −u Cα ) (3.24) By inserting Equation 3.24 into Equation 3.21, final structure of horizontal velocity can be written up to O(µ 2 ), u=u α + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T (3.25) + µ 2 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) + µ 2 Ψ b − 1 2 z 2 α −z 2 +h(z−z α ) + µ 2 (u C −u Cα )+O(µ 4 ) in which Ψ =τ b /{ρν v t (ζ +h)} and Ψ b =τ bb /{ρν v t (ζ +h)} are defined. Equation 3.25 revealsthatwehavethreesourcesofhorizontalvorticity,twoofwhicharebottomfriction and wave-current interaction while the remainder is the current’s external shear. Thehorizontalvelocityprofilecanbeutilizedforderivingdepth-integratedmomentum equations. The remaining steps are somewhat tedious, so detailed procedures to derive a depth-integrated equation are included in Appendix D. Asaresult,thedepth-integratedmomentumequationincludingtheeffectsofbottom- induced turbulence and wave-current interaction, can be written as, 74 (u α ) t + u α ·∇u α +∇ζ +µ 2 H p +H t +H c +ξ p +ξ t +ξ c − αµ ∇· ν h t ∇u α +βµν v t ∇S+βµ τ b ρ(ζ +h) −βµ τ bb ρ(ζ +h) − βµν v t {(u C ) z | z=ζ −(u C ) z | z=−h } = O(µ 4 ,αµ 3 ,βµ 3 ) (3.26) where H t = (ζ−h) 2 (Ψζ) t − ζ 2 −ζh+h 2 6 Ψ t + Ψ z 2 α 2 −ζz α t + (ζ−h) 2 ∇{u α ·(Ψζ)}− ζ 2 −ζh+h 2 6 ∇(u α ·Ψ) + ∇ u α · Ψ 1 2 z 2 α −ζz α −Ψ ( ζ 2 −ζh−2h 2 S 6 + (ζ +h)T 2 ) − (Ψ b ) t ( 1 2 z 2 α +hz α − ζ 2 +2ζh−2h 2 6 ) − ( 1 2 z 2 α − ζ 2 −ζh+h 2 6 ) ∇(u α ·Ψ b ) + ∇{u α ·(Ψ b h)} 1 2 (ζ−h)−z α − Ψ b ( 2ζ 2 +ζh−h 2 S 6 + (ζ +h)T 2 ) (3.27) 75 H p = 1 2 z 2 α ∇S t +z α ∇T t − 1 2 ∇ ζ 2 S t −∇(ζT t )+T∇T + 1 2 ∇ z 2 α u α ·∇S +∇(z α u α ·∇T)+ 1 2 ∇ ζ 2 S 2 − 1 2 ∇ ζ 2 u α ·∇S −∇(ζu α ·∇T)+∇(ζTS) (3.28) H c = 1 ζ +h ∇ u Wα · Z ζ −h u C dz −∇(u Wα ·u Cα ) − ζS+T ζ +h u C | z=ζ − hS−T ζ +h u C | z=−h + S ζ +h Z ζ −h u C dz (3.29) Alsoξ p = ξ p x ,ξ p y ,ξ t = ξ t x ,ξ t y andξ c = ξ c x ,ξ c y are described as ξ p x = −v α n (z α ) x (z α S y +T y )−(z α ) y (z α S x +T x ) o (3.30) − n (v α ) x −(u α ) y o "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) S y + z α − (ζ−h) 2 T y # ξ p y = u α n (z α ) x (z α S y +T y )−(z α ) y (z α S x +T x ) o (3.31) + n (v α ) x −(u α ) y o "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) S x + z α − (ζ−h) 2 T x # 76 ξ t x = −v α " ψ y 1 2 z 2 α −z α ζ x − ζ 2 −ζh+h 2 6 (ψ y ) x + (ζ−h) 2 (ψ y ζ) x − ψ x 1 2 z 2 α −z α ζ y + ζ 2 −ζh+h 2 6 (ψ x ) y − (ζ−h) 2 (ψ x ζ) y # − n (v α ) x −(u α ) y o ψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) + v α " ψ y b 1 2 z 2 α +z α h x + ζ 2 −ζh+h 2 6 ψ y b x + (ζ−h) 2 ψ y b h x − ψ x b 1 2 z 2 α +z α h y + ζ 2 −ζh+h 2 6 (ψ x b ) y − (ζ−h) 2 (ψ x b h) y # + n (v α ) x −(u α ) y o ψ y b ( z 2 α 2 +z α h− 2ζ 2 +2ζh−2h 2 6 ) (3.32) ξ t y = u α " ψ y 1 2 z 2 α −z α ζ x − ζ 2 −ζh+h 2 6 (ψ y ) x + (ζ−h) 2 (ψ y ζ) x − ψ x 1 2 z 2 α −z α ζ y + ζ 2 −ζh+h 2 6 (ψ x ) y − (ζ−h) 2 (ψ x ζ) y # + n (v α ) x −(u α ) y o ψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) − u α " ψ y b 1 2 z 2 α +z α h x + ζ 2 −ζh+h 2 6 ψ y b x + (ζ−h) 2 ψ y b h x − ψ x b 1 2 z 2 α +z α h y + ζ 2 −ζh+h 2 6 (ψ x b ) y − (ζ−h) 2 (ψ x b h) y # − n (v α ) x −(u α ) y o ψ x b ( z 2 α 2 +z α h− 2ζ 2 +2ζh−2h 2 6 ) (3.33) 77 ξ c x = v α ζ +h ( Z ζ −h u C dz y −ζ y u C | z=ζ −h y u C | z=−h − Z ζ −h v C dz x +ζ x v C | z=ζ +h x v C | z=−h − 1 ζ +h n (v α ) x −(u α ) y o Z ζ −h v C dz (3.34) ξ c y = u α ζ +h Z ζ −h v C dz x −ζ x v C | z=ζ −h x v C | z=−h − Z ζ −h u C dz y +ζ y u C | z=ζ +h y u C | z=−h ) + 1 ζ +h n (v α ) x −(u α ) y o Z ζ −h u C dz (3.35) Note that u α = (u α ,v α ), Ψ = (ψ x ,ψ y ) and Ψ b = ψ x b ,ψ y b . In Equation 3.34 and 3.35 Leibnitz rule had been applied, as well. 3.3.5 Depth-Integrated Continuity Equation for Waves and Currents SinceEquation3.5doesnotincludenonlinearterms,depthintegratedcontinuityequation can be obtained by fairly simple procedure. Integration of Equation 3.5 over the water column with applications of boundary conditions (Equation 3.8 and 3.9) gives ∇· Z ζ −h udz+ζ t =0 (3.36) 78 where we used Leibnitz’s rule for differentiating of integral terms. Substituting Equation 3.25 into above equation will approximate this equation into depth-integrated one, yielding ζ t +∇{(ζ +h)u α }+µ 2 N p +N t +N c =O µ 4 (3.37) where N p =−∇· (ζ +h) ζ 2 −ζh+h 2 6 − z 2 α 2 ∇S+ ζ−h 2 −z α ∇T (3.38) N t = ∇· " Ψ(ζ +h) ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# − ∇· " Ψ b (ζ +h) ( z 2 α 2 +z α h− 2ζ 2 +2ζh−2h 2 6 )# (3.39) N c =∇· Z ζ −h u C dz−(ζ +h)u Cα (3.40) 3.3.6 Modulation of Dispersion Properties by Currents Doppler shift is the effect of current on wave which should be included in wave-current models,andmaybeverifiedanalytically(Chenetal.(1998)). Inthissection,Dopplershift 79 will be investigated through the comparison of dispersion property embedded in the gov- erning equations, Equation 3.25 ,with the ‘analytical’ one. Assuming the steady current field and the linear wave system, the dispersion relationship can be derived analytically from the continuity Equation 3.37 and the momentum Equation 3.26. Subsequently, the linearized one dimensional depth-integrated governing equations in a constant depth under uniform current(u C =u Cα ) are written as, ζ t +h(u Wα ) x +u C ζ x +h 3 1 3 +γ (u Wα ) xxx =0 (3.41) (u Wα ) t +u C (u Wα ) x +gζ x + h 2 γ{(u Wα ) xxt +u C (u Wα ) xxx }=0 (3.42) where γ = (1/2)(z α /h) 2 +z α /h. All dimensions are recovered for the physical variables without using ∗ . To obtain the linearized dispersion characteristics, we imagine sinusoidal waves prop- agating with an angular frequency ω and wave number k, whose surface elevation and velocity are defined as, ζ =a m exp{i(kx−ωt)} (3.43) 80 u Wα =u m exp{i(kx−ωt)} (3.44) where a m is wave amplitude and u m is the maximum velocity. Inserting Equation 3.43 and 3.44 into Equations 3.41 and 3.42 presents linear dispersion relationship as follows. σ 2 =(ω−u C k) 2 =ghk 2 1− 1 3 +γ (kh) 2 1−γ(kh) 2 (3.45) AsmentionedbyNwogu(1993),thevalueofγ(=−2/5)wouldexactlyconfirmDoppler effectwhencomparedtoPade(2,2)approximantoflineardispersionrelationoffirst-order Stoke’s theory, while some errors have been found in Chen et al.(1998) in some other Boussinesq models. For long waves in a range of kh<π, above approximation will give overallgoodagreementtoexactlineardispersionwithanerrorlessthan5%,butthevalue of -0.39 forγ has been used for extended range of applicablekh and better accuracy. As a result, it is demonstrated that correct Doppler shift property retains in the equations. 3.3.7 Free parameter (b x ,b y ) Reasonable estimation of free parameter (b x ,b y ) is not yet achieved in spite of a com- prehensive effort by some researchers(e.g., Yang et al.(2006)). In this part, reasonable estimation for b x will be attempted as well as b y . Following the procedures in You(1996), with the use of Reynolds stress defined in Equation 3.10, an empirical formula for assessment of b x can be drawn out as, 81 (for following currents) b x = (Hω)kh u ∗ |u ∗ |sinhkh 0.004 Hω sinhkh −0.189|u ∗ | (3.46) (for opposing currents) b x = (Hω)kh u ∗ |u ∗ |sinhkh 0.002 Hω sinhkh −0.111|u ∗ | (3.47) Inthisestimation,theparameterb x hasaquadraticfunctionofwaveheightwithnon- zero value of H at vertex. By intuition, this may not be reasonable since the magnitude ofb x would not be proportional to the wave height and the application of this estimation resulted in poor agreements. Moreover, for phase-resolving model, direct application of this estimation would be difficult unless the local wave characteristics(i.e., H,k,ω) are known a prior. Therefore, an attempt to find an alternative way assessing b x is made here. The most two important factors governing wave-current interaction are considered as radiation stress of waves and bed shear stress by currents. Radiation stress S xx (or momentum flux due to the the presence of waves) under the underlying currents are clearly defined in Jonsson et al.(1978), and using the pressure and velocity expres- sions(Equation D.2 and 3.25), it will be given as, 82 S xx = Z ζ −h p+ρu 2 dz − 1 2 ρg(h+hζi) 2 −ρ Z ζ −h u 2 C dz = ρ 6 (h+ζ) 2 {3g−3(T t +u α ·∇T −TS) +(S t +u α ·∇S−S 2 )(h−2ζ) + ρ 60 15Y 2 (ζ +h)−30Y(T x )(ζ 2 −h 2 ) +10 2(T x ) 2 −(S x )Y (ζ 3 +h 3 )+15(S x )(T x )(ζ 4 −h 4 ) +3(S x ) 2 (ζ 5 +h 5 ) − 1 2 ρg(h+hζi) 2 −ρ Z ζ −h u 2 C dz (3.48) which is ended up with well-known form when applying Airy theory in the absence of background currents. S xx = 1 8 ρgH 2 2kh sinh2kh + 1 2 (3.49) In Equation 3.48, <> denotes phase-averaging and Y = 2u α +z 2 α S x + 2z α T x . As a reminder, S =∇·(u Wα +u Cα ),T =∇·{h(u Wα +u Cα )} and the velocity used here is ‘potential’partofEquation3.25. Radiationstressiny direction(stressnormaltoy-plane) can be obtained with the same procedure. Here, Equation 3.48 would be preferred in phase-resolving model since Equation 3.46 and 3.47) may not be appropriate in the time-dependent frame such as Boussinesq-type equations. Radiation stress scaled on currents stress will represent relative strength of 83 wave-average effects over the time-mean flow. Therefore, the linear relation between parameter b x and relative radiation stress of waves can be assumed as b x =C b S xx τ x Cb h (3.50) whereτ x Cb is bottom stress by currents. The value ofC b will be assessed in Chapter 3.4.2. 3.4 Validation In this section, three numerical tests will be performed to validate the present model. These numerical tests include the propagation of regular waves over uniform or depth- varying currents in a constant depth, and bichromatic waves traveling over spatially varying current field. 3.4.1 Waves over Uniform or Linearly-varying Currents Inthefirsttest,wavesdeformedbyuniformandshearedcurrentsareexamined. Swan(1990) investigated a modification of the wave motion due to the current effects, in particular the sheared currents. Four different cases within the relatively shallow or intermediate depth regime (kh ≤ 3.0) are chosen. Two of them include following currents while the rest include the opposing currents. The experimental parameters of waves and currents are listed in table 3.1. In CASE 1F and 1A, the current fields are essentially uniform throughout the water depth even though some shear may exist near the bottom. On the other hand, strong linear shear has been developed by honeycomb in CASE 2F and 2A, 84 Table 3.1: Wave and Current Conditions in Swan(1990)’s Experiments Case Wave height Water depth Wave period kh a Current Profile a m h T u C 1F 35.1mm 0.35m 1.412s 0.88 0.108m/s 1A 35.7mm 0.45m 0.877s 2.91 -0.120m/s 2F 31.5mm 0.35m 1.418s 0.85 following, linear shear b 2A 61.5mm 0.35m 1.420s 1.10 opposing, linear shear c a Calculated from linear dispersion relationship with currents, i.e., Equation 3.1 b See Figure 3.2 (a) c See Figure 3.2 (b) resulting in a current profiles shown in Figure 3.2. Such current conditions represent a low Froude number flow range (≈0.058∼0.091). Numerical tests are performed based on the actual conditions of experiments with a fixedgridsizeΔx=0.02mandavariabletimestepΔtaccordingtoCFLcondition(0.4in thepresenttest). Thecurrentvelocityispracticallyassumedtohavenobottomboundary layer as the channel bed is covered with plate glass in the experiment. This allows us to simulate this case in the model without bottom-induced turbulence. Therefore, the rotational effect is generated solely by the sheared current. Figure 3.3 through Figure 3.5 show the simulated results for each case along with the measured wave motions. The surface elevation and the oscillatory velocity in each figure have been measured at a location in the computational domain, at least 3 wave lengths farfromtheinternalsourceofwavegeneration, during2∼3waveperiodsasSwan(1990)’s attempt. The velocity profiles at different depth levels are calculated using Equation 3.25. In thedepth-uniformcurrentscases(Figure3.3and 3.4),thecomputedsurfaceelevationand the velocity at different depth levels are generally in a good agreement with experiments 85 −0.1 0 0.1 0.2 0.3 0.4 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 current(m/s) depth(m) (a) −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 current(m/s) depth(m) (b) Figure3.2: MeasuredshearcurrentsintheabsenceofwavesinSwan(1990); (a)CASE2F; (b) CASE2A in both magnitude and phase, even though some errors are found in Figure 3.4 (d) which are mostly owing to larger kh value(≈2.91). Computations agree well with the measurements for the case of opposing sheared currentasshowninFigure3.5. Thisdemonstratesthattherotationalbehaviorofcurrents needs to be taken into account when describing the velocity field in the combined flow of wave-current. The inclusion of the rotational terms due to the current’s vorticity in Equation 3.25 represents practically its effects on the flow field. Furthermore, the maximum values of oscillatory velocities given in Figure 3.6 are in good agreements with theStokes3rdordersolutionsbyKishidaandSobey(1988). Therefore, thissupportsthat the effects of sheared currents are appropriately included in the velocity profile as well as in the governing equations. As pointed out by Swan (1990), the conventional ‘potential’ 86 0 0.5 1 1.5 2 2.5 −0.05 0 0.05 ζ (m) (a) 0 0.5 1 1.5 2 2.5 −0.2 −0.1 0 0.1 0.2 u W (m/s) (b) z=−0.1m 0 0.5 1 1.5 2 2.5 −0.2 −0.1 0 0.1 0.2 u W (m/s) (c) z=−0.2m 0 0.5 1 1.5 2 2.5 −0.2 −0.1 0 0.1 0.2 t(sec) u W (m/s) (d) z=−0.3m Figure 3.3: Comparison of wave records under following, uniform current (Case 1F) between experiment(o) and numerical solution(-): surface elevation(a) and oscillating velocities at z=-0.1m(b); at z=-0.2m(c); at z=-0.3m(d) 87 0 0.5 1 1.5 2 −0.05 0 0.05 ζ (m) (a) 0 0.5 1 1.5 2 −0.1 −0.05 0 0.05 0.1 u W (m/s) (b) z=−0.15m 0 0.5 1 1.5 2 −0.1 −0.05 0 0.05 0.1 u W (m/s) (c) z=−0.25m 0 0.5 1 1.5 2 −0.1 −0.05 0 0.05 0.1 t(sec) u W (m/s) (d) z=−0.35m Figure 3.4: Comparison of wave records under opposing, uniform current (Case 1A) between experiment(o) and numerical solution(-): surface elevation(a) and oscillating velocities at z=-0.15m(b); at z=-0.25m(c); at z=-0.35m(d) 88 0 0.5 1 1.5 2 2.5 −0.1 −0.05 0 0.05 0.1 ζ (m) (a) 0 0.5 1 1.5 2 2.5 −0.5 0 0.5 u W (m/s) (b) z=−0.1m 0 0.5 1 1.5 2 2.5 −0.5 0 0.5 u W (m/s) (c) z=−0.2m 0 0.5 1 1.5 2 2.5 −0.5 0 0.5 t(sec) u W (m/s) (d) z=−0.3m Figure 3.5: Comparison of wave records under opposing, linear shear current (Case 2A) between experiment(o) and numerical solution(-): surface elevation(a) and oscillating velocities at z=-0.1m(b); at z=-0.2m(c); at z=-0.3m(d) 89 approach is less descriptive for the general effect of current vorticity, even though the Doppler effects are relatively well captured through it. 0 0.05 0.1 0.15 0.2 0.25 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 horizontal velocity(m/s) z (m) (a) Present Model Stoke’s Theory Experiment 0 0.1 0.2 0.3 0.4 0.5 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 horizontal velocity(m/s) z (m) (b) Present Model Stoke’s Theory Experiment Figure 3.6: Comparison of maximum wave velocity under following(a) and opposing(b) linear shear current; experiment(o), numerical solution(-) 3.4.2 Waves and Turbulent Currents with Bed Roughness For the turbulent current cases, reasonable estimation for the parameter b x is made to consider the wave-current interaction. Through experimental data of waves and turbu- lent currents, we can possibly assess an estimate for b x . Kemp and Simons(1982, 1983) are providing various types of physical measurements including turbulent characteristics. 90 Table 3.2: Wave and current characteristics in Kemp & Simons(1982, 1983) Case Wave Wave Current Current No. Height, H(m) Period, T(sec) Velocity, u Cα (m/s) Shear velocity, u ∗ C (m/s) WCA1 0.029 1.006 0.1983 0.00866 WCA3 0.0378 1.006 0.1983 0.00866 WCA4 0.0464 1.006 0.1983 0.00866 WCA5 0.0544 1.006 0.1983 0.00866 WDR1 0.0278 1.003 -0.1012 -0.00825 WDR3 0.0368 1.003 -0.1012 -0.00825 WDR4 0.045 1.003 -0.1012 -0.00825 WDR5 0.0487 1.003 -0.1012 -0.00825 They conducted experiments of the interaction between the wave and the turbulent cur- rent, either following or opposing direction in a 14.5m long, 0.457m wide and 0.69m deep channel. TestconditionslistedinTable3.2. Theygenerateddifferentbottomfrictionsfor the turbulent currents in the experiments by using either gloss-painted bed or triangular wooden strips placed at bed. Numerical simulations are set up using the actual scales of the experiments while using Δx = 0.01m and varying Δt for the stability(CFL ≈ 0.4). Such a small grid size(= 1/20h) will resolve the subdepth-scale motions of turbulence in depth-integrated model(Hinterberger et al.(2007)). Theprimaryissueinthesesimulationsistheapplicationofareasonablevalueforb x ,in order to ensure the turbulent effect of the interaction between waves and currents. Based on Yang et al.(2006), an appropriateb x value for an initial attempt in the simulation can be found. Then the final value of b x , which gives reasonable fit to the measurement in each case can be found as seen in Figure 3.7 and 3.8. As anticipated, general features of b x values are clearly seen from these results. The b x becomes negative (or positive) in following current (or opposing current) with its magnitude proportional to the wave 91 energy, i.e., radiation stress. It should be mentioned that the final values of b x are little different from those in Yang(2006), since they employed ‘idealized’ log-law profile for currents for the best fitting, instead of using real values from Kemp and Simons(1982, 1983). In the present study, however, ‘true’ shear stresses are directly applied to the simulation. Within this approach, the linear relationship between parameter b x and the relative radiation stress(S xx /τ x Cb h) in Equation 3.50 can be reasonably estimated. Figure 3.9 depicts the variation of the parameter b x with respect to the relative radiation stress. Indeed, there may exist two different ways to calculate S xx in wave-current co-existing condition. The first way is by using Equation 3.48 and the other by completely ignor- ing the presence of the underlying currents, i.e., letting u C = 0 in Equation 3.48. As seen in the figure, however, clear correlation between the parameter b x and the relative radiation stress can be found. This result is consistent with our assumption in Equa- tion 3.50. Moreover for both cases, almost identical regressed lines aboutb x values, each having either C b =-0.007674 or -0.00661, are obtained with minor errors. Meanwhile You(1996)’s empirical formulas(Equation 3.46 and 3.47) make inappropriate estimation as addressed before. Consequently, b x values can be dynamically determined from the following equation, using local characteristics of waves and currents. b x =−0.00661 S xx τ x Cb h (3.51) 92 0 0.05 0.1 0.15 0.2 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (a) b=−0.2 Model Experiment Current−only 0 0.05 0.1 0.15 0.2 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (b) b=−0.5 0 0.05 0.1 0.15 0.2 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (c) b=−0.8 0 0.05 0.1 0.15 0.2 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (d) b=−1.2 Figure 3.7: Mean-velocity profiles of combined waves and currents; (a)WCA1; (b)WCA3; (c)WCA4; (d)WCA5 93 −0.15 −0.1 −0.05 0 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (a) b=0.7 Model Experiment Current−only −0.25 −0.2 −0.15 −0.1 −0.05 0 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (b) b=0.9 −0.25 −0.2 −0.15 −0.1 −0.05 0 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (c) b=1.5 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 −0.2 −0.15 −0.1 −0.05 0 Time−averaged horizontal velocity(m/s) depth(m) (d) b=1.65 Figure3.8: Mean-velocityprofilesofcombinedwavesandcurrents; (a)WDR1; (b)WDR3; (c)WDR4; (d)WDR5 94 This equation opens, in the numerical implementation, a self-adaptive way to assess the parameter b x . Finally, the turbulent wave-current interaction which fully depend on the local wave and current properties can be considered in the model. 3.4.3 Bichromatic Waves and Uniform Currents in a Vayring Depth A complex hydraulic situation is expected when bichromatic waves interact with the current as explored in Dong et al.(2009)’ experiments. It is well known that bound long waves can be generated by the nonlinear wave-wave interactions in still water, but the influenceofunderlyingcurrentsontheirgenerationisyetnottotallyunderstood. Donget al.(2009)conductedlaboratoryexperimentsofbichromaticwavegroupspropagatingover the current in wave flume of 50m length, 3m width(reduced to 0.8m in the experiment) and 1m depth. A smooth bed of varying slope is also considered at the flume bottom as depicted in Figure 3.10. A bichromatic wave is composed of two sinusoidal waves of same amplitude with different frequencies as shown below. ζ =a 1 cos(k 1 x−2πf 1 t)+a 2 cos(k 2 x−2πf 2 t) (3.52) where a i ,f i and k i (i = 1,2) are the wave amplitudes, frequencies and wave numbers, respectively. Underlying currents will have either following or opposing directions to the wave propagation. It is also revealed from their measurements that the current variation along the changing depth can be effectively estimated from continuity and is almost depth-uniform. Thus, the current field is idealized through continuity in the simulation. The largest Froude number of currents in the experiments is 0.054 which ensures a weak 95 −300 −200 −100 0 100 200 300 −3 −2 −1 0 1 2 3 S xx /τ Cb h parameter b Following Currents Opposing Currents (1):S xx w/o u C (2):S xx w/ u C b by empirical formular Fitted line for (1) Fitted line for (2) Figure3.9: Distributionofparameterbaccordingtoradiationstress(S xx scaledbycurrent bed stress(τ Cb )); Fitted line is represented by solid and dashed line 96 current 1:60 1:20 1:40 1:10 x=8m 0.25m h=0.5m wave x=9m x=21m x=27m x=0m x=33m x=34m Figure3.10: ExperimentalsetupofbichromaticwavesunderanambientcurrentsbyDong et al.(2009) Table 3.3: Bichromatic waves and current characteristics in Dong et al.(2009) f 1 (Hz) f 2 (Hz) a 1 (cm) a 2 (cm) Current(cm/s) 1.1 0.9 1.75 1.75 6, -6, 12, -12 current condition without hydraulic jump. Table 3.3 lists the details of the bichromatic waves and the current conditions used in the experiments. Numerical simulation is set up based on the above information while using Δx = 0.01m and Δt corresponding to CFL (=0.4) conditions. Runtime is set to 200s and this ensuresatleast100bichromaticwavesand20sub-harmonic longwaves(f 1 −f 2 =0.2Hz) are generated. To analyze the sub-harmonic bound long waves and obtain local wave amplitude, FFT(Fast Fourier Transform) technique is adopted using the same frequency rate(40Hz) as Dong et al.(2009). Figure 3.11 to 3.13 show simulation results compared with experimental data, presenting the spatial variation of H rms along the flume. In Figure 3.11, which is the result of pure wave condition, calculatedH rms of both primary and bound long waves are in good agreement with experiments, as well as with the ana- lytical solution derived by Longuet-Higgins and Stewart(1962). Also, the overall shoaling process due to the varying depth is well represented. 97 5 10 15 20 25 0 1 2 3 4 5 6 7 x(m) Hrms(cm) (a) Numerical Measured 5 10 15 20 25 0 0.5 1 1.5 2 2.5 3 Long wave amplitude(mm) x(m) (b) Numerical Measured Analytic(Longuet−Higgins & Stewart(1962) Figure3.11: Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbounding long waves;(without currents) 98 For the following current conditions, both measured and calculated primary wave H rms are clearly decreased by about 10% due to the Doppler shift effect as shown in Figure 3.12 (a). On the other hand, increased H rms of primary waves are calculated for theopposingcurrentcase,andthisisalsoconsistentwiththemeasurementinFigure3.13. As seen from Figure 3.12 (b) and 3.13 (b), in both following and opposing current cases, Doppler effect is well characterized for sub-harmonic bound long waves. Finally, to see the current effect on wave energy modulation clearly, calculated ampli- tude spectra of bichromatic waves are presented in Figure 3.14. It is revealed that some of wave energy at higher frequency component is relocated to the lower frequency for all cases. In the opposing current case, such tendency becomes most dramatic and this is consistent with Dong et al.(2009)’s observations. Comparisonsbetweenwaveamplitudespectraatthreedifferentlocation(x=7m,18m, and24m)indicatealsothatasthewavespropagateovershoalingdepth,thesub-harmonic component energies develop for pure wave and for following current cases. Meanwhile, it is not obvious for the opposing current case because relatively strong current exists at shallower depth. 3.5 Summary A set of depth-integrated equations describing combined wave-current flows is derived in this chapter. It is well known that the presence of waves over the mean flow introduces additional stresses. With a parameterb, a linear shear profile can be defined to represent such stresses by wave-current interaction. 99 5 10 15 20 25 0 1 2 3 4 5 6 7 x(m) Hrms(cm) (a) Numerical Measured 5 10 15 20 25 0 0.5 1 1.5 2 2.5 3 Long wave amplitude(mm) x(m) (b) Numerical Measured Figure3.12: Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbounding long waves;(following currents) 100 5 10 15 20 25 0 1 2 3 4 5 6 7 x(m) H rms (cm) (a) Numerical Measured 5 10 15 20 25 0 0.5 1 1.5 2 2.5 3 Long wave amplitude(mm) x(m) (b) Numerical Measured Figure3.13: Spatialvariationof;(a)Hrmsofbicromaticwaves;(b)Amplitudeofbounding long waves;(opposing currents) 101 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 Frequency (Hz) Amplitude(cm) u C =0, x=7m u C =12cm/s, x=7m u C =−12cm/s, x=7m u C =0, x=18m u C =12cm/s, x=18m u C =−12cm/s, x=18m u C =0, x=24m u C =12cm/s, x=24m u C =−12cm/s, x=24m Figure 3.14: Calculated amplitude spectra of bichromatic waves with and without cur- rents 102 Doppler shift effect embedded in a system of equations is elucidated by limiting the equations into one dimensional linear case. It is confirmed that Doppler shift property retains in the equation set. Toimplementthepresentmodelnumerically,afourth-orderMUSCL-TVDschemein- corporatedwithapproximatedRiemannsolver,isadoptedforleadingordertermswhereas a cell-averaged finite volume method is utilized for higher-order terms. Accuracy of the model is then examined through three experimental data sets. The simulated results on long waves propagating over either uniform or sheared current, agree well with the measurements. It is found that the vorticity of the current needs to be taken into ac- count when describing the velocity of wave-current flows more accurately. Through the inclusion of rotational terms in the velocity profiles, we’re able to practically consider the effects of current’s vorticity on the flow field. Validity of the inclusion of additional stress term is also investigated through the simulation on turbulent current interacting with waves. Using Kemp and Simons(1982, 1983)’s experimental data, reasonable estimation of parameter b is provided. Based on provided values, b has been physically defined through linear relationship with relative radiation stress of waves. Asafinalattempt,themodelisappliedtoacomplexconfiguration,wherebichromatic waves interact with spatially varying currents. Subharmonic bound long waves generated by wave-wave interaction are well recreated in the presence of background currents and the comparisons with measurements show good agreements. Consequently, the results clearly indicate that the present model is capable of predicting both wave and current flows accurately. 103 Chapter 4 A Depth-Integrated Model for Free Surface Waves Propagating over Fluids with Weak Vertical and Horizontal Density Variation Inthischapter, weconsiderthedensitychangeoffluidsinthedepth-integratedlongwave model. By allowing horizontal and vertical variation of fluid density, a depth-integrated modelforlonggravitywavesovervariabledensityfluidhasbeendeveloped,wheredensity change effects are included as correction terms. In particular, two-layer fluid system is chosen as vertical density variations, where interfacial wave effects on the free surface is accounted for through direct inclusion of velocity component of the interfacial wave. For the numerical implementation of the model, a finite volume scheme coupled with approximate Riemann solver is adopted for leading order terms while cell-centered finite volume methods are utilized for others. Numerical tests have been performed to verify the model, in which the density field is configured to vary either horizontally or vertically. For horizontal variation of fluid density, a pneumatic breakwater system is simulated and fair agreement is observed 104 between computed and measured data, showing that the current induced by the upward bubblefluxisresponsibleforwaveattenuationtosomedegree. Toinvestigatetheeffectsof internalmotionsonthefreesurface,atwo-layerfluidsystemwithmonochromaticinternal wave motion is tested numerically. Comparison of simulated results with measured and analytical data shows good agreement. Lastly, nonlinear interaction between external- and internal-mode surface waves are studied numerical and analytically, and the model is shown to have nonlinear accuracy similar to published Boussinesq-type models. 4.1 Introduction Freshwater and seawater are often coexisting in the estuaries in the coastal regions. Thus the fluid densities are subject to change horizontally and vertically due to thermal and saline changes. This density variation will result in complex physical processes on the contrary to the uniform density fluids. For example, in the horizontally varying density fluids, the hydrostatic pressure will be imbalanced between the fluids and this pressure difference will cause the fluid to move. Armi and Farmer (1986) provide an example of the exchange of the two fluids through the contraction, and Wood (1970) provides an example of the classical lock-exchange problems. Moreover, in the density varying flow, the momentum needed to convey the lighter fluid is less, and that implies a significant change in the kinematics of the problem. On the other hand, the density variation inherently implies a three dimensional prob- lem in the natural oceanic system Mellor (1991). Therefore, the 3D hydrodynamic model 105 coupled with the density transport model is commonly adopted for the simulation of the variable density fluid flows (e.g., Kanarsk and Maderich 2003). Different approaches have been made to solve this problem. For instance, in the horizontallyvaryingdensityflows,theinitialsetupofthehorizontallyseparatedtwofluids with different densities is assumed to study the gravity current generated at the lock- exchange configurations (e.g., Lowe et al. (2005)), or in the bi-directional flows through contractions(e.g.,ArmiandFarmer(1986)). Thosestudiesfocusedontheinternalphysics of the problem such as the diffusion and the mixing caused by the barotropic forces. Recently, Leighton et al.(2010) examined the effect of the horizontally varying density conditions on the free surface physics. They presented one-phase, one-dimensional model for vertically well-mixed shallow flows with horizontal density gradient; however, they neglected the vertical variation of the fluid density. In the vertically stratified fluids, they are usually simplified to a two-layer system with different densities as in the estuarine system. In the stable stratification, there will be no driving force and therefore the two fluids will stay without any motion. However, if an initial disturbance is introduced at the interface, this disturbance will propagate as an oscillating form in order to keep the balance between the gravity and the buoyancy. Meanwhile, this internal motion will affect the hydrodynamics of the entire water column including the free surface. Topographical configuration could act as a source for internal wave generation at the interfaceofthetwo-layersystem(e.g.,FarmerandArmi(1999),HelfrichandMelville(2006)). Researchers are interested in the physical properties of these internal waves such as am- plitude and frequency dispersion without compromising the free surface waves. They 106 imposed the rigid-lid assumption on the upper boundary condition, thus they simplified the relevant physical problem in order to focus solely on the internal wave physics (e.g., Evans and Ford(1996)). Numerous models have been developed for investigating the internal wave propagation over the shallow and the deep waters applying the rigid-lid as- sumption(KoopandButler(1981),SegurandHammack(1982),ChoiandCamassa(1999), Debsarma et al.(2010)). However, there are some evidences showing that sub-surface movement can affect the freesurface. Forinstance,thetheoreticalstudiesofLamb(1932),Phillips(1977),Hwunget al.(2009)andtheexperimentbyUmeyama(2002)whichshowsthesurfacewaveexcitation by internal waves. Moreover, the observation by Elachi and Apel(1976) have confirmed the importance of including the free surface. Although some numerical models apply the free surface condition for long internal waves (Choi and Camassa(1996), P˘ ar˘ au and Dias(2001), Lynett and Liu(2002b), Nguyen and Dias(2008), Liu and Wang(2012)), they did not consider the free surface dynamics appropriately. In reality there are many occasions where the free surface gravity waves exist with in- ternal waves and interact with each other. Therefore, nonlinear interactions between internal and surface waves have also been explored in the literatures experimentally (Lewis et al.(1974), Joyce(1974)) and theoretically(Ma(1982), Donato(1999), Liu(2006), Selezov et al.(2010), Craig et al.(2011)). Triad resonance interaction is investigated in (Alam(2012), Hill and Foda(1998)). Most of the previously mentioned numerical models solvethetwosetsofequations, whicharegeneratedfromthetwo-layerfluidsystem. This will require more computational time and resources than the one equation set (Lynett and Liu(2004)). The depth-integrated approach, for the shallow water regime has been 107 well established for the past few decades. The vertical structure of some physical quan- tities such as velocity and pressure fields is represented by polynomials and successfully describes the long wave formulations. In the present study, for the purpose of efficiently but physically accurate modeling of surface waves over varying density, a new attempt based on depth-integration has been made. This attempt will include the effects on the surface waves by density variation as well as internal motions at the interface of two immiscible fluids with different den- sities. Density-variation effects on the free surface are included as correction terms in the Boussinesq-type models. The internal wave effect on the free surface is considered through direct inclusion of internal wave kinematics as weak components of the velocity profile. Long wave assumptions for both surface and internal waves are also implemented in the model derivation. Thus, the present model is designed to study the problem of long wave propagation over the density-varying fluids. This chapter is organized as follows. In the following section, depth-integrated model equations for long surface waves over variable density fluid are derived mathematically. Then, the model is applied to the surface wave propagation over either horizontally or vertically varying density fluids. An exercise to investigate nonlinear property of the designed model has been also performed as surface and internal wave interactions. The resultsarethencomparedwithexperimentalandanalyticaldataforthevalidation. Final conclusion is made in subsection 4.5. 108 4.2 Mathematical Formulation A mathematical model for weakly dispersive wave propagation over a varying-density fluid is derived in this section. A standard Boussinesq-type approach (e.g., Lynett and Liu(2002a), Kim et al.(2009)) is followed during the derivation to obtain a vertically- independent equation set. Worthy of note here is that we use the term ‘Boussinesq’ as related to weakly dispersive and nonlinear; we are not implying the same namesake commonly referenced in the oceanography field when including buoyancy effects. 4.2.1 Governing Physics and Boundary Conditions Toderiveamathematicalgoverningsystemofequations,thebasicphysicalparametersare definedasinFigure4.1. Atwo-layerfluidstablyseparatedbyasmall,butfinite,thickness pycnoclineisassumedtobeincompressibleandimmiscible. Appropriateparametricscales are now introduced to normalize the governing equations as well as boundary conditions. Typical for long wave scaling, ℓ 0 and h 0 , the wavelength and water depth, are used for horizontal and vertical coordinates, respectively. Additionally, ρ b is a characteristic fluid density of the entire fluid system. Using this set of scaling parameters, normalized variables are introduced as, (x,y)= (x ′ ,y ′ ) ℓ 0 , z = z ′ h 0 , ζ = ζ ′ h 0 , h= h ′ h 0 , t= t ′ √ gh 0 ℓ 0 , (u,v)= (u ′ ,v ′ ) √ gh 0 , w = w ′ µ √ gh 0 , p= p ′ ρ b gh 0 , ρ= ρ ′ ρ b , µ = h 0 ℓ 0 , ν h t = ν h ′ t αh 0 √ gh 0 , ν v t = ν v ′ t βh 0 √ gh 0 (4.1) 109 where a prime denotes dimensional variable, (u,v) and w represents instantaneous hor- izontal and vertical velocity, respectively. ζ is free-surface elevation which is a function of (x,y,t), p is pressure, g is gravitational acceleration and ρ is density. Dimensionless parameters α and β are used to scale the horizontal and vertical eddy viscosities ν h t and ν v t (Kim et al.(2009)), and velocities are scaled following long wave theory. Note that there is no wave amplitude scale, which would otherwise appear in the ζ, pressure, and velocity scalings; it is implicitly assumed that wave amplitude is the same order as the water depth. h (x , y ) x y z ] (x , y , t ) z = 0 z = -h z = z 0 ȡ l (x , y ) ȡ u (x , y ) Figure 4.1: Sketch of long wave propagation over two-layer fluids with horizontal density variation To examine the effects of spatially varaiable density on free surface waves, it is neces- sarytoparameterizethemagnitudeofthedensitychange. Theparameterγ isintroduced for quantifying the density variation, and is given as γ = Δρ ′ ρ b (4.2) 110 The term Δρ ′ is simply the characteristic density variation in the physical system to be examined. Following Kim and Lynett (2011), the dimensionless form of the spatially-filtered continuity and Navier-Stokes equations for incompressible flow are ∇·u+w z =0 (4.3) ρ(u t +u·∇u+wu z )+∇p=αµ ∇· ρν h t ∇u + β µ (ρν v t u z ) z (4.4) µ 2 ρ w t +µ 2 u·∇w+µ 2 ww z +p z +ρ=αµ 3 ∇· ρν h t ∇w +βµ (ρν v t w z ) z (4.5) where∇ is the horizontal derivative operator, subscriptsz andt function as vertical and time differentiations, respectively, andu=(u,v) is the horizontal velocity vector. Conditions applied at free surface and at bottom boundaries in Fig 4.1 are expressed in dimensionless form, as well: w =ζ t +u·∇ζ at z =ζ (4.6) w+u·∇h=0 at z =−h (4.7) 111 Equations 4.3 to 4.7 represent the primitive equation set needed to describe the fluid motion shown in Figure 4.1. 4.2.2 Derivation of Boussinesq-type Equations for Dispersive Waves over Variable Density Fluid With a few dimensionless parameters introduced, µ 2 is chosen to expand the physical variables. The assumption is made that the effect of density stratification and viscosity are as weak, and of a similar order as dispersion. That is, O(µ 2 )=O(µβ,γ )≪1 (4.8) The physical parameters, p,u,v,w are expanded as power series of µ 2 , f = ∞ X n=1 µ 2n f n . (4.9) The leading order terms of Equation 4.5 yields the hydrostatic condition, (p 0 ) z +ρ 0 =0 (4.10) which accordingly guarantees, u(x,y,z,t)=u 0 (x,y,t)+O(µ 2 ) (4.11) 112 Integrating the continuity equation over depth and applying the free surface and bottom boundary conditions, Equations 4.6 and 4.7, will give the relationship between horizontal and vertical velocities as, w 0 =−z(∇·u 0 )−∇·(hu 0 )=−zS−T (4.12) In order to determine the horizontal velocity profile, horizontal vorticity(ω) is exam- ined ω ′ =u ′ z −∇w ′ =µ 2 r g h 0 ω 1 +O(µ 4 ) (4.13) which demonstrates that horizontal vorticity can at most first appear at O(µ 2 ) within our scaling. Assuming thatω 1 is not zero permits rotational effects induced by bottom stress to be directly included in the velocity profile (Kim et al.(2009)). From the horizontal vorticity expression, the vertical profile of horizontal velocity can be approximated by vertically integratingω 1 from−h to z, u 1 = u 1 | z=−h − 1 2 z 2 ∇S+z∇T − 1 2 h 2 ∇S−h∇T + Z z −h ω 1 dz+O(µ 2 ) (4.14) The integral of vorticity term appearing above remains as a “residual” velocity compo- nent, and can be specified depending on the particular physics of the configuration at hand (e.g. bottom or free surface stress, or stratification effects). 113 The horizontal velocity up to O(µ 2 ) can be then expressed as, u=u 0 + µ 2 u 1 | z=−h − 1 2 z 2 ∇S−z∇T + 1 2 h 2 ∇S−h∇T (4.15) + Z z −h ω 1 dz +O(µ 4 ) and following Nwogo’s (1993) approach , we can defineu α evaluated at z =z α as, u α =u 0 + µ 2 u 1 | z=−h − 1 2 z 2 α ∇S−z α ∇T + 1 2 h 2 ∇S−h∇T (4.16) + Z zα −h ω 1 dz +O(µ 4 ) Subtracting Equation 4.16 from Equation 4.15 finalizes the expression of u in terms of u α as, u=u α + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T + Z z zα ω 1 dz +O(µ 4 ) (4.17) As a main purpose of this study is to extend the derivation to include internal motion due to stratification, the horizontal velocity of the background internal motion, u i (z), should be included. Here, this component is interpreted as part of the residual vorticity acting on the barotropic wave, i.e. Z z zα ω 1 dz =u i (z)+ Z z zα ω s 1 dz (4.18) 114 whereω s 1 represents the horizontal vorticity due to a bottom stress only. This decompo- sition yields u = u α +µ 2 u i (z) (4.19) + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T + Z z zα ω s 1 dz + O(µ 4 ) Note that in the above, it is implicitly admitted that O(γ) =O(µ 2 ), asu i (z) is usually scaled by reduced gravity (e.g., Lynett and Liu(2002b)). As a result, the direct inclusion of density-driven internal kinematics into the velocity structure allows for consideration internal motion effects on free surface waves. While the above expression is flexible in terms of theu i (z) that can be accommodated, the theory is only valid if the magnitude of this component is small. For a quantitative description of u i (z), specific cases will be chosen later for ver- ification. Regardless, the horizontal velocity with rotational terms( R z zα ω s 1 dz ′ ) remains undetermined; these rotational terms are assumed related only to bottom stress(τ) (i.e. Kim et al.(2009)). Imposing a linear stress profile from zero at the free surface toτ b at the bed, ω s 1 = τ b ρν v t ζ−z ζ +h , (4.20) 115 isalltheinformationrequired. Integrationofω s 1 producesarotationalvelocitycomponent of O(µ 2 ) as follows, Z z zα ω s 1 dz ′ =Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) (4.21) whereΨ=τ b /{ρν v t (ζ +h)}. With this, horizontal velocity up to O(µ 2 ) can be finally expressed as, u=u α + µ 2 u i +µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T (4.22) + µ 2 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) + O(µ 4 ) Utilizing the horizontal velocity profile above, the exact continuity equation can be given in a depth-integrated format. Integrating Equation 4.3 from−h toζ and applying boundary conditions(Equation 4.6 and 4.7) gives, ∇· Z ζ −h udz+ζ t =0 (4.23) To express this equation in terms ofu α , Equation 4.22 is substituted, and after manipu- lation yields ζ t +∇·{(ζ +h)u α }+µ 2 (N D +N B +N I )=O(µ 4 ) (4.24) 116 where the second order terms are, N D =−∇· (ζ +h) ζ 2 −ζh+h 2 6 − z 2 α 2 ∇S+ ζ−h 2 −z α ∇T (4.25) N B = ∇· " ψ(ζ +h) ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# (4.26) N I =∇· n (ζ +h)u i o (4.27) whereu i is the depth-averaged horizontal velocity vector due to internal motion. While integration of the continuity equation was straightforward, derivation of a depth-integrated momentum equation, however, needs a fairly complex procedure due to nonlinearity in the equation and the desired removal (through substitution) of the hy- drodynamic pressure term. To facilitate the final expression of our equations in a closed (non-integral) form, the first functional restriction will be placed on the internal den- sity field. Here a simple vertical density structure is assumed, however realistic enough to resemble a physical configuration found commonly in oceans and lakes (e.g. Kao et al.(1985)). A two-layer density profile is thereby specified as, ρ(x,y,z)=ρ 0 (x,y) 1−γtanh z−z 0 δ (4.28) 117 where ρ 0 is some average density of the two-layer system, z 0 is the midpoint of the tanh inflection separating the upper and lower layers, and δ is the finite transition pycnocline thickness. Note that all terms in the above density profile expression are dimensionless; density has been scaled by ρ b and the three distance terms inside the tanh operator have been scaled by the depth,h. A discrete two-layer density configuration thus can be accommodated using the tanh function, eliminating any density discontinuity issues. The steps to derive a final, depth-integrated equation are somewhat tedious, and are includedindetailinAppendix F. Theresultantformofthedepth-integratedmomentum equation appears as, (u α ) t + u α ·∇u α +∇ζ + ∇ρ 0 ρ 0 (ζ +h) 2 +γR v P + µ 2 R D +R B +R I +R h P +ξ − αµ 1 ρ 0 ∇· ρ 0 ν h t ∇u α +βµν v t ∇S + βµ 1 ρ 0 τ b ζ +h −βµν v t u i z | z=ζ − u i z | z=−h = O(µ 4 ,αµ 3 ,βµ 3 ) (4.29) in which the higher order terms are defined as, R D = 1 2 z 2 α ∇S t +z α ∇T t − 1 2 ∇ ζ 2 S t −∇(ζT t )+T∇T + 1 2 ∇ z 2 α u α ·∇S +∇(z α u α ·∇T)+ 1 2 ∇ ζ 2 S 2 − 1 2 ∇ ζ 2 u α ·∇S −∇(ζu α ·∇T)+∇(ζTS) (4.30) 118 R B = (ζ−h) 2 (Ψζ) t − ζ 2 −ζh+h 2 6 Ψ t + Ψ z 2 α 2 −ζz α t (4.31) + (ζ−h) 2 ∇{u α ·(Ψζ)}− ζ 2 −ζh+h 2 6 ∇(u α ·Ψ) + ∇ u α · Ψ 1 2 z 2 α −ζz α −Ψ ( ζ 2 −ζh−2h 2 S 6 + (ζ +h)T 2 ) R I = u i t +∇ u α ·u i (4.32) R h P = ∇ρ 0 ρ 0 −2ζ +h 6 S t +u α ·∇S−S 2 − 1 2 (T t +u α ·∇T −ST) (4.33) R v P = δ ζ +h {(u α ) t +u α ·∇u α }× ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ − 1 ρ 0 ∇ ρ 0 δlncosh ζ−z 0 δ + 1 ρ 0 1 ζ +h Z ζ −h ∇ ρ 0 δlncosh z−z 0 δ dz (4.34) 119 Andξ = ξ x ,ξ y are described as ξ x = −v α n (z α ) x (z α S y +T y )−(z α ) y (z α S x +T x ) o (4.35) − n (v α ) x −(u α ) y o "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) S y + z α − (ζ−h) 2 T y # − v α " ψ y 1 2 z 2 α −z α ζ x − ζ 2 −ζh+h 2 6 (ψ y ) x + (ζ−h) 2 (ψ y ζ) x − ψ x 1 2 z 2 α −z α ζ y + ζ 2 −ζh+h 2 6 (ψ x ) y − (ζ−h) 2 (ψ x ζ) y # − n (v α ) x −(u α ) y o ψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) + v i α n (v α ) x −(u α ) y o ξ y = u α n (z α ) x (z α S y +T y )−(z α ) y (z α S x +T x ) o (4.36) + n (v α ) x −(u α ) y o "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) S x + z α − (ζ−h) 2 T x # + u α " ψ y 1 2 z 2 α −z α ζ x − ζ 2 −ζh+h 2 6 (ψ y ) x + (ζ−h) 2 (ψ y ζ) x − ψ x 1 2 z 2 α −z α ζ y + ζ 2 −ζh+h 2 6 (ψ x ) y − (ζ−h) 2 (ψ x ζ) y # + n (v α ) x −(u α ) y o ψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) + u i α n (u α ) y −(v α ) x o whereu α =(u α ,v α ) andΨ=(ψ x ,ψ y ). In Equation 4.29, each higher order term implies weak effects by frequency dispersion(R D ), bottom stress(R B ), internal motion(R I ) and density variation(R h P andR v P ). 120 4.3 Limiting Cases of Derived Model 4.3.1 Linear, Inviscid Equations Eliminating all nonlinear terms in ζ andu α , assuming that ζ/h is very small, and con- verting to dimensional form yields the linear system: ζ t + ∇·(hu α )−∇· h 3 6 − hz 2 α 2 ∇S− h 2 2 +hz α ∇T + ∇· hu i = 0 (4.37) (u α ) t + g∇ζ + 1 2 z 2 α ∇S t +z α ∇T t + u i t +∇ u α ·u i + ∇ρ 0 ρ 0 1 2g + 1 6 S t − 1 2h T t +(R v P ) linear = 0 (4.38) where g is gravity and (R v P ) linear is constructed in the same manner as Equation 4.34, with the only difference being the nonlinear u α ·∇u α removed from the first line and all ζ set to zero. Both of the above equations are presented such that the standard, free-surface-only wave terms are shown in the first line, and the terms driven by the internal flow and stratification are given in the subsequent lines. In the second line of the linear momentum equation ( 4.38) above are the forcing terms associated with the 121 depth-averaged internal flow, while in the third line are the two terms due to horizontal and vertical density gradients, respectively. 4.3.2 First-Order Nonlinear, Inviscid Equations Here, we include only first order nonlinear terms;ζ/h is no longer negligible. The dimen- sional system becomes ζ t + ∇·(Hu α )−∇· h 3 6 − Hz 2 α 2 ∇S− h 2 2 +Hz α ∇T + ∇· Hu i = 0 (4.39) (u α ) t + u α ·∇u α +g∇ζ + 1 2 z 2 α ∇S t +z α ∇T t − ∇(ζT t )+T∇T + 1 2 ∇ z 2 α u α ·∇S +∇(z α u α ·∇T) + u i t +∇ u α ·u i +R v P + ∇ρ 0 ρ 0 H 2g + ζ 3h S t + 1 6 S t +u α ·∇S−S 2 − 1 2h (T t +u α ·∇T −ST) = 0 (4.40) whereH =h+ζ. Again, the standard Boussinesq-type terms are given first in the above equations. This nonlinear system will be analyzed later in this chapter to determine its analytical behavior. It is noted here that, in the absence of any vertical or horizontal density gradient, all “new” terms disappear, with the exception of the u i in both the 122 continuity and momentum equations. While these terms have been discussed in this study in the framework of internal motion, the effects of any type of imposed external flow could be accommodated throughu i . 4.4 Model Validation In this section, the model is applied to various types of problems in which the fluid density field is configured to vary either horizontally or vertically. Calculated results fromthesesimulationsarethencomparedtoexperimentaldataoranalyticalsolutionsfor verification. Completedetailsofthenumericalschemeusedinthefollowingexamplescan be found in G. It is noted that for the remainder of the chapter, as simplified analytical and laboratory cases are examined, all terms related to viscosity (with coefficients α or β) are neglected. 4.4.1 Horizontally Varying Fluid Density - Pneumatic Breakwater In this section, free surface wave propagation in a fluid with a horizontal density gradient is considered. The fluid is assumed to be well-mixed vertically such that stratification effects(R I and R v P ) can be ignored. While this physical setup is a reasonable represen- tation of a common coastal situation (e.g. wind waves near a river mouth or inside an estuary), published experimental data in this configuration is very limited. The reason for this lack of available data is likely due to the complexity of the experimental setup; a horizontal density gradient is not stable, and continuous forcing is required to maintain it. With this in mind, the model derived here will be compared with the available data such that the behavior of the model might be assessed. Data from a two-phase flow 123 (air/water), pneumatic breakwater experiment is chosen. The bubbly area around the pneumatic breakwater is regarded as low-density fluid while the surrounding clear water maintains a constant density. Most the studies of pneumatic breakwater focus on the surface currents generated by bottom-generated bubbly flows, which are believed to be the main drivers of wave attenuations(Taylor(1955), Lo(1991)). Nonetheless, a sudden drop of density at the bub- bly zone is responsible for some wave energy blockage, in that less momentum can be transported through the lower density, bubbly fluid. In this section the present model will be applied to pneumatic breakwater simulation, in which horizontal density gradient effects in wave attenuation are evaluated. Thereareafewlaboratorystudiesonpneumaticbreakwaters(e.g.,Straubetal.(1959), Bulson(1963), Zhang et al.(2010)). Among these, Zhang et al.(2010) carried out a labo- ratory test on the performance of pneumatic breakwaters at a 1:15 model scale. Using relatively long incident waves, they measured wave height behind the pneumatic break- water to quantify wave transmission through the bubble system. Figure 4.2 describes experimental setup of pneumatic breakwater. The wave tank has dimensions of 2m in width, 1.8m in depth and 69m in length, and the air discharge pipe is installed along the bottom. Through orifices of 0.8mm in diameter, spaced 0.01m along the pipe, air is discharged into the water to create a bubble curtain. To examine the breakwater performance as a function of different hydraulic conditions, two different wave periods and four different air fluxes are tested. Table 4.1 summarizes experimental conditions at laboratory scale. 124 h = 0.8m ȡ w H i = 0.236m H t air discharger air bubble ( ȡ a ) Figure 4.2: Experimental setup of pneumatic breakwater(Zhang et al.(2010)) Table 4.1: Experiment conditions of pneumatic breakwater(Zhang et al.(2010)) Water depth 0.8 m Wave height 0.236 m Wave period 1.29s, 1.55s Air amount 10 m 3 /h, 20 m 3 /h, 30 m 3 /h, 40 m 3 /h Forthenumericalsimulation,acomputationaldomainbasedonthewavetankdimen- sions is constructed using Δx = 0.025m. A flexible time step, Δt is adopted according to CFL stability condition. During the simulation, transmitted wave height is measured behind the low density region. One remaining issue in creating numerical configuration is to conceptualize the bubbly field into an “equivalent” fluid zone of low density. To this end, some simplifying assumptions are made. For shallow flows, the lateral distribution of bubble concentration at any depth is given by a top-hat profile. Moreover, as seen in Figure 4.3 (a), bubbles are assumed to be scattered, forming a linear slope (=1/10) under the wavy condition, as indicated in Zhang et al.(2010). Slip velocity(u s ) of the 125 ȡ w Air-bubble i t ( ) ( ) ȡ w mixture ȡ (z) S = 1/10 (a) W ȡ w Mean density ȡ (b) Figure 4.3: Simplified density field of bubbly fluids; (a)Depth-varying ρ(z) (b)Depth- constant ρ rising bubble is given empirically by(Clift et al.(1978)), u s =                4474m/s×r 1.357 b if 0≤r b ≤0.7mm 0.23m/s if 0.7mm<r b ≤5.1mm 4.202m/s×r 0.547 b if r b >5.1mm (4.41) wherer b is bubble radius and approximately the same as orifice radius(Ma et al.(2012)). Fluid velocity(u w ) is determined by (Milgram et al.(1983)), u w =1.5u s (4.42) 126 Table 4.2: Parameters for simplified density fields of bubbly area Air amount(q l ) Slip velocity(u s ) Fluid velocity(u w ) ρ Width(W) 5 m 3 /h/m 0.11 m/s 0.67 m/s 938.52 0.08m 10 m 3 /h/m 0.11 m/s 0.67 m/s 880.77 0.08m 15 m 3 /h/m 0.11 m/s 0.67 m/s 826.20 0.08m 20 m 3 /h/m 0.11 m/s 0.67 m/s 774.89 0.08m The total rising velocity of bubbles is then (u s +u w ). Guided by these assumptions, the density of the bubbly zone can be calculated as ρ(z)= 1− q a 2S(h+z)(u s +u w ) ρ w + q a 2S(h+z)(u s +u w ) ρ a (4.43) in which ρ a and ρ w represent densities of air(= 1.269kg/m 3 ) and water(= 1000kg/m 3 ), respectively. S denotes side slope of bubbly zone and q a represents injected air volume per unit width. For depth-integrated model use, ρ(z) is depth-averaged, as depicted in Figure 4.3 (b). Calculated quantities for density fields are listed in Table 4.2. The computed results of the transmission coefficient(K t = transmitted wave height divided by incident wave height) are plotted in Figure 4.4 and Figure 4.5 with measure- ment. As expected, the effects of the density transition can be found in both results, indicating greater air volume discharge yields less transmission through the pneumatic breakwater. Comparison between calculation and measurement shows good agreement although the numerical simulations have a clear over-prediction bias. The likely reason for this discrepancy is the numerical neglect of the bubble-induced currents. However, the model provides the proper parametric trends and dependencies, and the accuracy is reasonable. 127 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Air amount(m 3 /h/m) Transmission Coeff. (K t ) Measured Calculated Figure 4.4: Comparison of transmission coefficient(K t ) between measurement and calcu- lation; T=1.55 sec 4.4.2 Waves Excited by Internal Motion: Linear Surface Waves In this section two-layer internal flow is considered, in which an internal wave travels along the interface. Being concerned with vertical stratification only, terms associated with horizontal density variation (R h P ) are neglected. Itiswell-knownthatthepresenceofinterfacialmotioninatwo-layerfluidexcitesfree surface motion (Lamb(1932), Phillips(1977)). Lamb(1932) obtained an exact solution about the surface disturbance due to an internal wave using linear potential theory. Such 128 0 5 10 15 20 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Air amount(m 3 /h/m) Transmission Coeff. (K t ) Measured Calculated Figure 4.5: Comparison of transmission coefficient(K t ) between measurement and calcu- lation; T=1.29 sec disturbances at the free surface propagate along with internal waves, whose dispersive properties are locked to the internal motion. Withinthepresentedderivation, theprofileu i (z)hasuptothispointbeenconsidered arbitrary. To account for internal wave effects on the free surface modes, a reasonable description of internal wave velocity u i (z) is needed here. First, the simplified two-layer system as shown in Figure 4.6 is presented. Density and undisturbed thickness of upper 129 K s ȡ u z =- h u K i z =-(h u + h l ) ȡ l z h u Figure 4.6: Simplified two-layer system with interfacial wave a a Idea underlying the plot assumes internal-mode interfacial wave(η i ) generates internal-mode surface wave(η s ), For the definition of internal-mode, see P˘ ar˘ au and Dias(2001) layer are ρ u and h u , while the lower layer has ρ l and h l . The disturbances propagating along the surface(η s ) and interface(η i ) are defined as simply as η s =a s cosϑ i (4.44) η i =a i cosϑ i (4.45) in which a s and a i represent wave amplitudes at the surface and interface, respectively, and ϑ i =k i x−σ i t+ǫ is the phase function with arbitrary shift ǫ. Note that, k i and σ i are assumed to follow the internal wave dispersion relation(Lamb(1932)) as below, 130 ρ u ρ l tanh k i h u tanh k i h l +1 σ i 4 −gk i tanh k i h u +tanh k i h l σ i 2 +tanh k i h u tanh k i h l 1− ρ u ρ l g 2 k i 2 =0 (4.46) Following Liu(2006), who developed an analytic solution based on potential theory, the linear solution of u i can be expressed as, u i (z)=                Acosh k i (z+h l ) +Bsinh k i (z+h l ) cosϑ i (η i ≤z≤η s ) Ccosh (k i (z+h u +h l ) cosϑ i (−h l ≤z≤η i ) (4.47) where A,B,C are A= a i σ i " g 1− ρ l ρ u + ρ l σ i 2 ρ u k i tanh(k i h l ) # (4.48) B = a i σ i k i (4.49) C = a i σ i k i sinh(k i h l ) (4.50) 131 Thedepth-averagedinternalvelocity,u i ,isneededbythemodel,andcanbecalculated as u i = 1 h u +h l Z η s −(hu+h l ) u i (z)dz = σ i k i (h u +h l ) ( σ i (σ i ) 2 cosh(k i h u )−gk i sinh(k i h u ) ) a i cosϑ i (4.51) fromthelineartheory. Usingtherelationshipbetweena i anda s (Lamb(1932),Liu(2006)) a s = σ i (σ i ) 2 cosh(k i h u )−gk i sinh(k i h u ) a i , (4.52) Equation 4.51 can be further simplified to, u i = σ i k i (h u +h l ) η s (4.53) As only linear internal waves are examined, any non-linear correction for large am- plitude waves is neglected. With the above solution, however, it can be deduced that two-layer internal waves produce a mean instantaneous flux which excites linear sur- face motion. The above expression for depth-averaged internal velocity, when using in the derived model, will produce locked surface models that should reproduce the linear, dispersive solution above exactly. This can be seen using the linearized version of Equa- tion G.1. Assuming that water depth is constant, density has no gradient in horizontal plane, and only the internal-mode free surface wave exists, whose dispersion property is locked to the internal waves, Equaton 4.53 can be exactly obtained. Also note that for 132 other types of internal wave structures(e.g, solitary waves), u i can be straightforwardly obtained and applied to Equation 4.22, but will not be derived here. Now, surface movement excitation due to internal motion will be examined with experimental and numerical data. For laboratory data, Umeyama(2002) is used here. He used a 3m-long, 0.15m-wide and 0.22m-deep wave tank equipped with an oil-pressure- type wave maker on one end and wave-absorber on the other to generate internal waves. The density profile measured in the experiment is recreated and shown in Figure 4.7, with approximated value by Equation 4.28 for the comparison. −0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 −0.25 −0.2 −0.15 −0.1 −0.05 0 ρ l −ρ ρ u z(m) Figure 4.7: Measured(solid) and approximately calculated(dashed) density profile in Umeyama(2002) 133 Table 4.3: Experiment conditions of internal waves in two-layer(Umeyama(2002)) CASE I CASE II Wave amplitude(a i ) 0.005 m 0.005 m Wave period(T i ) 5.4 sec 6.4 sec Low layer thickness(h l ) 0.15m 0.15m Upper layer thickness(h u ) 0.15m 0.15m Pycnocline thickness(δ) 0.006m 0.006m Low layer density (ρ l ) 1050 kg/m 3 1050 kg/m 3 Upper layer density(ρ u ) 1000 kg/m 3 1000 kg/m 3 Sincethepycnoclinethicknessisfoundtobeextremelysmall(δ =δ ′ /h ′ =0.006m/0.30m= 0.02), the R v P term, which is O(µ 2 δ) can be neglected. Table 4.3 summarizes test con- ditions of the experiment. It should be noted here that internal waves propagate as the internal-modewavesofofEquation4.46. Thereisnoexternal-modewaveeitheratsurface or at the interface 1 . The physical layout of the experiment is mapped to the computational domain for numerical implementation, in which Δx = 0.02m and Δt = 0.002sec have been chosen. Simulated results of the free surface elevation induced by internal waves are presented in Figure 4.8 and 4.9, including both experimental data and analytical solutions (Equa- tion 4.52) for comparison. Expected agreement between numerical and analytical results is found. Some disagreement with measurement is seen in both in magnitude and phase of the free surface waves. Considering the small scale and measurement complexity in these laboratory tests, the calculation gives reasonable results. 1 Referring to the internal/external-mode, this study follows P˘ ar˘ au and Dias(2001)’s direction. 134 0 50 100 150 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x 10 −3 Time(sec) Free Surface Displacement(m) Measured Numerical Analytic Figure 4.8: Measured(dashed), numerical(solid) and analytical(dot) surface elevation in CASE I 4.4.3 Waves Excited by Internal Motion: Nonlinear Surface Waves As explored in the previous section, surface waves are generated by internal waves ac- cordingtothelinearsolutions. Thesewavesarelockedtothephaseoftheinternalwaves, and so are linearly governed by the internal wave dispersion relation. In a more general situation where the free gravity waves, such as wind waves or tsunamis, propagate on the surface, there may exist some nonlinear coupling between these two modes. Thus, in this section, the nonlinear interaction between internal- and external-mode surface waves 135 0 50 100 150 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x 10 −3 Time(sec) Free Surface Displacement(m) Measured Numerical Analytic Figure 4.9: Measured(dashed), numerical(solid) and analytical(dot) surface elevation in CASE II are investigated. This analysis will rely on quantifying the analytical behavior of the first-order nonlinear system, as given in Section 4.3.2. Previous nonlinear similar analysis for Boussinesq-type equations can be found in, for example, Kennedy et al.(2001) and Lynett and Liu(2004), among numerous others. Stokestheoryisawell-definedandeasily-usedmethodthatcanbeappropriatelyadopted 136 here. Following Kennedy et al.(2001), the first-order-nonlinear equation model in one- horizontal dimension, with constant depth, and weak horizontal and vertical density gradients is given as ζ [2] t + h u [2] α x +h 3 1 2 z α h 2 + z α h + 1 3 u [2] α xxx = − h ζ [1] u [1] α i x −h 2 1 2 z α h 2 + z α h h ζ [1] u [1] α xx i x − h ζ [1] u i i x (4.54) u [2] α t +g ζ [2] x +h 2 1 2 z α h 2 + z α h u [2] α xxt =−u [1] α u [1] α x +h h ζ [1] u [1] α xt i x −h 2 1 2 z α h 2 + z α h h u [1] α u [1] α xx i x − 1 2 h 2 u [1] α x 2 x − u [1] α u i x (4.55) where the superscript [1] and [2] represent first and second order solutions in the Stokes system, respectively. Suppose that there are internal- and external-mode surface waves expressed as, ζ [1] =ζ e +ζ i =a e cosϑ e +a i cosϑ i (4.56) u [1] α =b e cosϑ e (4.57) u i =b i cosϑ i (4.58) 137 where superscripts i and e imply internal- and external-mode components, respectively. Again, the internal-mode surface wave has been generated by internal interfacial waves while the external-mode surface wave is assumed a free gravity wave on the air-fluid surface. Theamplitudeoffreesurfacedisplacementandvelocityaresymbolizedbyaand b. Second-order solutions are then assumed as, ζ [2] =a ee cos(2ϑ e )+a ii cos 2ϑ i +a + cos ϑ e +ϑ i +a − cos ϑ e −ϑ i (4.59) u [2] α =b ee cos(2ϑ e )+b ii cos 2ϑ i +b + cos ϑ e +ϑ i +b − cos ϑ e −ϑ i (4.60) inwhichsuperscripts ii and ee areforself-interactions, whereas + and − areforthesuper- and sub-harmonic components, respectively. The phase function is ϑ (i,e) = k (i,e) x− σ (i,e) t+ǫ (i,e) . Substituting first and second order solutions into Equation 4.54 and 4.55, and collecting super- and sub-harmonic terms, a ± can be obtained as, a ± = M ± 22 L ± 1 −M ± 12 L ± 2 M ± 11 M ± 22 −M ± 12 M ± 21 (4.61) where M ± 11 =σ e ±σ i (4.62) 138 M ± 12 =−(k e ±k i )h 1− 1 2 z α h 2 + z α h + 1 3 (k e ±k i ) 2 h 2 (4.63) M ± 21 =−g(k e ±k i ) (4.64) M ± 22 = σ e ±σ i 1− 1 2 z α h 2 + z α h (k e ±k i ) 2 h 2 (4.65) L ± 1 = 1 2 (k e ±k i ) a i b e +a e b i − 1 2 z α h 2 + z α h a i b e (k e h) 2 (4.66) L ± 2 = 1 2 (k e ±k i ) b i b e −a i b e k e hσ e (4.67) This solution may be compared with an exact potential solution of Liu(2006), who derived second-order analytical expressions for internal and surface waves in a two-layer density-stratified fluid. By manipulating Liu’s solution, the amplitude of super- and sub-harmonic waves, a ± exact can be expressed (see Appendix H). Figure 4.10 presents the comparison of super- and sub-harmonic wave amplitudes between the presented model and potential solutions when h u = 0.5,h l = 0.5 and ρ u /ρ l = 0.9. The contours show model accuracy relative to potential theory; values near 1.0 are desired. As expected with the derived, long-wave-based model, accurate results are seen at small wave numbers and errors grow as wave number increases. Both 139 super-harmonic and sub-harmonic wave amplitudes tend to be underpredicted as wave dispersion increases. From this analysis, it is concluded that nonlinear interaction be- tween internal- and external-mode surface waves can be predicted reasonably within the range of k e h<0.3 and k i h<0.7. 0.5 0.5 0.5 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 0.8 0.8 0.8 0.8 0.8 0.85 0.85 0.85 0.85 0.85 0.9 0.9 0.9 0.9 0.9 0.95 0.95 0.95 0.95 0.95 k i h Subharmonic k e h Superharmonic 0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4 Figure 4.10: Second-order super-harmonic and sub-harmonic amplitude relative to the exact solution(a ± /a ± exact ) 4.5 Summary Byallowinghorizontalandverticalvariationoffluiddensity,adepth-integratedmodelfor long gravity waves over variable density fluid has been developed in the present chapter. 140 To include density change effects, correction terms to represent horizontal and vertical density change are added to typical Boussinesq model of uniform density. In particular, two-layerfluidsystemischosenasverticaldensityvariations,whereinterfacialwaveeffects onfreesurfaceisadditionallyconsideredthroughdirectinclusionofvelocitycomponentof interfacial wave. During the derivation of equation set, internal wave structure is kept as generic form within the long internal wave approximation for general use, To numerically solvederivedequations, finitevolumeschemecoupledwithapproximatedRiemannsolver is adopted for leading order terms while cell-centered finite volume methods are used for others. Then, numerical model is applied to various type of problems in which density field is configured to vary either horizontally or vertically. For horizontal variation of fluid density, pneumatic breakwater system is simulated to assess density drop effects and fair agreement is seen between computed and measured data, conceding that current effect by bubbly flux is responsible for wave attenuation in some degree. Interfacial wave of two-layerfluidsystem,asthecaseofverticaldensityvariationistestedtoinvestigatehow the internal motions affects free surface. As Lamb(1932) analyzed, surface disturbance by internal wave is observed in numerical results, which is in complete agreement with analytical solution and thus provides an insight on how the internal motions within two- layer fluid system affects surface motions. Lastly, nonlinear interaction between external- and internal-mode surface waves are studied numerically and analytically, which gives reasonable range of wave numbers for practical application. 141 Chapter 5 Conclusions and Future Works 5.1 Conclusions This dissertation studies long wave modeling using the Boussinesq-type approach, fo- cusing on the complex physics locally-induced at the nearshore such as bottom-induced turbulence, wave-current interactions and variable-density flows. The following summa- rizes the major conclusions from each section of this dissertation. In Chapter 2, a two-way coupled model for long wave simulation has been proposed for the purpose of seamlessly modeling long wave evolution from the deep ocean to the shoreline, with fine scale resolution, and without the loss of important physics in the nearshore. Thetwocomponentsaretheshallow-watersolverCOMCOTandadispersive, turbulent, and rotational Boussinesq model. A general benchmark test has been com- pleted with various conditions provided for validity of the coupled model application. As a further validation of the coupled model, long wave propagation onto a shallow shelf has been examined and compared with laboratory data. The results demonstrate that near coastalareas,dispersiveeffectsmaybelocallyimportant. Asafinaltest,arecenttsunami 142 event, the 2004 Sumatra tsunami, has been simulated with a far field focus on the Port of Salalah. The coupled model has successfully simulated various sizes of eddies generated by the tsunami through turbulent activity, and these local physics are confirmed by the observations presented in Okal et al.(2006). In Chapter 3, a set of depth-integrated equations describing combined wave-current flows are derived mathematically and discretized numerically. The Doppler shift effect, a basicpropertyinwave-currentmodels,isconfirmedanalyticallybylimitingtheequations into a one dimensional linear set. For the formulation of velocity component in deriva- tion, additional stresses introduced by nonlinear wave-current interaction are considered in addition to the bottom-induced stress. Using a parameter b, an additional stress is defined to represent the intensity of turbulent interaction between waves and currents. An appropriate estimation ofb is provided through Kemp and Simons(1982, 1983)’s ex- periments. Long wave propagation over either depth-uniform or linearly-sheared current fields are tested for validation. Comparison of numerical calculations with measurement reveals that the rotational behavior of currents needs to be taken into account when de- scribing the coupled wave-current velocity field. As a final attempt, the model is applied toacomplexconfiguration, whichcontainsbothwave-waveandwave-currentinteraction, as well as spatial variation of background currents. Subharmonic bound long waves gen- eratedbywave-waveinteractionarewellrecreatedinthepresenceofbackgroundcurrents and the comparison with the measurements shows good agreement. In Chapter 4, by allowing horizontal and vertical variation of fluid density, a depth- integrated model for long gravity waves over a variable density fluid has been devel- oped. Throughout the derivation, density change effects appear as correction terms to a 143 Boussinesq-typemodelofuniformdensitywhileinternalwaveeffectsonfreesurfacewaves in a two-layer fluid system are accounted for through direct inclusion of the internal wave velocity field. For general use of the model, the internal wave structure is kept arbitrary. Then, a numerical model is applied to various type of problems in which the density field is configured to vary either horizontally or vertically. A pneumatic breakwater system is tested by the model to examine horizontal density drop effects, and fair agreement is seenbetweencalculatedandmeasureddata. Atwo-layerinternalwavewithweakdensity stratification is simulated to investigate how internal motions affect the free surface. Nu- mericalresultsshowcompleteagreementwithananalyticalsolution, providinginsighton how the internal motions within a two-layer fluid system affect surface motions. Lastly, nonlinear interaction between external- and internal-mode surface waves is studied nu- merically and analytically, which gives reasonable range of wave conditions for practical application. 5.2 Future Works The following are future research ideas to improve the accuracy and applicability of presentedmodels. Forfurtherunderstandingofthemixingandtransportpropertiesofthe nearshore flows, an unsteady density and current field, which should be two-way coupled with a hydrodynamic model, will be considered. With this, more practical application of the model can be expected. Regarding the numerical technique, a non-diffusive numerical scheme is in need since numerical diffusion, as we explored in the first topic, may affect the current patterns and 144 thusleadtonumericalerrors. Atthesametime,amorerobustandaccuratesolvercanbe soughtfortheapplicationofthemodeltomorechallengingproblems. Fine-resolution(less thanO(1m)) modeling at the large vicinity of coastal areas may be a valuable effort. To capture the local physics and to reduce numerical errors which are increasing with grid size, fine-scale simulations are regarded as necessary measures. To this end, present model will be parallelized by MPI for practical implementation. Lastly, comparison of the wave-current model with field data will be performed, for the rigorous validation of wave-current model. 145 Bibliography [1] Alam, M.R. (2012) A new triad resonance between co-propagating surface and inter- facial waves, Journal of Fluid Mechanics, 691, p. 267-278 [2] Armi, L. and Farmer, D.(1986) Maximal two-layer exchange through a contraction with barotropic net flow, Journal of Fluid Mechanics, 164, p. 27-51 [3] Benjamin, T. B. (1962) The solitary wave on a stream with an arbitrary distribution of vorticity, Journal of Fluid Mechanics, 12, p. 97-116. [4] Bulson, P.S.(1963) Large scale bubble breakwater experiments, Dock and Harbour Authority, London, XII(516), p. 191-197 [5] Chen,Q.(2006)FullynonlinearBoussinesq-typeequationsforwavesandcurrentsover porous beds, Journal of Engineering Mechanics, 132(2), p. 220-230. [6] Chen,Q.,Kirby,J.T.,Dalrymple,R.A.,Shi,F.andThornton,E.B.(2003)Boussinesq modeling of longshore currents, Journal of Geophysical Research, 108(C11), p. 3362- 3380. [7] Choi,W.(2003)Stronglynonlinearlonggravitywavesinuniformshearflows,Physical Review E, 68, p. 026305 [8] Choi, W. (2009) Nonlinear surface waves interacting with a linear shear current, Mathematics and Computers in Simulation, 80, p. 29-36 [9] Choi, W. and Camassa, R. (1996) Weakly nonlinear internal waves in a two-fluid system, Journal of Fluid Mechanics, 313, p. 83-103 146 [10] Choi, W. and Camassa, R. (1999) Fully nonlinear internal waves in a two-fluid system, Journal of Fluid Mechanics, 384, p. 27-58 [11] Christoffersen, J.B. and Jonsson, I.G. (1985) Bed friction and dissipation in a com- bined current and wave motion, Ocean Engineering, 12(5), p. 387-423 [12] Clift, R., Grace, J., Weber, M. (1978) Bubbles, drops and particles, Academic Press. [13] Craig, W., Guyenne, P. and Sulem, C. (2011) Coupling between internal and surface waves, Natural Hazards, 57, p. 617-642 [14] Davies, A.G., Soulsby, R.L. and King, H.L. (1988) A numerical model of the com- bined wave and current bottom boundary layer, Journal of Geophysical Research, 93(C1), p. 491-508 [15] Dean, R.G. and Dalrymple, R.A. (1984) Water Wave Mechanics for Engineers and Scientists, Prentice-Hall [16] Dingemans, M. W., Van Kester, J. A. Th. M., Radder, A. C. and Uittenbogaard, R. E. (1996) The effect of the CL-vortex force in 3D wave-current interaction. In Proc. 25th Intl Conf. on Coastal Eng., Orlando, p. 4821-4832 [17] Dong, G., Ma, X., Xu, J., Ma, Y. and Wang, G. (2009) Experimental study of the transformation of bound long waves over a mild slope with ambient currents, Coastal Engineering, 56, pp. 1035-1042 [18] Donato,A.N.,Peregrine,D.H.andStocker,J.R.(1999)Thefocusingofsurfacewaves by internal waves, Journal of Fluid Mechanics, 384, p. 27-58 [19] Elachi, C. and Apel, J.R. (1976) Internal wave observatins made with an airborne synthetic aperture imaging radar, Geophysical Research Letters, 3(11), p. 647-650 [20] Elder, J.W. (1959) The dispersion of marked fluid in turbulent shear flow, Journal of Fluid Mechanics, 5, p. 544-560 147 [21] Erduran, K.S., Ilic, S., Kutija, V. (2005) Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations. Int. J. Numer. Methods Fluids, 49, p. 1213-1232 [22] Evans, W.A.B. and Ford, M.J. (1996) An integral equation approach to internal (2-layer) solitary waves, Phys. Fluids, 8, p. 2032-2047 [23] Farmer, D. and Armi, L. (1999) The generation and trapping of solitary waves over topography, Science, 283, p. 188-190 [24] Fernando, R., Guo, J. and Lin, P. (2011) Wave-current interaction at an angle 1: experiment, Journal of Hydraulic Research, 49(4), p. 424-436. [25] Freeman, N.C. and Johnson R.S.(1970) Shallow water waves on shear flows, Journal of Fluid Mechanics, 42, p. 401-409 [26] Fredsoe,J.,Andersen,K.H.andSumer,B.M.(1999)Wavepluscurrentoveraripple- covered bed, Coastal Engineering, 38, p. 177-221 [27] Gargett, A.E. and Hughes, B.A. (1972) On the interactin of surface and internal waves, J. Fluid Mech., 52, p. 179-191 [28] George, D.L. and LeVeque, R. J.,(2006) Finite volume methods and adaptive re- finement for global tsunami propagation and local inundation, Science of Tsunami Hazards, vol.24, p. 319-328 [29] Gobbi, M. F. and Kirby, J. T., (1999) Wave evolution over submerged sills: Tests of a high-order Boussinesq model, Coastal Engineering, 37, p. 57-96 [30] Grant, W.D. and Madsen, O.S. (1979) Combined wave and current interaction with a rough bottom, Journal of Geophysical Research, 84(C4), p. 1797-1808 [31] Grilli, S. T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J.,and Watts, P.(2007) Source constraints and model simulation of the December 26, 2004 Indian Ocean 148 tsunami, Journal of Waterway, Port, Coastal and Ocean Engineering, 133(6), p. 414- 428 [32] Groeneweg, J. and Klopman, G. (1998) Changes in the mean velocity profiles in the combined wave-current motion described in GLM formulation, Journal of Fluid Mechanics, 370, p. 271-296 [33] Harig, S., Chaeroni, Pranowo, W. S. and Behrens, J. (2008) Tsunami simulations on several scales: comparison of approaches with unstructured meshes and nested grids. Ocean Dynamics, 58, p. 429-440 [34] Helfrich,K.andMelville,W.(2006)Longnonlinearinternalwaves,Annu. Rev. Fluid Mech., 38, p. 395-425 [35] Hill,D.F.andFoda,M.A.(1998)Subharmonicresonanceofobliqueinterfacialwaves by a progressive surface waves, Proc. R. Soc. Lond. A, 454, p. 1129-1144 [36] Hinterberger, C., Frohlich, J. and Rodi, W. (2007) Three-dimensional and depth- averaged large-eddy simulations of some shallow water flows, Journal of Hydraulic Engineering, 133(8), p. 857-872 [37] Huang, Z. and Mei, C.C. (2003) Effects of surface waves on a turbulent current over a smooth or rough seabed, Journal of Fluid Mechanics, 497, p. 253-287 [38] Huang, Z. and Mei, C.C. (2006) Wave-induced longitudinal vortices in a shallow current, Journal of Fluid Mechanics, 551, p. 323-356 [39] Johnson, R.S. (2003) The Camassa-Holm equation for water waves moving over a shear flow, Fluid Dynamics Research, 33, p. 97-111 [40] Jonsson, I.G., Brink-Kjer, O. and Thomas, G.P. (1978) Wave action and set-down for waves on a shear current, Journal of Fluid Mechanics, 87(3), p. 401-416 [41] Joyce,T.M.(1974)Nonlinearinteractionsamongstandingsurfaceandinternalgrav- ity waves, J. Fluid Mech., 63, p. 801-825 149 [42] Kanarsk, Y. and Maderich, V. (2003) A non-hydrostatic numerical moel for calcula- tion free-surface stratified flows, Ocean Dynamics, 53, p. 176-185 [43] Kao, T., Pan, F.-S. and Renouard, D. (1985) Internal solitons on the pycnocline: generation, propagation, and shoaling and breaking over a slope, J. Fluid Mech., 159, p. 19-53 [44] Kemp, P. H. and Simons, R. R. (1982) The interaction of waves and a turbulent current: waves propagating with the current, Journal of Fluid Mechanics, 116, p. 227- 250 [45] Kemp, P. H. and Simons, R. R. (1983) The interaction of waves and a turbulent current: waves propagating against the current, Journal of Fluid Mechanics, 130, p. 73-89 [46] Kennedy, A.B., Kirby, J.T., Chen, Q. and Dalrymple, R.A. (2001) Boussinesq-type equations with improved nonlinear performance, Wave Motion, 33, p. 225-243 [47] Kim, D.-H., Lynett, P., and Socolofsky, S. (2009) A Depth-Integrated Model for Weakly Dispersive, Turbulent, and Rotational Fluid Flows, Ocean Modelling, 27(3-4), p. 198-214 [48] Kim, H., O’Connor, B.A., Park, I., Lee, Y. (2001) Modeling effect of intersection angle on near-bed flows for waves and currents, Journal of Waterway, Port, Coastal and Ocean Engineering, 127, p. 308-318 [49] Koop, C.G. and Butler, G. (1981) An investigation of internal solitary waves in a two-fluid system, J. Fluid Mech., 112, p. 225-251 [50] Kowalik, Z., Knight, W., Logan, T., Whitmore, P., (2005) Numerical modeling of the global tsunami: Indonesian tsunami of 26 December 2004, Science of Tsunami Hazards, 23(1) p. 40-56 150 [51] Lacor,C.A.,Smirnov,S.A.,Baelmans,M.(2004)Afinitevolumeformulationofcom- pact central schemes on arbitrary structured grids, Journal of Computational Physics, 198, p. 535-566 [52] Lamb, H. (1932) Hydrodynamics, Cambridge Univ. Press, New York [53] Leighton, F., Borthwick, A., and Taylor, P. (2010) 1-D numerical modelling of shal- lowflowswithvariablehorizontaldensity,Int.J.Numer.Meth.Fluids,62,p.1209-1231 [54] Lewis, J.E., Lake, B.M. and Ko, D.R.S.(1974) On the interaction of internal waves and surface gravity waves, J. Fluid Mech., 64, p. 773-800 [55] Lin, P. (2008) Numerical Modeling of Water Waves, Taylor & Francis Group [56] Liu, C.-M. (2006) Second-order random internal and surface waves in a two-fluid system, Geophy. Res. Letter, 33, L06610 [57] Liu, P. L.-F., Woo, S.B., Cho, Y.S. (1998) Computer Programs for Tsunami Propa- gation and Inundation, Technical Report, Cornell University [58] Liu, P.L.-F. and Wang, X. (2012) A multi-layer model for nonlinear internal wave propagation in shallow water, J. Fluid Mech., 695, p. 341-365 [59] Lo, J.-M. (1991) Air bubble barrier effect on neutrally buoyant objects, J. Hydraulic Res., 29(4), p. 437-455 [60] Longuet-Higgins,M.S.andStewart,R.W.(1960)Changesintheformofshortgravity waves on long waves and tidal currents, J. Fluid Mech., 8, p. 565-583 [61] Longuet-Higgins, M.S. and Stewart, R.W., (1962) Radiation stress and mass trans- port in gravity waves, with applications to ‘surf beats’, Journal of Fluid Mechanics, 13, p. 481-504 [62] Lowe, R.J., Rottman, J.W. and Linden, P.F. (2005) The non-Boussinesq lock- exchange problem. Part 1. Theory and experiments, J. Fluid Mech., 537, p. 101-124 151 [63] Lynett, P. (2006) Nearshore wave modeling with high-order Boussinesq-type equa- tions, Journal of Waterway, Port, Coastal and Ocean Engineering, 132(5) p. 348-357 [64] Lynett, P., Borrero, J., Liu, P. L.-F., and Synolakis, C.E. (2003) Field Survey and Numerical Simulations: A Review of the 1998 Papua New Guinea Tsunami, Pure and Applied Geophysics, v.160, p. 2119-2146 [65] Lynett, P. and Liu, P. L.-F. (2002a) A Numerical Study of Submarine Landslide Generated Waves and Runup, Proc. Royal Society of London A., 458, p. 2885-2910 [66] Lynett, P. and Liu, P. L.-F. (2002b) A two-dimensional, depth-integrated model for internal wave propagation, Wave Motion, 36, p. 221-240 [67] Lynett, P. and Liu, P. L.-F. (2004) A two-layer approach to wave modelling, Proc. R. Soc. Lond. A, 460, p. 2637-2669 [68] Ma, G., Shi, F. and Kirby, J. T. (2011) A polydisperse two-fluid model for surf zone bubble simulation, J. Geophys. Res., 116, C05010 [69] Ma, Y.-C. (1982) Effect of long waves on the evolution of deep-water surface gravity waves, Phys. Fluids, 25(3), p. 411-419 [70] Madsen, O.S. and Mei, C.C. (1969) The transformation of a solitary wave over an uneven bottom, Journal of Fluid Mechanics, 39, p. 781-791 [71] Matsuyama, M., Ikeno, M., Sakakiyama, T. and Takeda, T. (2007) A study of tsunami wave fission in an undistorted experiment, Pure and Applied Geophysics, 164, p. 617-631 [72] Mellor,G.L.(1991)Anequationofstatefornumericalmodelsofoceansandestuaries, J. Atm. Ocean. Tech., 8, p. 609-611 [73] Milgram, J. H. (1983) Mean flow in round bubble plumes. J. Fluid Mech., 133, p. 345-376 152 [74] Nguyen, H.Y. and Dias, F. (2008) A Boussinesq system for two-way propagation of interfacial waves, Physica D, 237, p. 2365-238 [75] Nwogu, O.G. (2009) Interaction of finite-amplitude waves with vertically sheared current fields, Journal of Fluid Mechanics, 627, p. 179-213 [76] Okada, Y. (1985) Surface deformation to shear and tensile faults in a half-space, Bull. Seismol. Soc. Am., 75, p. 1135-1154 [77] Okal, E.A., Fritz, H.M., Raad, P.E., Synolakis, C.E., Al-Shijbi,Y. and Al-Saifi, M., (2006)OmanFieldSurveyaftertheDecember2004IndianOceanTsunami,Earthquake Spectra, 22, S203-S218 [78] Olabarieta, M., Medina, R. and Castanedo, M. (2010) Effects of wave-current inter- action on the current profile, Coastal Engineering, 57, p. 643-655 [79] Phillips, O. M. (1977) Dynamics of the Upper Ocean, second ed., Cambridge Univ. Press, New York [80] Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P. (1992) Numerical Recipes in Fortran, Cambridge University Press [81] P˘ ar˘ au, E. and Dias, F. (2001) Interfacial periodic waves of permanent form with free-surface boundary conditions, J. Fluid Mech., 437, p. 325-336 [82] Saito, T. and Furumura, T., (2009) Three-dimensional simulation of tsunami gen- eration and propagation: application to Intraplate events, Journal of Geophysical Re- search, 114, B02307 [83] Sch¨ affer, H.A. (1996) Second-order wavemaker theory for irregular waves, Ocean Eng., 23(1), p. 47-88 [84] Segur, H. and Hammack, J.L. (1982) Soliton models of long internal waves, J. Fluid Mech., 118, p. 285-305 153 [85] Selezov,I.T.,Avrameko,O.V.,Gurtovyi,Y.V.andNaradovyi,V.V.(2010)Nonlinear interaction of internal and surface gravity waves in a two-layer fluid with free surface, J. Math. Sci., 168(4), p. 590-602 [86] Sharma, J.N. and Dean, R.G. (1981) Second-order directional seas and associated wave forces, Soc. Pet. Eng. J., 4, p. 129-140. [87] Shi, F., Kirby, J. T., Harris, J. C., Geiman, J. D. and Grilli, S. T., (2012) A high- order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation, Ocean Modelling, 43-44, p. 36-51 [88] Son,S.andLynett,P.(2012)Interactionofshallowwaterwaveswithweaklysheared currents of arbitrary profile, (in preparation) [89] Straub, L.G., Bowers, C.E. and Tarapore, Z.S. (1959) Experimental studies of pneu- matic and hydraulic breakwaters, Technical Paper, No. 25, Series B, St. Anthony Falls Hyd. Lab. University of Minnesota, Minneapolis, USA [90] Swan, C., Cummins, I.P. and James, R.L. (2001) An experimental study of two- dimensional surface water waves propagating on depth-varying currents. Part 1. Reg- ular waves, Journal of Fluid Mechanics, 428, p. 273-304 [91] Swan, C. and James, R.L. (2001) A simple analytical model for surface water waves on a depth-varying current, Applied Ocean Research, 22, p. 331-347 [92] Swan,C.(1990)Anexperimentalstudyofwavesonastronglyshearedcurrentprofile. In Proc. 22nd Conf. on Coastal Engng., Delft, The Netherlands, p. 489-502 [93] Shapiro,R.(1970)Smoothing,filtering,andboundaryeffects,Reviews of Geophysics and Space Physics, 8(2), p. 359-387 [94] Shen,C.Y.(2001)ConstituentBoussinesqequationsforwavesandcurrents,Journal of Physical Oceanography, 31, p.850-859 154 [95] Smith, J.A. (2006) Observed variability of ocean wave Stokes drift, and the Eulerian response to passing groups, Journal of Physical Oceanography, 36, p. 1381-1402 [96] Taylor, Sir Geoffrey (1955) The action of a surface current used as a breakwater, Proc. Roy. Soc. Lond. A, 231, p. 446-478 [97] Titov, V.V., and Synolakis, C.E. (1998) Numerical modeling of tidal wave runup, Journal of Waterway, Port, Coastal and Ocean Engineering, 124(4), p. 157-171 [98] Thomas, G.P. (1981) Wave-current interactions: an experimental and numerical study. Part 1. Linear waves, Journal of Fluid Mechanics, 110, p. 457-474 [99] Thomas, G.P. (1990) Wave-current interactions: an experimental and numerical study. Part 2. Nonlinear waves, Journal of Fluid Mechanics, 216, p. 505-536 [100] Tonelli, M. and Petti, M. (2009) Hybrid finite volume-finite difference scheme for 2DH improved Boussinesq equations. Coastal Eng., 56, p. 609-620 [101] Toro, E. F. (2002) Shock-capturing methods for free-surface shallow flows, Wiley, Chichester, UK [102] Umeyama, M. (2002) Experimental and theoretical analyses of internal waves of finite amplitude, J. Waterway, Port, Coastal and Ocean Eng., 128(3), p. 133-141 [103] Umeyama, M.(2005)Reynoldsstressesandvelocitydistributionsinawave-current coexisting environment, Journal of Waterway, Port, Coastal and Ocean Engineering, 131(5), p. 203-212 [104] Vinokur, M. (1974) Conservation equations of gas dynamics in curvilinear coordi- nate systems, Journal of Computational Physics, 14(2), p. 105-125 [105] Wang, X., Liu, P. L.-F. (2006) An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami, Journal of Hydraulic Research, 44(2), p. 147- 154 155 [106] Wei, G., Kirby, J. T. (1995) A time-dependent numerical code for extended Boussi- nesq equations, Journal of Waterway, Port, Coastal and Ocean Engineering, 121, p. 251-261 [107] Wei, G., Kirby, J. T., Grilli, S. T. and Subramanya, R. (1995) A fully nonlinear Boussinesqmodelforsurfacewaves.I.Highlynonlinearunsteadywaves.J.FluidMech., 294, p. 71-92 [108] Wood, I.R. (1970) A lock exchange flow, J. Fluid Mech., 42, p. 671-687 [109] Yamamoto, S. and Daiguji, H. (1993) Higher-order-accurate upwind schemes for solving the compressible Euler and Navier-Stokes equations, Comput. Fluids, 22, p. 259-270 [110] Yang, S.Q., Tan, S.K., Lim, S.Y. and Zhang, S.F. (2006) Velocity distribution in combined wave-current flows, Advances in Water Resources, 29(8), p. 1196-1208 [111] Yao, A. and Wu, C.H.(2005) Incipient breaking of unsteady waves on sheared cur- rents. Physics of Fluids, 17, p. 082104 [112] Yoon,S.B.,(2002)Propagationofdistanttsunamisoverslowlyvaryingtopography. Journal of Geophysical Research, 107(C10), 3140 [113] You, Z. (1996) The effect of wave-induced stress on current profiles, Ocean Engi- neering, 23(7), p. 619-628 [114] Zhang,C.-X.,Wang,Y.-X.,Wang,G.-Y.andYu,L.-M.etal.(2010)Wavedissipat- ing performance of air bubble breakwaters with different layouts, J. Hydrodynamics, 22(5), p. 671-680 [115] Zhou, J.G., Causon, D.M., Mingham, C.G. and Ingram, D.M. (2001) The surface gradient method for the treatment of source terms in the shallow water equations, Journal of Computational Physics, 168(2), p. 1-25 156 Appendix A Numerical Scheme of COMCOT The numerical solutionscheme employed by COMCOT isthe explicitleap-frogdifference method. Nonlinear terms in the model are approximated with upwind finite differences. The final forms for the continuity and momentum equations are described below: ζ n+1/2 i,j =ζ n−1/2 i,j − Δt Δx M n i+1/2,j −M n i−1/2,j − Δt Δy N n i,j+1/2 −N n i,j−1/2 (A.1) M n+1 i+1/2,j = 1 1+ν x Δt (1−ν x Δt)M n i+1/2,j −g Δt Δx H n+1/2 i+1/2,j ζ n+1/2 i+1,j −ζ n+1/2 i,j − Δt Δx(1+ν x Δt)   λ 11 M n i+3/2,j 2 H n i+3/2,j +λ 12 M n i+1/2,j 2 H n i+1/2,j +λ 13 M n i−1/2,j 2 H n i−1/2,j    (A.2) − Δt Δy(1+ν x Δt) " λ 21 (MN) n i+1/2,j+1 H n i+1/2,j+1 +λ 22 (MN) n i+1/2,j H n i+1/2,j +λ 23 (MN) n i+1/2,j−1 H n i+1/2,j−1 # 157 N n+1 i,j+1/2 = 1 1+ν y Δt (1−ν y Δt)N n i,j+1/2 −g Δt Δy H n+1/2 i,j+1/2 ζ n+1/2 i,j+1 −ζ n+1/2 i,j − Δt Δx(1+ν y Δt) " λ 31 (MN) n i+1,j+1/2 H n i+1,j+1/2 +λ 32 (MN) n i,j+1/2 H n i,j+1/2 +λ 33 (MN) n i−1,j+1/2 H n i−1,j+1/2 # − Δt Δy(1+ν y Δt)   λ 41 N n i,j+3/2 2 H n i,j+3/2 +λ 42 N n i,j+1/2 2 H n i,j+1/2 +λ 43 N n i,j−1/2 2 H n i,j−1/2    (A.3) where the coefficients of the upwind scheme are obtained by    λ 11 =0, λ 12 =1, λ 13 =−1, if M n i+1/2,j ≥0 λ 11 =1, λ 12 =−1, λ 13 =0, if M n i+1/2,j <0    λ 21 =0, λ 22 =1, λ 23 =−1, if N n i+1/2,j ≥0 λ 21 =1, λ 22 =−1, λ 23 =0, if N n i+1/2,j <0    λ 31 =0, λ 32 =1, λ 33 =−1, if M n i,j+1/2 ≥0 λ 31 =1, λ 32 =−1, λ 33 =0, if M n i,j+1/2 <0    λ 41 =0, λ 42 =1, λ 43 =−1, if N n i,j+1/2 ≥0 λ 41 =1, λ 42 =−1, λ 43 =0, if N n i,j+1/2 <0 Bottom friction terms are given as ν x = 1 2 gm 2 H n i+1/2,j 7/3 M n i+1/2,j 2 + N n i+1/2,j 2 1/2 (A.4) ν y = 1 2 gm 2 H n i,j+1/2 7/3 M n i,j+1/2 2 + N n i,j+1/2 2 1/2 (A.5) 158 Appendix B Second order terms in Boussinesq Equation D c =−∇· " (ζ +h) ( ζ 2 −ζh+h 2 6 − z 2 α 2 ! ∇S+ (ζ−h) 2 −z α ∇T )# +∇· " ψ(ζ +h) ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# (B.1) (D x m ,D y m )=D+D ν +ξ+ξ ν −∇· ν h t ∇U α +ν v t ∇S+ τ b ρH (B.2) D= 1 2 ∇ z 2 α U α ·∇S +∇(z α U α ·∇T)+(T∇T) − 1 2 ∇ ζ 2 ∂S ∂t −∇ ζ ∂T ∂t + 1 2 z 2 α ∂∇S ∂t +z α ∂∇T ∂t − 1 2 ∇ ζ 2 U α ·∇S −∇(ζU α ·∇T)+∇ 1 2 ζ 2 S 2 +∇(ζTS) (B.3) 159 D ν = (ζ−h) 2 ∂ψζ ∂t − ζ 2 −ζh+h 2 6 ∂ψ ∂t + ∂ ∂t ψ z 2 α 2 −ζz α + (ζ−h) 2 ∇{U α ·(ψζ)}− ζ 2 −ζh+h 2 6 ∇(U α ·ψ) +∇ U α · ψ z 2 α 2 −ζz α −ψ ( ζ 2 +ζh−2h 2 S 6 + (ζ +h)T 2 ) (B.4) where∇ = (∂/∂x,∂/∂y), S =∇·U α , T =∇·(hU α ) andψ =τ b /{ρν v t (ζ +h)}. Also ξ = ξ x ,ξ y andξ ν = ξ ν x ,ξ ν y are given by ξ x = −V α ∂z α ∂x z α ∂S ∂y + ∂T ∂y − ∂z α ∂y z α ∂S ∂x + ∂T ∂x (B.5) − ∂V α ∂x − ∂U α ∂y "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) ∂S ∂y + z α − (ζ−h) 2 ∂T ∂y # ξ y = U α ∂z α ∂x z α ∂S ∂y + ∂T ∂y − ∂z α ∂y z α ∂S ∂x + ∂T ∂x (B.6) + ∂V α ∂x − ∂U α ∂y "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) ∂S ∂x + z α − (ζ−h) 2 ∂T ∂x # ξ ν x = −V α " ∂ ∂x ψ y 1 2 z 2 α −z α ζ − ζ 2 −ζh+h 2 6 ∂ψ y ∂x + (ζ−h) 2 ∂ψ y ζ ∂x − ∂ ∂y ψ x 1 2 z 2 α −z α ζ + ζ 2 −ζh+h 2 6 ∂ψ x ∂y − (ζ−h) 2 ∂ψ x ζ ∂y # − ∂V α ∂x − ∂U α ∂y ψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) (B.7) 160 ξ ν y = U α " ∂ ∂x ψ y 1 2 z 2 α −z α ζ − ζ 2 −ζh+h 2 6 ∂ψ y ∂x + (ζ−h) 2 ∂ψ y ζ ∂x − ∂ ∂y ψ x 1 2 z 2 α −z α ζ + ζ 2 −ζh+h 2 6 ∂ψ x ∂y − (ζ−h) 2 ∂ψ x ζ ∂y # + ∂V α ∂x − ∂U α ∂y ψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 ) (B.8) and (ψ x ,ψ y )=ψ. 161 Appendix C Variables in Numerical Scheme of Boussinesq Model E =E LO +E D +E V (C.1) F =F LO +F D +U α (E D +E V ) (C.2) G=G LO +G D +V α (E D +E V ) (C.3) E LO , F LO , and G LO are rewritten by E LO =− ∂HU α ∂x − ∂HV α ∂y (C.4) F LO =− ∂ ∂x HU 2 α + 1 2 gH 2 − ∂HU α V α ∂y +gH ∂h ∂x (C.5) G LO =− ∂HU α V α ∂x − ∂ ∂y HV 2 α + 1 2 gH 2 +gH ∂h ∂y (C.6) and E D , E V , F D , G D , F 1 and G 1 are defined as 162 E D = H 1 6 ζ 2 −ζh+h 2 − 1 2 z 2 α ∇S+ 1 2 (ζ−h)−z α ∇T x + H 1 6 ζ 2 −ζh+h 2 − 1 2 z 2 α ∇S+ 1 2 (ζ−h)−z α ∇T y (C.7) E V = − " Hψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# x − " Hψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# y (C.8) (F D ,G D ) = H 1 2 ∇ ζ 2 U α ·∇S +∇(ζU α ·∇T)− 1 2 ∇ ζ 2 S 2 − 1 2 ∇ z α 2 U α ·∇S −∇(z α U α ·∇T)−(T∇T)−∇(ζTS) − (ζ−h) 2 ∇{U α ·(ψζ)}+ ζ 2 −ζh+h 2 6 ∇(U α ·ψ) − ∇ U α · ψ z 2 α 2 −ζz α + ψ ( ζ 2 +ζh−2h 2 S 6 + HT 2 ) −ξ−ξ ν + ∇· ν h t ∇U α −ν v t ∇S− τ b ρH (C.9) F 1 = H 2 ζ 2 −z 2 α V αxy −H(z α −ζ)(hV α ) xy +Hζ x n ζV αy +(hV α ) y o (C.10) G 1 = H 2 ζ 2 −z 2 α U αxy −H(z α −ζ)(hU α ) xy +Hζ y {ζU αx +(hU α ) x } (C.11) 163 F p v , G p v , F c v and G c v are rewritten by F p v = H n ζ 2 −ζh+h 2 +3z 2 α n 6 n 2(ψ x ) n −3(ψ x ) n−1 +(ψ x ) n−2 o − H n (ζ−h−2z α ) n 2 n 2(ψ x ζ) n −3(ψ x ζ) n−1 +(ψ x ζ) n−2 o (C.12) G p v = H n ζ 2 −ζh+h 2 +3z 2 α n 6 n 2(ψ y ) n −3(ψ y ) n−1 +(ψ y ) n−2 o − H n (ζ−h−2z α ) n 2 n 2(ψ y ζ) n −3(ψ y ζ) n−1 +(ψ y ζ) n−2 o (C.13) F c v = H n+1 ζ 2 −ζh+h 2 +3z 2 α n+1 6 n (ψ x ) n+1 −(ψ x ) n o − H n+1 (ζ−h−2z α ) n+1 2 n (ψ x ζ) n+1 −(ψ x ζ) n o (C.14) G c v = H n+1 ζ 2 −ζh+h 2 +3z 2 α n+1 6 n (ψ y ) n+1 −(ψ y ) n o − H n+1 (ζ−h−2z α ) n+1 2 n (ψ y ζ) n+1 −(ψ y ζ) n o (C.15) 164 Appendix D Derivation of Momentum Equation in Wave-Current Model Vertical momentum equation, Equation 3.7 with the use of Equation 3.25 will give, µ 2 {(w 0 ) t +u α ·∇w 0 +w 0 (w 0 ) z }+p z +1=O µ 4 ,αµ 3 ,βµ 3 (D.1) Integrating this equation in z direction from z to ζ then gives an vertical expression of pressure field including hydrodynamic terms as, p = ζ−z +µ 2 1 2 z 2 −ζ 2 S t +(z−ζ)T t + 1 2 z 2 −ζ 2 u α ·∇S +(z−ζ)u α ·∇T − 1 2 z 2 −ζ 2 S 2 −(z−ζ)TS +O µ 4 ,αµ 3 ,βµ 3 (D.2) 165 Note that we used Equation 3.15. To derive a horizontal momentum equation of depth-integrated form, we first substitute Equation 3.25 and D.2 into Equation 3.6. Accordingly, unsteady term can be expressed as, u t =(u α ) t + µ 2 1 2 z 2 α −z 2 ∇S+(z α −z)∇T t + µ 2 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) t + µ 2 Ψ b − 1 2 z 2 α −z 2 +h(z−z α ) t + O(µ 4 ) (D.3) Horizontal convection term is given as, u·∇u = u α ·∇u α +µ 2 ∇ u α · 1 2 z 2 α −z 2 ∇S+(z α −z)∇T +µ 2 ∇ u α ·Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) +µ 2 ∇ u α ·Ψ b − 1 2 z 2 α −z 2 +h(z−z α ) +µ 2 ∇[u Wα ·{u C −u Cα }] +µ 2 ξ p +ξ t +ξ c +O(µ 4 ) (D.4) by virtue of vector identity,∇(A·B)=A·∇B+B·∇A+A×(∇×B)+B×(∇×A), which spawnsξ p ,ξ t andξ c terms defined respectively as, ξ p = − u α × ∇× 1 2 z 2 α −z 2 ∇S+(z α −z)∇T − 1 2 z 2 α −z 2 ∇S+(z α −z)∇T ×(∇×u α ) (D.5) 166 ξ t = − u α ×(∇×Ψ) 1 2 z 2 α −z 2 +ζ(z−z α ) − Ψ×(∇×u α ) 1 2 z 2 α −z 2 +ζ(z−z α ) − u α ×(∇×Ψ b ) − 1 2 z 2 α −z 2 +h(z−z α ) − Ψ b ×(∇×u α ) − 1 2 z 2 α −z 2 +h(z−z α ) (D.6) ξ c =−u Wα ×{∇×(u C −u Cα )}−(u C −u Cα )×(∇×u Wα ) (D.7) Ofnote, thesetermsareassociatedtoverticalvorticitywhichcanbeincurredfromeither horizontally rotational or irrotational flows(Chen et al.(2003)). Vertical convective term becomes, w(u) z =µ 2 z 2 S∇S+zT∇S+zS∇T +T∇T +w 0 Ω 1 +O(µ 4 ) (D.8) And pressure gradient term on the right hand side of Equation 3.6 produces, ∇p = ∇ζ + µ 2 1 2 z 2 ∇S t − 1 2 ∇ ζ 2 S t +z∇T t −∇(ζT t )+ 1 2 z 2 ∇(u α ·∇S) − 1 2 ∇ ζ 2 u α ·∇S +z∇(u α ·∇T)−∇(ζu α ·∇T)− 1 2 z 2 ∇S 2 + 1 2 ∇ ζ 2 S 2 −z∇(TS)+∇(ζTS) + O(µ 4 ,αµ 3 ,βµ 3 ) (D.9) Then, horizontal and vertical viscosity terms are written, 167 αµ ∇· ν h t ∇u =αµ ∇· ν h t ∇u α +O(αµ 3 ) (D.10) β µ (ν v t u z ) z = −βµν v t ∇S−βµ τ b ρ(ζ +h) +βµ τ bb ρ(ζ +h) + βµν v t (u C ) zz +O(βµ 3 ) (D.11) Addingupallthesecomponentssubsequentlymakesthedepth-integratedmomentum equations as, (u α ) t + u α ·∇u α +∇ζ + µ 2 1 2 z 2 α ∇S t +µ 2 z α ∇T t −µ 2 1 2 ∇ ζ 2 S t −µ 2 ∇(ζT t )+µ 2 T∇T + µ 2 1 2 ∇ z 2 α u α ·∇S +µ 2 ∇(z α u α ·∇T)+µ 2 1 2 ∇ ζ 2 S 2 − µ 2 1 2 ∇ ζ 2 u α ·∇S −µ 2 ∇(ζu α ·∇T)+µ 2 ∇(ζTS) + µ 2 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) t + µ 2 Ψ b − 1 2 z 2 α −z 2 +h(z−z α ) t + µ 2 ∇ u α ·Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) + µ 2 ∇ u α ·Ψ b − 1 2 z 2 α −z 2 +h(z−z α ) + µ 2 ∇{u Wα ·(u C −u Cα )}+µ 2 (w 0 Ω 1 )+µ 2 ξ p +ξ t +ξ c − αµ ∇· ν h t ∇u α +βµν v t ∇S+βµ τ b ρ(ζ +h) −βµ τ bb ρ(ζ +h) − βµν v t (u C ) zz = O(µ 4 ,αµ 3 ,βµ 3 ) (D.12) 168 Now our concern turns to elimination of z dependency, as a final effort. Among others, depth-averaging procedure proposed by Chen(2006) is employed to get the depth- integratedmomentumequationwithoutanyzdependency. Exceptdepth-invariantterms, i.e., terms from irrotational velocity, all others including ξ p should be depth-averaged defined as, for example, 1 ζ +h Z ζ −h z 2 dz = 1 3 ζ 2 −ζh+h 2 (D.13) The final equation is provided in the main body of this chapter 169 Appendix E Numerical Scheme in Wave-Current Model The system of equations are discretized to develop numerical solutions. All the spatial derivative terms are discretized by using a conservative finite volume method, while the temporal derivative terms by third-order Adams-Bashforth predictor and the fourth- order Adams-Moulton corrector scheme. For the application of finite volume method, governing equations are in a conservative form to effectively represent conservation laws. Subsequently, a complete set of equations in a conservative form can be obtained by equating each term using an assumption of h t =0. (ζ +h) t +{(ζ +h)u α } x +{(ζ +h)v α } y +C =0 (E.1) {(ζ +h)u α } t + (ζ +h)u 2 α + 1 2 g(ζ +h) 2 x +{(ζ +h)u α v α } y − g(ζ +h)h x +(ζ +h)M x +u α C =0 (E.2) {(ζ +h)v α } t + {(ζ +h)u α v α } x + (ζ +h)v 2 α + 1 2 g(ζ +h) 2 y − g(ζ +h)h y +(ζ +h)M y +v α C =0 (E.3) 170 where C = N p +N t +N c (E.4) (M x ,M y ) = H p +H t +H c +ξ p +ξ t +ξ c − ∇· ν h t ∇u α +ν v t ∇S+ τ b (ζ +h)ρ − τ bb (ζ +h)ρ − ν v t {(u C ) z | z=ζ −(u C ) z | z=−h } (E.5) Note that all physical parameters are valued in units without ∗ . Strictly speaking, above governing system is not yet casted in strong conservative form(thus weak-conservative form) due to the bottom slope term as well as higher order terms, which are thought of assourcesorsinksintheequations. Thepresenceofsuchtermspreventstheachievement of overall conservation of mass and momentum, which may thus lead to inferior shock capturing capability in regions of very large shock curvature(Vinokur(1974)). Bottom slope term of leading order, therefore, needs some special treatment to physically and numerically satisfy the conservation law. In this study, surface gradient method intro- duced by Zhou et al.(2001) is employed to preserve conservative property of mass and momentum even in stationary flow condition. E.1 Finite Volume Method GodunovschemecoupledwithaRiemannsolverisoneofthemostsuccessfultechniqueto solve hyperbolic system of conservation laws. As a practical advantage for hydrodynamic modeling, it is capable to capture the shock-like discontinuities such as bores by wave breaking or hydraulic jumps by abrupt bathymetric changes. Since the higher order discretization is required for high degree of accuracy, the high frequency oscillation errors 171 due to over- or under-shoot are anticipated when using a classic Godunov scheme. So- calledGodunovtheoremprovedthatlinearschemecannotbefreeofspuriousoscillations to get higher order accuracy and this urges the needs for TVD treatment by means of reconstructing. Therefore, fourth-order compact MUSCL-TVD(Monotone Upstream- centeredSchemeforConservationLaw-TVD)schemeisusedinthisstudytosolveleading order conservative part of the governing equations(Kim et al.(2009)). Reconstructed values at each interface of piecewise volume, which usually are limited by slope limiter functions, can be then solved by approximated HLL Riemann solver(Toro (2002)). On the other hand, higer-order terms are discretized by a cell averaged finite volume method given by Lacor et al.(2003). E.2 Time Marching The governing equations are marched through time by a third-order Adams-Bashforth predictor and a fourth-order Adams-Moulton corrector scheme. At the predictor step, values at new time step (n + 1) are evaluated explicitly by known values at three time steps, (n−2), (n−1) and n, which can be described as ζ n+1 =ζ n + Δt 12 23E n −16E n−1 +5E n−2 (E.6) P n+1 = P n + Δt 12 23F n −16F n−1 +5F n−2 +2F n 1 −3F n−1 1 +F n−2 1 +F p v (E.7) Q n+1 = Q n + Δt 12 23G n −16G n−1 +5G n−2 +2G n 1 −3G n−1 1 +G n−2 1 +G p v (E.8) where the superscript n denotes time step. And P, Q are defined as 172 P =(ζ +h)u α + (ζ +h) 2 z 2 α −ζ 2 (u α ) xx +(ζ +h)(z α −ζ)(hu α ) xx − (ζ +h)ζ x {ζ(u α ) x +(hu α ) x } (E.9) Q=(ζ +h)v α + (ζ +h) 2 z 2 α −ζ 2 (v α ) yy +(ζ +h)(z α −ζ)(hv α ) yy − (ζ +h)ζ y n ζ(v α ) y +(hv α ) y o (E.10) E, F, G, F 1 , G 1 , F p v , G p v , F c v , G c v in the above equations are defined respectively, as E =E LO +E HO (E.11) F =F LO +F HO +u α E HO (E.12) G=G LO +G HO +v α E HO (E.13) E LO , F LO , and G LO are rewritten by E LO =−{(h+ζ)u α } x −{(h+ζ)v α } y (E.14) F LO =− (h+ζ)u 2 α + 1 2 g(h+ζ) 2 x −{(h+ζ)u α v α } y +g(h+ζ)h x (E.15) G LO =−{(h+ζ)u α v α } x − (h+ζ)v 2 α + 1 2 g(h+ζ) 2 y +g(h+ζ)h y (E.16) 173 and E HO , F HO , G HO , F 1 and G 1 are defined as E HO = (h+ζ) ζ 2 −ζh+h 2 6 − 1 2 z 2 α ∇S+ ζ−h 2 −z α ∇T x + (h+ζ) ζ 2 −ζh+h 2 6 − 1 2 z 2 α ∇S+ ζ−h 2 −z α ∇T y − " (h+ζ)ψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# x − " (h+ζ)ψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# y + " (h+ζ)ψ x b ( z 2 α 2 +z α ζ− 2ζ 2 +2ζh−2h 2 6 )# x + " (h+ζ)ψ y b ( z 2 α 2 +z α ζ− 2ζ 2 +2ζh−2h 2 6 )# y − Z ζ −h u C dz−(h+ζ)u Cα x − Z ζ −h v C dz−(h+ζ)v Cα y (E.17) 174 (F HO ,G HO ) = (h+ζ) 1 2 ∇ ζ 2 u α ·∇S +∇(ζu α ·∇T)− 1 2 ∇ ζ 2 S 2 − 1 2 ∇ z 2 α u α ·∇S −∇(z α u α ·∇T)−∇(ζTS) − (T∇T)−∇{E(ζS+T)}−ξ p i − E(ζS+T)∇ζ− 1 2 ζ 2 −z 2 α E∇S−(ζ−z α )E∇T + (h+ζ) " ζ 2 −ζh+h 2 6 ∇(u α ·Ψ)− (ζ−h) 2 ∇{u α ·(Ψζ)} + Ψ ( ζ 2 +ζh−2h 2 S 6 + (h+ζ)T 2 ) − ∇ u α · Ψ z 2 α 2 −ζz α −ξ t + (h+ζ) "( z 2 α 2 − ζ 2 −ζh+h 2 6 ) ∇(u α ·Ψ b ) − (ζ−h) 2 −z α ∇{u α ·(Ψ b h)} + Ψ b ( 2ζ 2 +ζh−h 2 S 6 + (h+ζ)T 2 )# + (ζ +h) n ∇(u Wα ·u Cα )−ξ c o −∇ u Wα · Z ζ −h u C dz + (ζS+T)u C | z=ζ +(hS−T)u C | z=−h −S Z ζ −h u C dz + (h+ζ) ∇· ν h t ∇u α −ν v t ∇S− τ b −τ bb (ζ +h)ρ + ν v t {(u C ) z | z=ζ −(u C ) z | z=−h }] (E.18) F 1 = 1 2 (h+ζ) ζ 2 −z 2 α (v α ) xy −(h+ζ)(z α −ζ)(hv α ) xy + (h+ζ)ζ x n ζ(v α ) y +(hv α ) y o (E.19) 175 G 1 = 1 2 (h+ζ) ζ 2 −z 2 α (u α ) xy −(h+ζ)(z α −ζ)(hu α ) xy + (h+ζ)ζ y {ζ(u α ) x +(hu α ) x } (E.20) F p v , G p v , F c v and G c v are rewritten by F p v = (h+ζ) n ζ 2 −ζh+h 2 −3z 2 α n 6 n 2(ψ x ) n −3(ψ x ) n−1 +(ψ x ) n−2 o − (h+ζ) n (ζ−h−2z α ) n 2 n 2(ψ x ζ) n −3(ψ x ζ) n−1 +(ψ x ζ) n−2 o + (h+ζ) n 3z 2 a +6hz α −ζ 2 −2ζh+2h 2 n 6 n 2(ψ x b ) n −3(ψ x b ) n−1 +(ψ x b ) n−2 o (E.21) G p v = (h+ζ) n ζ 2 −ζh+h 2 −3z 2 α n 6 n 2(ψ y ) n −3(ψ y ) n−1 +(ψ y ) n−2 o − (h+ζ) n (ζ−h−2z α ) n 2 n 2(ψ y ζ) n −3(ψ y ζ) n−1 +(ψ y ζ) n−2 o + (h+ζ) n 3z 2 a +6hz α −ζ 2 −2ζh+2h 2 n 6 n 2 ψ y b n −3 ψ y b n−1 + ψ y b n−2 o (E.22) F c v = (h+ζ) n+1 ζ 2 −ζh+h 2 −3z 2 α n+1 6 n (ψ x ) n+1 −(ψ x ) n o − (h+ζ) n+1 (ζ−h−2z α ) n+1 2 n (ψ x ζ) n+1 −(ψ x ζ) n o + (h+ζ) n+1 3z 2 a +6hz α −ζ 2 −2ζh+2h 2 n+1 6 n (ψ x b ) n+1 −(ψ x b ) n o (E.23) 176 G c v = (h+ζ) n+1 ζ 2 −ζh+h 2 −3z 2 α n+1 6 n (ψ y ) n+1 −(ψ y ) n o − (h+ζ) n+1 (ζ−h−2z α ) n+1 2 n (ψ y ζ) n+1 −(ψ y ζ) n o + (h+ζ) n+1 3z 2 a +6hz α −ζ 2 −2ζh+2h 2 n+1 6 n ψ y b n+1 − ψ y b n o (E.24) With values at all four time steps including (n + 1), whose values are resulted by predictor step, another new values at time step (n + 1) are obtained through implicit corrector step. ζ n+1 =ζ n + Δt 24 9E n+1 +19E n −5E n−1 +E n−2 (E.25) P n+1 = P n + Δt 24 9F n+1 +19F n −5F n−1 +F n−2 +F n+1 1 −F n 1 +F c v (E.26) Q n+1 = Q n + Δt 24 9G n+1 +19G n −5G n−1 +G n−2 +G n+1 1 −G n 1 +G c v (E.27) Once this interative step ensures the relative error between the values at predictor and at corrector is less then 0.0001, the velocity components,u α andv α can be extracted from P and Q using tridiagonal matrix, which come from the discretized version of Equation E.9 and E.10, 177 P = (u α ) i−1,j (ζ +h i,j ) z 2 α −ζ 2 2Δx 2 + (z α −ζ)h i−1,j Δx 2 + ζ x ζ 2Δx + ζ x h i−1,j 2Δx + (u α ) i,j (ζ +h i,j ) 1− z 2 α −ζ 2 Δx 2 − 2(z α −ζ)h i,j Δx 2 (E.28) + (u α ) i+1,j (ζ +h i,j ) z 2 α −ζ 2 2Δx 2 + (z α −ζ)h i+1,j Δx 2 − ζ x ζ 2Δx − ζ x h i−1,j 2Δx Q = (v α ) i,j−1 (ζ +h i,j ) z 2 α −ζ 2 2Δy 2 + (z α −ζ)h i,j−1 Δy 2 + ζ y ζ 2Δy + ζ y h i,j−1 2Δy + (v α ) i,j (ζ +h i,j ) 1− z 2 α −ζ 2 Δy 2 − 2(z α −ζ)h i,j Δy 2 (E.29) + (v α ) i,j+1 (ζ +h i,j ) z 2 α −ζ 2 2Δy 2 + (z α −ζ)h i,j+1 Δy 2 − ζ y ζ 2Δy − ζ y h i,j−1 2Δy where subscript (i, j) identifies cell location. 178 Appendix F Derivation of Momentum Equation of Boussinesq Model for Variable Density Fluid Flows From the vertical momentum equation, the pressure field will be extracted. Substituting Equation 4.22 into Equation 4.5 yields µ 2 ρ 0 {(w 0 ) t +u α ·∇w 0 +w 0 (w 0 ) z }+p z +ρ 0 1−γtanh z−z 0 δ =O µ 4 ,αµ 3 ,βµ 3 ,γµ 2 (F.1) Integrating the equation above with respect to z provides the expression for pressure, p = ρ 0 (ζ−z)−γδρ 0 ln cosh ζ−z 0 δ −ln cosh z−z 0 δ +µ 2 ρ 0 1 2 z 2 −ζ 2 S t +(z−ζ)T t + 1 2 z 2 −ζ 2 u α ·∇S +(z−ζ)u α ·∇T − 1 2 z 2 −ζ 2 S 2 −(z−ζ)TS +O µ 4 ,αµ 3 ,βµ 3 ,γµ 2 (F.2) It is noted that Equation 4.12 has also be applied in the above equation. 179 To derive a depth-integrated momentum equation foru α , Equation 4.22 and F.2 can be applied to Equation 4.4 where each term is written as, ρu t = ρ 0 (u α ) t −γtanh z−z 0 δ ρ 0 (u α ) t +µ 2 ρ 0 u i t +µ 2 ρ 0 1 2 z 2 α −z 2 ∇S+(z α −z)∇T t +µ 2 ρ 0 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) t +O(µ 4 ,γµ 2 ) (F.3) ρu·∇u = ρ 0 u α ·∇u α −γtanh z−z 0 δ ρ 0 u α ·∇u α +µ 2 ρ 0 ∇ u α ·u i +µ 2 ρ 0 ∇ u α · 1 2 z 2 α −z 2 ∇S+(z α −z)∇T +µ 2 ρ 0 ∇ u α ·Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) +µ 2 ρ 0 ξ+O(µ 4 ,γµ 2 ) (F.4) ρw(u) z =µ 2 ρ 0 z 2 S∇S+zT∇S+zS∇T +T∇T +w 0 ω 1 +O(µ 4 ,γµ 2 ) (F.5) 180 ∇p = ∇{ρ 0 (ζ−z)} − γ∇ ρ 0 δ ln cosh ζ−z 0 δ −ln cosh z−z 0 δ + µ 2 ρ 0 1 2 z 2 ∇S t − 1 2 ∇ ζ 2 S t +z∇T t −∇(ζT t )+ 1 2 z 2 ∇(u α ·∇S) − 1 2 ∇ ζ 2 u α ·∇S +z∇(u α ·∇T)−∇(ζu α ·∇T)− 1 2 z 2 ∇S 2 + 1 2 ∇ ζ 2 S 2 −z∇(TS)+∇(ζTS) + µ 2 (∇ρ 0 ) 1 2 z 2 −ζ 2 S t +(z−ζ)T t + 1 2 z 2 −ζ 2 u α ·∇S + (z−ζ)u α ·∇T − 1 2 z 2 −ζ 2 S 2 −(z−ζ)TS + O(µ 4 ,αµ 3 ,βµ 3 ,γµ 2 ) (F.6) αµ ∇· ρν h t ∇u =αµ ∇· ρ 0 ν h t ∇u α +O(αµ 3 ,αγµ ) (F.7) β µ (ρν v t u z ) z =−βµρ 0 ν v t ∇S−βµ τ b ζ +h +βµρ 0 ν v t u i zz +O(βµ 3 ,βγµ ) (F.8) In Equation F.4,ξ is defined as, ξ = − u α × ∇× 1 2 z 2 α −z 2 ∇S+(z α −z)∇T − 1 2 z 2 α −z 2 ∇S+(z α −z)∇T ×(∇×u α ) − u α × ∇×u i −u i ×(∇×u α ) − u α ×(∇×Ψ) 1 2 z 2 α −z 2 +ζ(z−z α ) − Ψ×(∇×u α ) 1 2 z 2 α −z 2 +ζ(z−z α ) (F.9) 181 Equation F.3 to F.8 are substituted into Equation 4.4 to produce the horizontal mo- mentum equation written in terms ofu α . Thus, (u α ) t + u α ·∇u α +∇ζ + ∇ρ 0 ρ 0 (ζ−z) + µ 2 − 1 2 ∇ ζ 2 S t −∇(ζT t )+ 1 2 ζ 2 ∇S t +z α ∇T t − 1 2 ∇ ζ 2 u α ·∇S −∇(ζu α ·∇T)+ 1 2 ∇ ζ 2 S 2 + 1 2 ∇ z 2 α u α ·∇S +∇(z α u·∇T)+T∇T +∇(ζTS) + µ 2 Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) t + µ 2 ∇ u α ·Ψ 1 2 z 2 α −z 2 +ζ(z−z α ) − µ 2 Ψ(ζ−z)(zS+T)+µ 2 ξ+µ 2 u i t +∇u α ·u i + µ 2 ∇ρ 0 ρ 0 1 2 z 2 −ζ 2 S t +(z−ζ)T t + 1 2 z 2 −ζ 2 u α ·∇S + (z−ζ)u α ·∇T − 1 2 z 2 −ζ 2 S 2 −(z−ζ)TS − γ tanh z−z 0 δ (u α ) t −tanh z−z 0 δ u α ·∇u α − 1 ρ 0 ∇ ρ 0 δ ln cosh ζ−z 0 δ −ln cosh z−z 0 δ = αµ 1 ρ 0 ∇· ρ 0 ν h t ∇u α −βµν v t ∇S−βµ 1 ρ 0 τ b ζ +h +βµν v t u i zz (F.10) Now, the remaining procedure is to eliminate z dependency in the above equation; depth-averaging is employed (i.e. Chen(2006)) over the entire equation. This final equa- tion is provided in the main body of the manuscript. 182 Appendix G Numerical Formulation The derived Equations 4.24 and 4.29 will be discretized to find numerical solutions. In the present work, a conservative-form finite volume method is adopted for spatial derivatives while 3rd order Adams-Bashforth predictor and 4th order Adams-Moulton corrector scheme is used for time integration. Prior to discretization of the governing system, Equation 4.24 and 4.29 are converted to the conservative form before applying the finite volume method. In this section, all dimensionsarerecoveredwithprimes( ′ )omittedforconvenience. Utilizingafixedbottom assumption(h t = 0), the conservative form of the continuity and momentum equations can be obtained as, H t +(Hu α ) x +(Hv α ) y +(N D +N B +N I )=0 (G.1) (Hu α ) t + Hu 2 α + 1 2 gH 2 x +(Hu α v α ) y −gHh x + 1 2 (ρ 0 ) x ρ 0 gH 2 + HM x +u α (N D +N B +N I )=0 (G.2) 183 (Hv α ) t + (Hu α v α ) x + Hv 2 α + 1 2 gH 2 y −gHh y + 1 2 (ρ 0 ) y ρ 0 gH 2 + HM y +v α (N D +N B +N I )=0 (G.3) where H =ζ +h and terms of O(µ 2 ,γ,αµ,βµ ) are given by (M x ,M y ) = R D +R B +R I +R h P +R v P +ξ − 1 ρ 0 ∇· ρ 0 ν h t ∇u α +ν v t ∇S+ τ b Hρ 0 − ν v t u i z | z=ζ −u i z | z=−h (G.4) G.1 Time Integration Timederivativetermsintheaboveequationsaresolvedbyathird-orderAdams-Bashforth predictor and a fourth-order Adams-Moulton corrector scheme(Wei et al.(1995), Lynett and Liu(2002a)) to minimize truncation error to order of O(Δt 3 )(Liu and Wang(2012)). Through an iterative predictor-corrector time-marching scheme, the solution at the next time step, (n+1) can be found. The explicit predictor step is given by, ζ n+1 =ζ n + Δt 12 23E n −16E n−1 +5E n−2 (G.5) P n+1 = P n + Δt 12 23F n −16F n−1 +5F n−2 +2F n 3 −3F n−1 3 +F n−2 3 +F p 4 (G.6) Q n+1 = Q n + Δt 12 23G n −16G n−1 +5G n−2 +2G n 3 −3G n−1 3 +G n−2 3 +G p 4 (G.7) 184 and the implicit corrector step is written as, ζ n+1 =ζ n + Δt 24 9E n+1 +19E n −5E n−1 +E n−2 (G.8) P n+1 = P n + Δt 24 9F n+1 +19F n −5F n−1 +F n−2 +F n+1 3 −F n 3 +F c 4 (G.9) Q n+1 = Q n + Δt 24 9G n+1 +19G n −5G n−1 +G n−2 +G n+1 3 −G n 3 +G c 4 (G.10) where the superscript n denotes time step and P, Q are defined numerically as (Kim et al.(2009)) P = (u α ) i−1,j H i,j z 2 α −ζ 2 2Δx 2 + (z α −ζ)h i−1,j Δx 2 + ζ x ζ 2Δx + ζ x h i−1,j 2Δx + (u α ) i,j H i,j 1− z 2 α −ζ 2 Δx 2 − 2(z α −ζ)h i,j Δx 2 (G.11) + (u α ) i+1,j H i,j z 2 α −ζ 2 2Δx 2 + (z α −ζ)h i+1,j Δx 2 − ζ x ζ 2Δx − ζ x h i−1,j 2Δx Q = (v α ) i,j−1 H i,j z 2 α −ζ 2 2Δy 2 + (z α −ζ)h i,j−1 Δy 2 + ζ y ζ 2Δy + ζ y h i,j−1 2Δy + (v α ) i,j H i,j 1− z 2 α −ζ 2 Δy 2 − 2(z α −ζ)h i,j Δy 2 (G.12) + (v α ) i,j+1 H i,j z 2 α −ζ 2 2Δy 2 + (z α −ζ)h i,j+1 Δy 2 − ζ y ζ 2Δy − ζ y h i,j−1 2Δy Subscript (i, j) in P and Q identifies cell location. All other terms included in Equa- tion G.5 to G.10 are given below: 185 E =E 1 +E 2 (G.13) F =F 1 +F 2 +u α E 2 (G.14) G=G 1 +G 1 +v α E 2 (G.15) E 1 , F 1 , and G 1 are rewritten by E 1 =−{Hu α } x −{Hv α } y (G.16) F 1 =− Hu 2 α + 1 2 gH 2 x −{Hu α v α } y +gHh x − 1 2 (ρ 0 ) x ρ 0 gH 2 (G.17) G 1 =−{Hu α v α } x − Hv 2 α + 1 2 gH 2 y +gHh y − 1 2 (ρ 0 ) y ρ 0 gH 2 (G.18) and E 2 , F 2 , G 2 , F 3 , G 3 are expressed as E 2 = H ζ 2 −ζh+h 2 6 − 1 2 z 2 α ∇S+ ζ−h 2 −z α ∇T x + H ζ 2 −ζh+h 2 6 − 1 2 z 2 α ∇S+ ζ−h 2 −z α ∇T y − " Hψ x ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# x − " Hψ y ( z 2 α 2 −z α ζ + 2ζ 2 −2ζh−h 2 6 )# y − h Hu i i x − h Hv i i y (G.19) 186 (F 2 ,G 2 ) = H 1 2 ∇ ζ 2 u α ·∇S +∇(ζu α ·∇T)− 1 2 ∇ ζ 2 S 2 − 1 2 ∇ z 2 α u α ·∇S −∇(z α u α ·∇T)−∇(ζTS) − (T∇T)−∇{E(ζS+T)}]−Hξ − E(ζS+T)∇ζ− 1 2 ζ 2 −z 2 α E∇S−(ζ−z α )E∇T + H " ζ 2 −ζh+h 2 6 ∇(u α ·Ψ)− (ζ−h) 2 ∇{u α ·(Ψζ)} + Ψ ( ζ 2 +ζh−2h 2 S 6 + HT 2 ) − ∇ u α · Ψ z 2 α 2 −ζz α −H∇ u α ·u i − ∇ρ 0 ρ 0 H(−2ζ +h) 6 (u α ·∇S−S 2 )− H 2 (u α ·∇T −ST) − δ(u α ·∇u α ) ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ + H ρ 0 ∇ ρ 0 δlncosh ζ−z 0 δ − 1 ρ 0 Z ζ −h ∇ ρ 0 δlncosh z−z 0 δ dz + H 1 ρ 0 ∇· ρ 0 ν h t ∇u α −ν v t ∇S− τ b Hρ 0 + ν v t u i z | z=ζ −u i z | z=−h (G.20) F 3 = 1 2 H ζ 2 −z 2 α (v α ) xy −H(z α −ζ)(hv α ) xy + Hζ x n ζ(v α ) y +(hv α ) y o (G.21) 187 G 3 = 1 2 H ζ 2 −z 2 α (u α ) xy −H(z α −ζ)(hu α ) xy + Hζ y {ζ(u α ) x +(hu α ) x } (G.22) F p 4 , G p 4 , F c 4 and G c 4 are rewritten by F p 4 = H n ζ 2 −ζh+h 2 −3z 2 α n 6 Σ p (ψ x )− H n (ζ−h−2z α ) n 2 Σ p (ψ x ζ) + H n Σ p (u i )− H n (h−2ζ) n 6 (ρ 0 ) x ρ 0 Σ p (S)+ H n 2 (ρ 0 ) x ρ 0 Σ p (T) − δ ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ n Σ p (u α ) (G.23) G p 4 = H n ζ 2 −ζh+h 2 −3z 2 α n 6 Σ p (ψ y )− H n (ζ−h−2z α ) n 2 Σ p (ψ y ζ) + H n Σ p v i − H n (h−2ζ) n 6 (ρ 0 ) y ρ 0 Σ p (S)+ H n 2 (ρ 0 ) y ρ 0 Σ p (T) − δ ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ n Σ p (v α ) (G.24) F c 4 = H n+1 ζ 2 −ζh+h 2 −3z 2 α n+1 6 Σ c (ψ x )− H n+1 (ζ−h−2z α ) n+1 2 Σ c (ψ x ζ) + H n+1 Σ c u i − H n+1 (h−2ζ) n+1 6 (ρ 0 ) x ρ 0 Σ c (S)+ H n+1 2 (ρ 0 ) x ρ 0 Σ c (T) − δ ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ n+1 Σ c (u α ) (G.25) 188 G c 4 = H n+1 ζ 2 −ζh+h 2 −3z 2 α n+1 6 Σ c (ψ y )− H n+1 (ζ−h−2z α ) n+1 2 Σ c (ψ y ζ) + H n+1 Σ c v i − H n+1 (h−2ζ) n+1 6 (ρ 0 ) y ρ 0 Σ c (S)+ H n+1 2 (ρ 0 ) y ρ 0 Σ c (T) − δ ln cosh −h−z 0 δ −ln cosh ζ−z 0 δ n+1 Σ c (v α ) (G.26) where Σ p (φ)=2φ n −3φ n−1 +φ n−2 and Σ c (φ)=φ n+1 −φ n . G.2 Spatial Discretization : Finite Volume Method Recently, finite volume schemes coupled with Riemann solvers have been successfully ap- pliedtoshallowwater(Erduranetal.(2005))andBoussinesq-type(TonelliandPetti(2009), Kim et al.(2009), Shi et al.(2012)) equations, and have shown relatively robust per- formance. For the shallow water terms embedded in Equation G.1, G.2 and G.3, a 4th order compact MUSCL-TVD scheme has been applied and combined with the HLL Riemann solver(Kim et al.(2009)). The fractional volume flux determination procedure used in this study is described here. The interface value of any conservative quantity (q =H,Hu α ,Hv α ) is constructed by, q L i+1/2 =q i + 1 6 n Δ ∗ q i−1/2 +2Δ ∗ e q i+1/2 o (G.27) q R i+1/2 =q i+1 − 1 6 n 2Δ ∗ q i+1/2 +Δ ∗ e q i+3/2 o (G.28) 189 where superscriptL andR denotes the values at the left-hand and the right-hand side of the interface(i+1/2), respectively. Δ ∗ q and Δ ∗ e q in the above can be calculated from Δ ∗ q i−1/2 =minmod Δ ∗ q i−1/2 ,bΔ ∗ q i+1/2 Δ ∗ e q i+1/2 =minmod Δ ∗ q i+1/2 ,bΔ ∗ q i−1/2 Δ ∗ q i+1/2 =minmod Δ ∗ q i+1/2 ,bΔ ∗ q i+3/2 Δ ∗ e q i+3/2 =minmod Δ ∗ q i+3/2 ,bΔ ∗ q i+1/2                                      (G.29) in which Δ ∗ q i+1/2 =Δq i+1/2 − 1 6 Δ 3 q i+1/2 (G.30) Δ 3 q i+1/2 =Δq i−1/2 −2Δq i+1/2 +Δq i+3/2 (G.31) Δq i−1/2 =minmod Δq i−1/2 ,b 1 Δq i+1/2 ,b 1 Δq i+3/2 Δq i+1/2 =minmod Δq i+1/2 ,b 1 Δq i+3/2 ,b 1 Δq i−1/2 Δq i+3/2 =minmod Δq i+3/2 ,b 1 Δq i−1/2 ,b 1 Δq i+1/2                        (G.32) The Minmod limiter function used in Equation G.29 and G.32 is defined by minmod(x,y)=sign(x)max[0,min{|x|,y sign(x)}] minmod(x,y,z)=sign(x)max[0,min{|x|,y sign(x),z sign(x)}]            (G.33) 190 The coefficients b and b 1 in Equation G.29 and G.32 are set to 2 and 4, respectively. It is noted that ignoring the 3rd order component(Δ 3 q) in Equation G.30 will degrade the scheme from 4th to 3rd order. Constructed interface values(q L ,q R ) poses local Riemann problem, so are input for computing numerical fluxes through a HLL approximate Riemann solver, given by (Toro, 2002) Θ q L ,q R =                        Θ q L if s L ≥0 Θ ∗ q L ,q R if s L <0<s R Θ q R if s R ≤0 (G.34) where Θ ∗ q L ,q R = s R Θ q L −s L Θ q R +s R s L q R −q L s R −s L (G.35) The wave speeds used here are given by s L =u L − p gH L ϕ L , s R =u R + p gH R ϕ R (G.36) in which flux ϕ (L,R) is given by ϕ (L,R) =              r 1 2 (H∗+H (L,R) )H ∗ H (L,R) 2 , H ∗ >H (L,R) (shock) 1 , H ∗ ≤H (L,R) (rarefaction) (G.37) H ∗ = 1 g 1 2 p gH L + p gH R + 1 4 u L −u R 2 (G.38) 191 The remaining terms, including higher order spatial derivatives, are differenced by the cell averaged finite volume method proposed by Lacor et al.(2003). 192 Appendix H Second-order Sub- and Super-harmonic Solution in Two-layer Fluids a ± exact = P 5 j=1 Π ± j Λ ± j Det ± (H.1) where Π ± 1 = −K ± 25 K ± 33 K ± 42 K ± 51 +K ± 23 K ± 35 K ± 42 K ± 51 +K ± 25 K ± 33 K ± 41 K ± 52 − K ± 23 K ± 35 K ± 41 K ± 52 (H.2) Π ± 2 =K ± 15 K ± 33 K ± 42 K ± 51 −K ± 15 K ± 33 K ± 41 K ± 52 (H.3) Π ± 3 =−K ± 15 K ± 23 K ± 42 K ± 51 +K ± 15 K ± 23 K ± 41 K ± 52 (H.4) Π ± 4 =K ± 25 K ± 33 K ± 51 −K ± 23 K ± 35 K ± 51 −K ± 15 K ± 23 K ± 31 K ± 52 (H.5) 193 Π ± 5 =−K ± 25 K ± 33 K ± 41 +K ± 23 K ± 35 K ± 41 +K ± 15 K ± 23 K ± 31 K ± 42 (H.6) Λ ± 1 =− 1 2 k e A e a i0 +k i A i a e0 (H.7) Λ ± 2 =− 1 2 k e cosh(k e h l )C e a i0 +k i cosh k i h l C i a e0 (H.8) Λ ± 3 = 1 2 ρ u ρ l k e σ e B e a i0 +k i σ i B i a e0 +k e k i (A e A i ±B e B i ) − 1 2 k e σ e sinh(k e h l )C e a i0 +k i σ i sinh k i h l C i a e0 −k e k i C e C i cosh k ∓ h l (H.9) Λ ± 4 =− 1 2 k e a i {A e cosh(k e h u )+B e sinh(k e h u )} +k i a e A i cosh k i h u +B i sinh k i h u (H.10) Λ ± 5 = 1 2 k e σ e a i {A e sinh(k e h u )+B e cosh(k e h u )} +k i σ i a e A i sinh k i h u +B i cosh k i h u −k e k i a e (A e A i ∓B e B i )cosh k ∓ h u +(B e A i ∓B i A e )sinh k ∓ h u (H.11) Det ± =K ± 25 K ± 33 K ± 44 K ± 51 −K ± 23 K ± 35 K ± 44 K ± 51 −K ± 15 K ± 23 K ± 31 K ± 44 K ± 52 −K ± 25 K ± 33 K ± 41 K ± 54 +K ± 23 K ± 35 K ± 41 K ± 54 +K ± 15 K ± 23 K ± 31 K ± 42 K ± 54 (H.12) 194 Variables A i ,B i ,C i ,a i and K ± ij are, A i = a i0 σ i " g 1− ρ l ρ u + ρ l σ i 2 ρ u k i tanh(k i h l ) # (H.13) B i = σ i a i0 k i (H.14) C i = σ i a i0 k i sinh(k i h l ) (H.15) a i = σ i 2 a i0 (σ i ) 2 cosh(k i h u )−gk i sinh(k i h u ) (H.16) and K ± 15 =σ ± /k ± , K ± 23 =sinh k ± h l , K ± 25 =σ ± /k ± , K ± 31 =−σ ± ρ l /ρ u , K ± 33 =σ ± cosh k ± h l , K ± 35 =g(ρ l /ρ u −1), K ± 41 =sinh k ± h u , K ± 42 =cosh k ± h u , K ± 44 =σ ± /k ± K ± 51 =−σ ± cosh k ± h u , K ± 52 =−σ ± sinh k ± h u , K ± 54 =g (H.17) Similarly, A e ,B e ,C e ,a e can be expressed. 195

Asset Metadata

Creator Son, Sangyoung (author)

Core Title Wave induced hydrodynamic complexity and transport in the nearshore

Contributor Electronically uploaded by the author (provenance)

School Andrew and Erna Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 07/25/2014

Defense Date 06/11/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag Boussinesq equations,coherent structures,eddies,internal waves,long waves,nearshore,OAI-PMH Harvest,shallow water,stratified flows,Transport,tsunami,tsunami whirlpool,turbulence,two-way coupling,wave-current interactions

Language English

Advisor Lynett, Patrick J. (committee chair), Lee, Jiin-Jen (committee member), Redekopp, Larry G. (committee member)

Creator Email sangyous@usc.edu,sonsiest@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-66462

Unique identifier UC11289369

Identifier usctheses-c3-66462 (legacy record id)

Legacy Identifier etd-SonSangyou-1005.pdf

Dmrecord 66462

Document Type Dissertation

Rights Son, Sangyoung

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Abstract (if available)

Abstract In the coastal area, defined as the region between the shoreline and some offshore limit where the depth can no longer influence the waves, complex behavior of waves is anticipated due to various physical effects such as turbulence, wave-structure interaction, wave-current interaction, wave breaking and fluid-density variations. For modeling of nearshore hydrodynamics, many numerical models have been developed so far, but many of such effects are not yet considered appropriately. ❧ In this dissertation, depth-integrated numerical models used in long wave simulation are developed for better understanding of complicated hydrodynamics at the nearshore. First, a non-dispersive shallow water equation model and dispersive Boussinesq model are two-way coupled. The fundamental purpose of the coupling effort is to develop the capability to seamlessly model long wave evolution from deep to shallow water with fine scale resolution, without the loss of locally important physics. Second, a set of depthintegrated equations describing combined wave-current flows are derived mathematically and discretized numerically. To account for the effect of turbulent interaction between waves and underlying currents with arbitrary profile, new additional stresses are introduced, which represent radiation stress of waves over the ambient current field. Finally, a numerical model for gravity waves propagating over variable density fluids is developed by allowing horizontal and vertical variation of fluid density. Throughout the derivation, density change effects appear as correction terms while the internal wave effects on the free surface waves in a two-layer system are accounted for through direct inclusion of the internal wave velocity component. For each of the studied topics, numerical tests are performed to support accuracy and applicability. Consequently, we have developed a comprehensive tool for numerical simulation of complex nearshore hydrodynamics.