University of Southern California Dissertations and Theses

Tsunami-induced turbulent coherent structures

(USC Thesis Other)

Tsunami-induced turbulent coherent structures

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content TSUNAMI-INDUCED TURBULENT COHERENT STRUCTURES by Nikos Kalligeris A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL AND ENVIRONMENTAL ENGINEERING) May 2017 Copyright 2017 Nikos Kalligeris Acknowledgements I would like to acknowledge the funding provided through the Viterbi and Myronis fellowships offered by the University of Southern California. The work presented in this thesis was funded by the National Science Foundation NEES Research program, with award number 1135026. The experimental work in Oregon State University’s Hinsdale Wave Research Laboratory was a collaborative effort between a group of people. Aykut Ayca was a fellow investigator in the experiments and his big contribution to this work is greatly appreciated. So is Adam Ryan’s, who besides being great company, he helped a lot setting up the lab instrumentation and also provided input during the preliminary data-processing stage. My thanks also go to Hinsdale lab personnel Tim Maddux, Alicia Lyman-Holt, Jeff Gent and Cooper for all their help and the happy times. It was a delightful summer in Corvalis. During my graduate studies, I had the privilege and honor to collaborate with many great scientists and mentors. Professor Costas Synolakis introduced me to tsunamis and believed in me from the start of our collaboration. The rest of this great journey is a result of his support, encouragement and guidance. My field-work mentor cannot be other than Professor Hermann Fritz, who welcomed me to my first tsunami field survey in Java, Indonesia. It was the start of many other field surveys we did together. I thank him for training me in the field and trusting me in so many occasions to be his field collaborator. I also like to thank Professor Emile Okal for being such a great mentor. It has been a privilege and a pleasure working with him. Last, but not least, I want to thank Dr. Jose Borrero for doing some great work together, and I look forward to more collaborations in the future. This part is dedicated to all my teachers who enlightened me by passing on their knowledge, and gave me the incentive to explore my own ideas. My biggest gratitudes go to my PhD advisors, Professors Lynett and Synolakis, for their support and encouragement throughout my time at USC. Professor Lynett gave me the opportunity to work on this project and encouraged me from the start to pursuit my ideas. He has been such a great mentor, always offering constructive comments and ideas which shaped the outcome of this work. He trusted me with the task at hand, and provided the support for the duration of my studies. Costas, if I may, it’s been a great pleasure working with you. You have been more than a mentor and a teacher for me, and I will always be grateful for your guidance in education and beyond. A big thank you to all my friends, both in Greece and LA whos companionship has shaped me and kept me busy outside of work. Guys we made it! I want to thank Aykut Ayca, Vassilios Skanavis and Antonios Papantoniou for being such great roommates in LA. I also want to thank my lab mates, the fellow worshipers of coastal engineering, with whom we shared good and tough times. Discussions and collaborations with them led to some of the ideas presented in this work. A special thanks goes to Sasan Tavakkol for all the great discussions and suggestions. I’m leaving my parents and my brother to the end, who are my life’s example to follow. I can’t thank you enough for your love and support. I’ve made it this far because of you. iii Contents The table of contents is hyperlinked List of Tables vi List of Figures vii Abstract xv 1 Introduction 1 1.1 Objectives of Study and Proposed Approach . . . . . . . . . . . . . . . . . . . . . 9 2 Experimental Setup 12 2.1 Boundary conditions, water level and experiment scaling . . . . . . . . . . . . . . 14 2.2 Cameras and light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 ADV’s and wave gauges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3 Post Processing 22 3.1 Image pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1.1 Camera calibration - intrinsic parametrization . . . . . . . . . . . . . . . . 23 3.2 Particle identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 2D Particle Tracking Velocimetry (PTV) . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 3D Particle Tracking Velocimetry (PTV) . . . . . . . . . . . . . . . . . . . . . . . 30 3.4.1 V orono¨ ı particle tracking method . . . . . . . . . . . . . . . . . . . . . . . 32 3.4.2 Ray intersection and matching algorithm . . . . . . . . . . . . . . . . . . 35 3.5 Extrinsic parametrization - coordinate transformation . . . . . . . . . . . . . . . . 37 3.5.1 2D rectification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5.2 3D rectification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5.3 3D ray parametric equation . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.6 Surface tracer’s horizontal position correction . . . . . . . . . . . . . . . . . . . . 45 3.7 2D surface velocity extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.8 3D surface velocity extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.9 Conversion to polar coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.9.1 Identification of the TCS center . . . . . . . . . . . . . . . . . . . . . . . 51 3.9.2 TCS-center paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.9.3 TCS-center propagation velocity . . . . . . . . . . . . . . . . . . . . . . . 54 3.10 Lagrangian to Eulerian - gradients, interpolation, ensemble averaging, and azimuthal averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.10.1 Interpolation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.10.2 Space-averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.10.3 Spatial derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.11 V orticity measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.12 ADV data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.13 Free surface elevation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 iv 4 TCS generation and circulation growth 79 4.1 Theoretical background on circulation growth . . . . . . . . . . . . . . . . . . . . 79 4.2 Circulation growth theoretical framework for the experimental layout . . . . . . . 82 4.2.1 Local coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 V ortex identification and quantification . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.1 V ortex identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3.2 Circulation quantification . . . . . . . . . . . . . . . . . . . . . . . . . . 86 4.3.3 Upwelling/downwelling quantification . . . . . . . . . . . . . . . . . . . . 88 4.4 Inshore TCS generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.4.1 Circulation growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.5 Main (offshore) TCS generation . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.5.1 Circulation growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.5.2 Other characteristics of the starting-jet vortex . . . . . . . . . . . . . . . . 104 5 TCS characteristics, evolution and flow transition to Q2D 108 5.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.2 Monopolar vortex theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.3 2D flow field evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.4 Azimuthal-averaged properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.5 Decay of mean vortex properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.5.1 V ortex decay model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.5.2 TCS growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.6 Flow transition to Q2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.6.1 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 5.6.2 Experimental observations and secondary flow quantification . . . . . . . . 137 6 Field measurements of vortical motions in Ventura harbor during the 2015 Chilean tsunami 143 6.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.1.1 Optical particle tracking. . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.1.2 GPS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.1.3 Numerical Modelling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 7 Conclusions 154 Reference List 158 v List of Tables 2.1 Experimental trial short description. . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Froude scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1 Mean and max coordinate transformation error for each camera. . . . . . . . . . . 41 3.2 Accuracy of 3D coordinate transformation. . . . . . . . . . . . . . . . . . . . . . 44 3.3 Mean SNR and correlation values for all ADV beams of trial 3. . . . . . . . . . . . 73 vi List of Figures 1.1 The path freighter Maersk Mandraki followed after it was caught in the turbulent eddies that broke its moorings during the the 2004 Indonesian tsunami. The path was approximately hand-drawn by an eyewitness. Source: Okal et al. (2006a) . . . 2 1.2 Near-shore rotational flow structures observed during the 2011 Japan tsunami in Oarai harbor (left) and Iwaki city (right). Source: Reuters. . . . . . . . . . . . . . 4 1.3 Left: Tidal flow leaving a sharp channel may induce the formation of a dipole by advection of vorticity created at the boundaries (Source: Wells and van Heijst (2003)). Right: Nozzle type piston-cylinder configuration to generate a vortex ring by advection of vorticity created at the cylinder boundary layer (Source: Krueger (2005)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 The experiment in three phases: 1) The wavemaker forward stroke creates a coher- ent structure on the inshore basin side; 2) the wave gets reflected on the back wall of the basin and the wavemaker retracts to create a strong current through the chan- nel that generates the offshore TCS; 3) The TCS detaches from the breakwater and gets advected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1 The experimental set-up. Dots denote acoustic Doppler velocimetry (ADV) instru- ment locations, squares denote the location of the overhead HD cameras used for the 2D PTV and dye experiments, blue dashed box denotes the study area of the basin and the red dashed box denotes the area shown in Figure 2.4b. The local coordinate system y-axis starts from the center-line and is positive upwards. The local coordinate system x-axis starts from the wavemaker mean position and is positive towards the basin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Wavemaker boundary conditions. (a) The piston displacement time history (0 displacement corresponds to x= 0 in basin coordinates), starting from the fully retreated position and reaching maximum stroke; (b) the recorded surface eleva- tion near the face of the wavemaker. . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Camera fields of view (FOVs) used for the 2D PTV experiments. . . . . . . . . . . 17 2.4 3D PTV experimental setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Basin light adjustments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 ADV and wave gauge used for the experiments. . . . . . . . . . . . . . . . . . . . 20 2.7 Positions of gauges recording water surface elevation. The inset shows the position of the wave gauges during the last two trials (red dots). . . . . . . . . . . . . . . . 21 3.1 Image pre-processing. (a) The original RGB still frame, as extracted from the video; (b) lens distortion is removed; (c) converted to grayscale and applied mask; (d) subtracted from mean video grayscale intensity and enhanced surfaced tracer intensity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Camera calibration images used to obtain camera intrinsic parameters. Example shown for one zoom level. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Particle center detection using a local peak above threshold algorithm (red crosses) and by applying the brightness weighted centroid offset of Equation 3.3 (green crosses) using w= 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 vii 3.4 Inter-frame displacement fluctuation using both the filtered and unfiltered images to detect the center of a random surface tracer (upper two figures): blue curves show the raw inter-frame displacement and red are smoothed using a low-pass Butterworth filter. The bottom figure shows the absolute difference in position (x- and y-coordinates) of the detected particle center between the two input images for each frame number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.5 The feature detection algorithm identifies single particles (green circles) and rejects groups of particles (red dots). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.6 Selected surface tracer paths projected on the mean grayscale FOV intensity. . . . . 31 3.7 Symbols and definitions for V orono¨ ı polygons. . . . . . . . . . . . . . . . . . . . 33 3.8 The particle screening (left plots) and selection (right plots) processes of the V orono¨ ı PTV method. (a) candidate matches P j;2 () for particle P i;1 (O) are sought in an area defined by a radius L around P i;1 . (b) remaining candidate particles P j;2 () after screening. (c) the V orono¨ ı 1-stars of a candidate S j;2 and its pair S i;1 . (d) V orono¨ ı 1-star S j;2 is translated to match the central vertex of S i;1 and apply Equa- tion 3.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.9 Intersection of two rays associated with the left and right camera viewpoints. Source: Douxchamps et al. (2005) . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.10 FOV from camera B showing the GCPs layed out on the (dry) basin prior to the 2D PTV experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.11 Image in Figure 3.10 rectified using the calibration coefficients obtained from Equation 3.22. The image acquires the right coordinates for any given point on the plane z= 0:55m. GCP basin coordinates are marked with the red crosses and the transformed coordinates from Equation 3.21 are shown with blue circles. . . . . 40 3.12 Calibration table used to obtain camera extrinsic parameters for the 3D PTV exper- iments. The bottom row GCPs correspond to the intersection points of the table background grid, while the top row GCPs are laid on the 4 4’s. . . . . . . . . . . 43 3.13 Geometrical depiction of the flat surface simplification: the horizontal coordinates of the surface tracer, that is elevated from the still water level, are erroneously calculated. Knowledge of its vertical displacement can reduce the error. . . . . . . 46 3.14 Left: horizontal distances between the overhead cameras and the water surface, covered by the respective FOVs. Right: error in estimation of horizontal coordi- nates of tracers due to the flat surface assumption (see Equation 3.31). The camera height was set to z c = 10m (a mean value of the elevation of the four cameras) and mean water depth was set to h 0 = 55cm. . . . . . . . . . . . . . . . . . . . . . . . 47 3.15 The mean error of tracer’s horizontal position in the overlapping camera FOVs without applying the water level correction. Legends correspond to the overlapping cameras FOV’s shown in Figure 2.3. . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.16 The mean error of tracer’s horizontal position in the overlapping camera FOVs after applying the water level correction. Legends correspond to the overlapping cameras FOVs shown in Figure 2.3. . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.17 The four camera FOVs merged together for a single time step. Tracer velocity vectors are overlaid. Each camera is assigned a different vector color. . . . . . . . 49 viii 3.18 Surface velocity vectors as obtained by the close-up stereo PTV analysis, 45 sec after the initiation of the experiment, showing how the flow separates at the tip of the breakwater and the the TCS is generated. Insets display a comparison between the surface elevation measured by two wave gauges (black curves) and the corre- sponding obtained by the stereo PTV analysis (red curves). The location of the wave gauges is shown with the red crosses. . . . . . . . . . . . . . . . . . . . . . 50 3.19 Tracer density around the TCS center at different times. Flow convergence/divergence either draws the tracers towards the TCS center or away from it. . . . . . . . . . . 52 3.20 Examples of detecting the edge of the tracer-conglomerate at the TCS center through image processing. At the early stages of TCS development, flow convergence keeps the tracer-conglomerate compact and in one piece. As time progresses, the flow becomes divergent, and the tracer-conglomerate expands and breaks up into smaller fractions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.21 Methodology for the extraction of the TCS center. In the early stages of the TCS development (up to point c1) - when the density of tracers around the TCS center is sufficient to interpolate the velocity field - the center is extracted from vorticity maps. In the next stage of TCS development, the edges of the blob are detected to compute the center of the TCS as the center of mass of the blob perimeter. The vortex path is subsequently filtered using different cutoff frequencies: 0.5s up to point f1, 0.2s up to point f2, and 0.015s until the end of the record. The filtered path is colored according to the corresponding vortex center speed. All images/subplots are to scale except the blob images showing the edge detection, which are scaled by a factor of two. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.22 TCS center path from one trial. As time progresses, the TCS moves slower and the frequency of the oscillations decreases. The path was filtered using multiple low-pass filters of decreasing cutoff frequencies. . . . . . . . . . . . . . . . . . . . 56 3.23 Paths of off-center tracers in an idealized shielded-Gaussian vortex flow withw 0 = 10 s 1 , R= 1 m, and vortex-center translation velocity u= 0:25 m=s. . . . . . . . 57 3.24 TCS center paths from all trials. The mean path is shown with the bold curve. . . . 57 3.25 TCS center velocity realizations (u;v) in x- and y-directions. The mean is shown with the bold curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.26 Stationary ensemble creation from the surface velocity vectors of the individual trials. This ensemble is used to extract the mean flow velocities during TCS gen- eration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.27 TCS-centered ensemble creation from the surface velocity vectors of the individual trials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.28 Ensemble surface velocity vectors and resulting mean surface velocity field near the breakwater tip at t= 27s, using W R = 0:35 m andDx=Dy= W R . Frame shown is after the arrival of the leading depression wave. . . . . . . . . . . . . . . . . . . 62 3.29 Number of tracers per interrogation window and standard error of the mean veloc- ities (Eqn. 3.47) for the stationary ensemble. Metrics evaluated in the time interval t=10-73 s using W R = 0:35 m andDx=Dy= W R . . . . . . . . . . . . . . . . . . . 64 3.30 Minimum distance of the TCS-center to a vertical boundary. The mean of all trials is shown with the solid line, whereas dashed lines correspond to the min/max envelope for the individual trials. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 ix 3.31 Statistics of the TCS-centered ensemble for the duration of the experiment. (a) The number of trials. (b) Mean tracer spacing in the domain extending to x;y2 [d i min ;d i min ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.32 Uncertainty of the mean velocity field for the TCS-centered ensemble. (a) Standard deviation of the u r and u q velocity components, evaluated at the TCS-center. (b) Standard error of the mean. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.33 Azimuthal-average sampling points for (a) variable spacing and constant dq, and (b) for constant spacing and variable dq; (c) interrogation windows in the form of annuli. The dashed red lines (b) correspond to interrogation windows of radius W R = dr and blue dashed lines (c) correspond to the annuli used as interrogation windows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.34 Azimuthal-averaging of the TCS azimuthal (u q ) and radial (u r ) velocity compo- nents - example shown for time t 333s. Black vectors and dots correspond to the scattered velocities and red curves to the azimuthal-averaged profiles. Error bars indicate the standard deviation of the sample velocities in each annulus of width W = 0:5m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.35 V orticity maps created using the least-squares and Richardson differential schemes at time t = 166s. Maximum vorticity (w max ) and circulation (G, evaluated inside the green circle) values are noted. The mean velocity field is evaluated on a grid withDx=Dy= 0:4m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.36 (a), (b) The two tracer configurations used to infer the vortex–center vorticity. The left configuration uses a 12.7 cm radius, whereas the right configuration uses a 4.9 cm radius; (c) Schematic of measuring the angular velocity of a square-tracer. . . . 70 3.37 The paths of the tracer configurations following the TCS center, overlaid on recti- fied images of the trial in which a square-tracer was deployed. . . . . . . . . . . . 72 3.38 Raw (a) and filtered (b) vorticity decay curves near the TCS center, as measured using the square- and cross-tracers. Colors match the paths shown in Figure 3.37. . 72 3.39 ADV raw (left column) and processed data (right column). . . . . . . . . . . . . . 73 3.40 Comparison between horizontal velocities extracted from ADV 1 (located at x= 4:52m, y= 0:00m), ADV 3 (located at x= 20:69m, y= 13:18m) and 2D PTV . . . . 74 3.41 Delaunay triangulation mesh (blue lines) to interpolate the wave gauge surface elevation data (red dots). The breakwater perimeter is set as a hard boundary. . . . 75 3.42 Interpolated surface elevation map in the study area 5 s after wavemaker start. The inset shows the time relative to the wameker motion. Wave gauge locations shown with black dots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.43 Comparison of surface elevation extracted from the surface tracers (3D PTV , red dashed curves) and from the wave gauges (black curves). . . . . . . . . . . . . . . 77 3.44 Comparison of surface elevation extracted from the surface tracers (3D PTV) and from the wave gauges (continued from Figure 3.43). . . . . . . . . . . . . . . . . . 78 3.45 Wave gauge numbering corresponding to Figures 3.43-3.44. . . . . . . . . . . . . 78 4.1 Schematic of (a) nozzle and (b) orifice type piston-cylinder vortex ring generators. Source: Krueger (2005) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Bounding contour C to determine the circulation influx to the offshore basin during 1 st flow reversal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 x 4.3 Channel-average velocity evaluated using the offshore arc of length 3:2m shown in Figure 4.4 as U = Q=3:2, where Q is the flow rate per unit depth: Q= R l 0 ru r dq. . . 83 4.4 Polar coordinate system (r;q) used to compute circulation influx through the exper- imental channel. Coordinate z is the vertical (out of the page). The blue polygons define the areas over which the total circulation is evaluated in each side of the basin. 85 4.5 Radius vortex definition through an example for an axisymmetric vortex centered at(0;0). (a) Azimuthal velocity profile. (b) Swirl-strength map and velocity vec- tors. (c) V orticity map. Yellow circles are plotted every 1m in the radial direction. The magenda circle in (b) shows thel ci = 0 contour, and the red circle in (c) shows the w = 0 contour. The radii corresponding to the two contours are plotted in (a) with vertical dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.6 Swirl strength (b) and vorticity (c) maps corresponding to an ensemble flow field (a) during TCS generation at the inshore side of the breakwater. Swirl strength and vorticity contours are plotted with a step of 0.1 and 0.3 s 1 respectively. The lowest contours correspond tol ci = 0 andjwj= 0:15 s 1 . . . . . . . . . . . . . . 88 4.7 Velocity vector maps at various times during the inshore TCS generation. . . . . . 90 4.8 V orticity (w) maps at various times during the inshore TCS generation. Green polygons correspond to the l ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.9 Swirl strength (l ci sign(w)) maps at various times during the inshore TCS gen- eration. Green polygons correspond to thel ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.10 (a) Flow profiles across the channel for different times during the generation of the inshore TCS. (b,c) The corresponding velocity components in the polar coordinate system. Profiles are color-matched for time, and the aspect ratio of the flow profiles in (a) is skewed for presentation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.11 Velocity gradient and vorticity profiles in the channel for different times during the generation of the inshore TCS. Profiles are color-matched according to the times shown in Figure 4.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.12 Circulation growth during the inshore TCS generation. (a) Total circulation com- pared against the vorticity influx parameters of Eqn. 4.13. (b) Total circulation compared with the TCS-enclosed circulation, as defined in Eqn. 4.15. Arrow denotes the pinch-off time of the TCS. . . . . . . . . . . . . . . . . . . . . . . . . 95 4.13 Flow profiles across the channel for different times (a) during the generation of the offshore TCS. The corresponding velocity components in the curvilinear coor- dinate system are plotted in (b) and (c) with continuous lines. Profiles are color- matched for time, and the aspect ratio of the flow profiles in (a) is skewed for presentation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.14 Velocity gradient and vorticity profiles in the channel for different times during the generation of the offshore TCS. Profiles are color-matched according to the times shown in Figure 4.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 xi 4.15 Circulation growth during the offshore TCS generation. (a) Total circulation com- pared against the vorticity influx parameters of Eqn. 4.13. (b) Total circulation compared with the TCS-enclosed circulation, as defined in Eqn. 4.15. Arrows denote the pinch-off time of the TCS and the time of flow reversal. . . . . . . . . . 99 4.16 Velocity vector maps at various times during the offshore TCS generation. . . . . . 100 4.17 V orticity (w) maps at various times during the offshore TCS generation. Green polygons correspond to the l ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.18 Swirl strength (l ci sign(w)) maps at various times during the offshore TCS gen- eration. Green polygons correspond to thel ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.19 Upwelling/downwelling maps at various times during the offshore TCS genera- tion. Units shown in the colorbar correspond to 10 3 cm 3 =s. Upwelling has negative value. Red circles indicate the TCS center. . . . . . . . . . . . . . . . . . . . . . . 103 4.20 Images from a dye visualization experiment captured during the offshore TCS gen- eration. Fluorescent green dye is released from the breakwater tip (inside the sepa- ration zone) and fluorescent red dye is released just inshore of the tip. The images are taken from an oblique angle and have not been undistorted/rectified. . . . . . . 104 4.21 Figure 4.20 continued. Images taken from another camera and angle. . . . . . . . . 105 4.22 Metrics describing the offshore starting-jet vortex dynamics. (a) V ortex-equivalent- radius growth as a function of formation timet. (b) V ortex horizontal displacement Vs equivalent piston stroke. (c) Horizontal displacement. (d) Vertical displacement. 107 4.23 Trajectories of the centers of the starting-jet vortices for the isnhore and offshore flow phases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.1 TCS-centered vorticity maps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.2 TCS-centered flow divergence (¶w=¶z=(¶u=¶x+¶v=¶y)) maps. Downwelling is positive. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.3 TCS-centered vorticity maps. Continued from Figure 5.1. . . . . . . . . . . . . . . 117 5.4 TCS-centered flow divergence maps. Continued from Figure 5.2. . . . . . . . . . . 118 5.5 TCS azimuthal-averaged u q and w profiles. (a) Azimuthal velocity (u q ) normal- ized with u q;max and its corresponding radius R vmax . (b) V orticity normalized with w max . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.6 TCS azimuthal-averaged u q velocity profiles for times t = 100 500s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 122 5.7 TCS azimuthal-averaged u q velocity profiles for times t = 550 950s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 122 5.8 TCS azimuthal-averaged u r velocity profiles for times t = 100 500s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 123 5.9 TCS azimuthal-averaged u r velocity profiles for times t = 550 950s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 123 5.10 TCS azimuthal-averaged vorticity profiles for times t = 100 500s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 124 xii 5.11 TCS azimuthal-averaged vorticity profiles for times t = 550 950s. Error bars represent the standard deviation of the sample velocities. . . . . . . . . . . . . . . 124 5.12 a-profile fit for times t= 100 500s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . . . . 125 5.13 a-profile fit for times t= 550 950s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . . . . 125 5.14 a-profile fit for times t = 1000 1400s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . 126 5.15 a-profile fit for times t = 1450 1850s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . 126 5.16 Whittaker-profile fit for times t= 100500s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . 127 5.17 Whittaker-profile fit for times t= 550950s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. . . . . . . . . . 127 5.18 Decay of TCS mean properties. (a) Maximum azimuthal (black line) and minimum radial (blue circles) velocity decay. (b) maximum vorticity decay. The 1 st order decay prediction is plotted with the dashed red lines. . . . . . . . . . . . . . . . . 130 5.19 (a) TCS size growth with time. Dashed lines correspond to viscous diffusion rates as functions of p t, red crosses and green squares to the TCS-core growth (from the azimuthal-averaged u q profiles and the u q profile fits), black circles to TCS bulk radius growth, and the cyan curve to the mean-minimum distance to the closest vertical boundary. (b) Local depth-based Reynolds number decay (Eqn. 5.38). . . . 132 5.20 TCS center paths from all trials. The mean path is shown with the bold curve. The dashed line corresponds to the maximum circle fitting the polygon defined by the solid boundaries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.21 Quantification of the secondary flow magnitude. (a) Raw (black) and smoothed (dashed red) kinetic energy decay of the radial velocity. (b) Spectral energy plot of the radial velocity fluctuations. (c) Kinetic energy decay of the azimuthal velocity. (d) Kinetic energy ratio E r =E q decay. . . . . . . . . . . . . . . . . . . . . . . . . 140 5.22 Basin resonant periods and modes. . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.23 Basin resonant periods and modes - continued from Figure 5.22. . . . . . . . . . . 142 6.1 (a) The camera’s FOV and the ground control points used to obtain the camera’s extrinsic parameters are shown in the inlet; (b) The Ventura entrance channel, cam- era position and field of view (FOV), and local coordinate systems (labelled as CS 1 and CS 2). The coordinate system y axes are positioned in the centre-line of the entrance channel and its two branches, whereas the x axes define the distance normal to the centre-line. Three zones are defined by distance x: zone 1 with jxj<= 25m, zone 2 with 50m<jxj<= 25m, and zone 3 withjxj> 75m. . . . . . . 146 xiii 6.2 (a)-(d) Time history of channel-parallel velocity in Ventura Harbour, measured (points) versus predicted (curves). Crosses (+) represent GPS measurements, whereas circles (o) and squares () represent optical measurements referenced to coordi- nate systems CS 1 and CS 2 respectively. All GPS measurements are referenced to CS 1. The flow field has been separated in three zones, defined by the distance (x) normal to the channel center-line. The mean flow channel speed u y predicted by MOST is overlaid in all subplots as continuous curves, whereas the dashed curves correspond to the minimum and maximum MOST-predicted values in each chan- nel zone. All results are shown together in (d). Along-channel distance (y) shown with colormap. (e) Ventura’s tide gauge record compared to the surface elevation time history extracted from MOST. The tide gauge record has been high-pass fil- tered using a cutoff period of 100 min. The location of the tide gauge is shown in Figure 6.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 6.3 Rectified camera frames captured at (a) 8:29 and (b) 8:43 AM local time, over- laid over ortho-rectified USGS aerial imagery (times are shown in Figure 6.2 with vertical dashed lines). (a) transitional-jet and TCS forming when flow reversed towards the harbour; (b) two transitional-jets expelled from the channel bifurca- tion corners. Vectors correspond to tracer flow speeds. TCS streamlines and the fronts of transitional-jets are approximately shown with dashed curves. . . . . . . . 149 6.4 Projected tracer paths in the harbor entrance channel. Continuous (-) and dash- dot (-.) curves correspond to paths extracted from optical tracking, referenced to coordinate systems CS 1 and CS 2 respectively. Dashed (- -) curves correspond to GPS tracks, all referenced to CS 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.5 Scatter plot of velocities measured with the DEX GPS unit versus velocities mea- sured optically. Squares correspond to velocity magnitude, crosses to the velocity in the eastern direction, and stars to the velocity in the north direction. The R 2 coefficients of determination are noted. . . . . . . . . . . . . . . . . . . . . . . . . 153 xiv Abstract The aim of this thesis is to study the generation and evolution of tsunami-induced coherent struc- tures in a well-controlled environment using realistic scaling. The results of this study are intended to provide a step forward towards unraveling the complexities involved in how tsunami-induced coherent structures affect our coastal infrastructure. A physical configuration was created in the tsunami wave basin of Oregon State University, in the image of a port entrance. A small-amplitude, long period wave is generated, and the flow is accelerated through the harbor channel. A separated region forms, which when coupled with the transient flow, leads to the formation of a turbulent coherent structure (TCS). The wave forcing translates to realistic tsunami prototype scales. The fully turbulent and subcritical flow conditions suggest that the laboratory results are scale independent and applicable to geophysical flows. Surface velocities are inferred from surface tracers using 2D and 3D PTV techniques. Surface elevation data was collected throughout the basin to provide an overview of the long wave motion. Acoustic Doppler Velocimeter data were also collected to provide external forcing information and data on coupling between the external tsunami forcing and the generated TCS. The experimental layout was optimized to allow a large coherent structure to form in the off- shore basin. The experimental data collection extended long-enough to capture the growth and evolution of the TCS until the TCS radial growth became limited by the vertical boundaries. The circulation growth of the TCS is compared with studies on tidal flushing in near-shore inlets. The structure of the resulting monopolar coherent structure is assessed through its velocity components in cylindrical coordinates, so as to characterize its profile through other known types of geophysi- cal vortices. Kinematic energy decay due to bottom friction and the growth of the vortex structure are compared with analytical models. The three-dimensional motions are examined through the secondary flow components to characterize the two-dimensionality of the shallow flow. xv Chapter 1 Introduction Tsunami science has evolved rapidly during the past three decades, and now accurate real-time forecasting of transoceanic tsunamis is possible (e.g. Titov et al. (2005)). Reaching this point has been a collective effort involving the development of analytical and numerical tools and testing them against well-controlled experimental data (e.g. Synolakis et al. (2008)). Experimental data provide the basis of benchmarking before applying engineering tools to the real world. On the other hand, field survey observations provide incentives to study new aspects of tsunami hydrodynamics and have raised new questions about our current physical understanding of tsunamis. Field surveys, studies of the population and civil defense exercises allow for diachronic evaluation of progress towards effective tsunami hazed mitigation, as discussed in Okal (2015). Among other vexing and still unresolved questions are the sequencing of waves in a tsunami train (Okal and Synolakis, 2015), the late arrival of tsunami waves (Tsai et al., 2013), the overland flow velocities (Fritz et al., 2012) and various harbor excitation phenomena which can not always be explained with existing paradigms for harbor resonance. Earthquake-generated tsunami waves have a long wavelength compared to the local water depth, allowing to describe their dynamics using the shallow water equations. The validity of this approximation has been tested numerous times against both laboratory and field measure- ments of water surface excitation and on-land inundation (e.g. Titov and Synolakis (1997)). Given appropriate initial conditions, it has been shown that tsunami models can predict on-land runup 1 and inundation to satisfactory levels (e.g. Montoya et al. (2016)). This allows to utilize tsunami models to produce tsunami hazard maps (e.g. Barberopoulou et al. (2011a)) that can be used for hazard zoning and evacuation planning in populated coastal areas. Field observations Having established a working methodology to accurately predict tsunami amplitudes and inun- dation, predicting tsunami-induced currents in the near-shore is not equally well understood and lacks extensive research. Reports of unusual excitations in ports and harbors have emerged for numerous tsunami events. Perhaps the best-known example was reported by Okal et al. (2006a) for the port of Salalah in Oman, following the 2004 Indonesian mega-tsunami. Approximately 90 min after the first-wave arrival, the moorings of the 285m freighter Maersk Mandraki broke as it was caught inside a series of tsunami-induced turbulent eddies. The ship was spinning out of con- trol for several hours (miraculously without colliding with other ships or harbor structures) before being left stranded on a nearby sand bar (Figure 1.1). The inferred tsunami-induced currents from the observations were disproportional to the runup measurements in the area, which did not exceed 1:8m. Figure 1.1: The path freighter Maersk Mandraki followed after it was caught in the turbulent eddies that broke its moorings during the the 2004 Indonesian tsunami. The path was approximately hand- drawn by an eyewitness. Source: Okal et al. (2006a) 2 A similar incident was reported by Okal et al. (2006c) about a 196m container ship named MSC Uruguay becoming uncontrollable for several hours in the harbor of Le Port, R´ eunion Island during the same tsunami event. The strong tsunami-induced currents broke its moorings twice, the first time was 1.5 hrs after the period of major vertical fluctuations, and the second about 2.5 hrs later. Possibly the most vexing observation from the 2004 tsunami field survey expeditions was reported by Okal et al. (2006b) for Toamasina harbor in Madagascar. Several hours after the arrival of the waves (approximately 6.5 hrs), of what appeared to be oscillations of the order of 0:6m, considerable amplification of the currents in the harbor were reported, resulting in another incident involving a 50m freighter ship by the name Soavina III. This incident, which is likely related to the late arrival of higher-frequency waves and harbor resonant amplification (Okal et al., 2006b), highlighted the fact that currents in ports and harbors can persist for much longer than on-land inundation. The notion of tsunami amplification in ports and harbors is not new. As the word tsunami (translated to “harbor wave” from Japanese) suggests, tsunamis have created havoc in ports and harbors well before these recent observations. Archaeological evidence suggests that the bronze age eruption and caldera collapse of the Santorini volcano destroyed important naval infrastructure of the Minoans in Crete Island, Greece (Bruins et al., 2008) and it is believed that it had a direct impact on the demise of what was the first western civilization (Driessen and MacDonald, 1997). Pre-historical tsunamis aside, one of the earliest reported examples of a tsunami causing damage to ports and harbors is that of the 1868 Arica, Peru tsunami described in detail by Borrero et al. (2015). The description of other notable tsunami-current events in the modern era are is included in the comprehensive review of Borrero et al. (2015). The last event mentioned in more detail here is the 2011 Tohoku, Japan tsunami. Apart from the extensive inundation - most prominently in the Sendai plain (Mori and Tomoyuki, 2012) - eddies of big proportions were reported in both near and far-field harbors (Lynett et al. (2012, 2014) and ref- erences therein). In the near-field, aerial footage of Oarai harbor taken from a helicopter shows the emergence of large rotational structures approximately 3 hr after the earthquake (Figure 1.2). The structures formed by the interaction of the tsunami-induced currents with the man-made coastline. 3 Topographic forcing from coastal infrastructure features (such as breakwaters and jetties) creates sharp velocity gradients to form a flow instability that develops into large-scale coherent structures that persist for long periods of time. The coherent structures observed in Oarai occupied the whole harbor basin and featured distinctive zones of flow convergence near the center and a zone of diver- gence in the outer vortex region (Lynett et al., 2012). Boats and ships entrained in the flow would converge towards the eye of the eddy that experiences a local surface elevation depression due to a pressure minimum, giving rise to descriptions of the eddy as a “whirpool” that sucks items in its wake. Similar flow patterns were observed in other near-shore areas in Japan, such as in Iwaki city. The picture shown in Figure 1.2 (right), features a dipole structure with two counter-rotating circulations. Figure 1.2: Near-shore rotational flow structures observed during the 2011 Japan tsunami in Oarai harbor (left) and Iwaki city (right). Source: Reuters. Coherent structures in shallow water flows Although the observations of tsunami-induced currents and the emergence of coherent structures during tsunami events are relatively new, such phenomena are well documented in geophysical flows. The term “turbulent coherent structures”, henceforth abbreviated as TCS, refers to Two-dimensional, connected, large-scale turbulent fluid masses that extend uniformly over the full water depth and contain a phase-correlated vorticity, with the exception of a thin near-bottom boundary layer (Hussain, 1983; Jirka, 2001). 4 Coherent structures form in uni-directional fully turbulent shallow flows, with large Reynolds num- ber Re h = UH=n >> 10 3 (where U is characteristic velocity, H is the local depth, and n is the kinematic viscosity) (Jirka, 2001). They are generated by local disturbances that undergo internal oscillations and grow in size, following the concept of the inverse energy cascade of 2D turbulence through which energy is transferred from smaller to larger wavelengths (Kraichnan, 1967). The most common generation mechanism is topographic forcing (in the presence of features such as islands, breakwaters, groynes etc) that leads to flow separation, formation of a shear layer and return velocities in the lee of the feature (Jirka, 2001). Shallow (2D) flows are layered turbulent flows in a domain where the two horizontal dimen- sions - including the dimension along the flow direction - greatly exceed the third (vertical) dimen- sion (Jirka and Uijttewaal, 2004). Such flows are found in nature in rivers, lakes, estuaries, stratified water bodies and coastal regions. It should be stressed at this point that term “shallow flow” refers here to the analogy of the length of a characteristic vortex structure being larger than the water depth and should not be confused with the water-wave shallow-water approximation for which the depth is much smaller than the surface wavelength. Defining a characteristic length scale of the flow structures L and vertical a scale from the local water depth H, it follows that for shallow flows L H >> 1: (1.1) The second requirement for shallow flows is the presence of at least one shear-supporting boundary, while the others may be shear-free or shear-supporting (Jirka and Uijttewaal, 2004). The base flow is governed by wall turbulence produced from the no-slip boundary condition gives rise to small-scale 3D flow features (L 3D < H), while the mean velocity profile can be characterized using the logarithmic-law (Jirka and Uijttewaal, 2004). The small-scale 3D features constitute what is called the “secondary flow” and their energy is suppressed by increasing flow shallowness (Paret and Tabeling, 1997; Sous et al., 2005) but also depends on the Reynolds number (Dolzhanskii et al., 1992; Satijn et al., 2001; Duran-Matute et al., 2010). The flow will eventually grow large- scale structures with 2D character (L 2D >> H) and vertically-aligned vorticity vectors. In their study on plane turbulent jets, Dracos et al. (1992) have found that the flow development can be 5 categorized to three stages according to the distance x from the source, where in the last stage (far-field, x=H 10) the structure of the flow has grown into scales greater than the local depth. At the point where the secondary flow can be considered insignificant, the flow can be characterized as quasi-two-dimensional (Q2D) (Sous et al., 2004). Satellite imagery has revealed a wealth of horizontal structures on various scales in the Earth’s atmosphere and oceans (Fl´ or and Van Heijst, 1994). The observed large-scale vortices of long lifes- pan share the commonality of forming in two-dimensional flow fields (Fl´ or and Van Heijst, 1994). Depending on the initial conditions or forcing, geophysical-scale vortices vary in size, amplitude and radial profile (Mcwilliams, 1990). Two particular types of vortices of interest to this study, which are found in nature in abundance, are the monopolar and dipolar (e.g. Figure 1.2). Other types of vortices such as the tripole (or triangular) may emerge as a perturbation of an unstable, isolated circular vortex (Carnevale and Kloosterziel, 1994), but they are not further discussed here. A dipolar vortex is a compact vortex pair of counter-rotating vortices with axially aligned columns that possesses a self-propelling mechanism (van Heijst and Clercx, 2009). They have been observed to form in the near-shore zone in tidal channels and in the head of rip-currents (Sous et al., 2004). Tidal flows create strong currents at inlet systems and are responsible for the exchange of important nutrients and sediment in the coastal zone (Whilden et al. (2014) and references therein). The flow may separate through topographic forcing at sharp coastal features and lead to the formation of a starting-jet dipole or re-circulating eddies (Figure 1.3). Their formation and propagation characteristics have been studied experimentally by Wells and van Heijst (2003) and Nicolau del Roure et al. (2009) for oscillatory flow, by Afanasyev (2006) for impulsive flow, and their mean turbulent vortex properties were studied by Bryant et al. (2009). Wells and van Heijst (2003) and Afanasyev (2006) used a a two-layer stratification was used to minimize the bottom friction effects. Kashiwai (1984) and Wells and van Heijst (2003) provide theoretical models, based on potential flow theory, to predict under which conditions the dipole escapes the tidal inlet. The circulation growth of the vortex system is estimated using a model developed for vortex ring formation (Shariff and Leonard, 1992). 6 V ortex rings are toroidal vortices typically generated in the lab through piston-cylinder ring generators using either a nozzle or orifice geometry: fluid is impulsively pushed through the cylin- der, forming a boundary layer at the cylinder wall (Figure 1.3). The vortex sheet separates at the nozzle/orifice to form a vortex ring (Didden, 1979). The generation mechanism of vortex rings has striking similarities with the formation of starting-jet vortices in tidal inlets, and in essence vortex rings correspond to the three-dimensional analog of a vortex dipole (Afanasyev, 2006). Figure 1.3: Left: Tidal flow leaving a sharp channel may induce the formation of a dipole by advection of vorticity created at the boundaries (Source: Wells and van Heijst (2003)). Right: Nozzle type piston-cylinder configuration to generate a vortex ring by advection of vorticity created at the cylinder boundary layer (Source: Krueger (2005)). Monopolar two-dimensional vortices are the most common type of vortex found in nature and usually have a circular or elliptical shape (Fl´ or and Van Heijst, 1994). Different generation methods have been developed to generate monopolar vortices in the laboratory. van Heijst and Clercx (2009) provide an extensive review of all the methods. Perhaps the most relevant generation method to this study is the “stirring method”. It involves confining fluid inside a rotating cylinder (with no- slip boundaries) which is lifted once a purely azimuthal flow is achieved inside the cylinder. The initially still ambient fluid interacts with the rotating fluid to create a zone of vorticity with an opposite sign than the vortex core. This technique was used by Seol and Jirka (2010) to generate a monopolar vortex with high initial Reynolds number, in what is one of the few studies on monopo- lar vortices that is extended to fully turbulent flows (to the authors knowledge). Another similar technique to generate stirring-type monopolar vortices is the “injection method” which consists of injecting fluid along the inner wall of an open thin-walled cylinder submerged in the fluid (Flor 7 and Van Heijst, 1996). The process of injecting fluid stops and after allowing for the flow to evolve into regular rotational motion the cylinder is lifted. Both techniques create isolated vortices with Gaussian profiles. This profile description was found to match the experimental observations for the monopolar vortex of this study. Numerical modeling and field measurements of tsunami-induced currents Turbulent phenomena as the flow features observed during the recent tsunamis would ideally be modeled using 3D models that are able to resolve the spatial and temporal scales involved (e.g. Choi et al. (2008); Kim et al. (2013)). On the contrary, due to computational constraints, tsunami currents are typically modeled either through the depth-averaged non-linear shallow water equa- tions (NSWE, e.g. Ayca and Lynett (2016)) or using higher-order Boussinesq-type (BT) models (e.g. Tavakkol and Lynett (2016)). Kim and Lynett (2011) have added additional terms to the standard Boussinesq equations to better resolve horizontal and vertical turbulent mixing, bottom stress and depth averaging effects. Son et al. (2011) have developed a multi-physics approach to nest a BT with a NSWE model, switching to the BT model in the near-shore to better resolve the complex flow dynamics. The non-linear interaction between tsunami- and tidal-induced currents has recently been studied by Ayca and Lynett (2016) for ports and harbors in California. In a recent modeling workshop in February 2015 (Lynett et al. (2017), in preparation), BT models and NSWE models were validated against field observations. It is thus far clear, that for non-research level analyses, eddies that impact ports can be simulated at least as accurately with NSWE as with BT, at least as compared to field observations. The bathymetric grid resolution was found to be of great importance, and that a grid cell-size as small as 2 m is required to resolve the eddies in both BT or NSWE simulations. The estimation of tsunami-induced current statistics takes a much longer computational time than estimation of the maximum runup or inundation. In addition to the fine grid resolution needed to resolve eddies, strong currents can persist for days, and this makes computations more tem- peramental, as numerical errors accumulate. To examine the build-up of of numerical errors in long simulations, Lynett et al. (2014) simulated the 2011 Japan tsunami in Crescent City harbor, 8 California, for 60 hrs (real-time) using a NSWE model and demonstrated that the energy decay of wave amplitude is accurately reproduced. The tsunami current predictions from numerical models over long simulation periods should be approached in a statistical framework (e.g. Lynett et al. (2014)). By combining the modeling results for past tsunami events with recorded damage observations in ports and harbors (e.g. Wilson et al. (2013)) current speeds can be related to damage levels. Using this correlation for ports and harbors for which damage reports exist, it is possible to produce tsunami fragility curves that can be used for mitigation and evacuation planning for future tsunami events (Wilson et al., 2016). Similar methodology has been applied to develop tsunami fragility curves for buildings after the 2011 Tohoku tsunami (e.g. Suppasri et al. (2012)). In terms of field measurements, which are important for physical understanding and numerical model validation, a small but increasing number of current measurements have been reported from acoustic current meters (Lacy et al., 2012; Borrero et al., 2013; Admire et al., 2014) and by track- ing floating objects in eyewitness videos (Fritz et al., 2012; Hayashi and Koshimura, 2013; Admire et al., 2014). In a recent study, Kalligeris et al. (2016) (also included in this thesis as a separate chapter) report Lagrangian flow measurements of tsunami-induced currents in Ventura harbor, Cal- ifornia, during the 2015 Chilean tsunami. They provide a comprehensive description of the flow field, based on quantitative data of tsunami current effects and large scale coherent flow structures unfolding during a real event. Apart from the handful reported recordings of tsunami-induced currents, data to validate numerical models remain sparse. 1.1 Objectives of Study and Proposed Approach This study is concerned with the generation and evolution of tsunami-induced coherent struc- tures in the laboratory using realistic scaling. Coherent structures control, to a great extent, the tsunami-induced hazard and sediment transport in ports and harbors, particularly for small ampli- tude tsunami events. They have a significant impact on the energy cascade as well as on the dynamics affecting the neighboring boundaries. 9 What distinguishes this effort is that we seek to examine the transfer of momentum and energy through the long wave generation, to the coherent structure generation, to the force applied to boundaries (e.g. bed and harbor sidewalls) by asking general questions involving the cascade of energy through multiple scales. To achieve these goals, a physical configuration was created in the tsunami wave basin of Oregon State University, in the image of a port entrance (Figure 1.4). A small-amplitude, long period wave is generated, and the flow is accelerated through the gap between the breakwater and the basin side wall (phases 1&2 in Figure 1.4). The flow rate through the gap is strong enough, such that a separated region forms, which when coupled with the transient flow leads to the formation of a TCS on either side of the breakwater (inshore and offshore). The TCS circulation grows as a result of vorticity influx and there is no influence from an expelled boundary layer. After the channel return flow, the offshore vortex eventually detaches from the trailing jet resulting in a self-propagating coherent structure that grows in size through viscous diffusion (phase 3 in Figure 1.4). Figure 1.4: The experiment in three phases: 1) The wavemaker forward stroke creates a coherent structure on the inshore basin side; 2) the wave gets reflected on the back wall of the basin and the wavemaker retracts to create a strong current through the channel that generates the offshore TCS; 3) The TCS detaches from the breakwater and gets advected. The experimental layout was optimized to allow a large coherent structure to form in the off- shore basin. The experimental data collection extended long-enough to capture the growth and evolution of the TCS until its radial growth became limited by the vertical boundaries of the off- shore basin. The structure of the monopolar vortex is assessed through its velocity components in cylindrical coordinates, so as to characterize its profile through known types of geophysical vortices. The kinetic energy decay timescales are examined through a force-balance model and 10 the bottom layer effects are examined through the secondary flow components to characterize the two-dimensionality of the shallow flow conditions. This particular application has the rare property, at least in tsunami experimentation, of using both an undistorted model spatial scale with a reasonable time scale. A typical port entrance channel is designed for depths on the order of 15 m, implying a 1/27 scale with the 55 cm water depth used in the TWB. The Reynolds number resulting from the flow conditions is in the fully- turbulent range (Re= O(10 5 )) throughout the duration of the experiment, the Froude number is sub-critical (Fr 0:38) and the Rossby number is large (Ro>> 1). Previous studies (Jirka, 2001; Socolofsky and Jirka, 2004; Bryant et al., 2012) have shown that for Fr < 0:4 and Re h > 10 3 the laboratory results are scale independent, and therefore the results of this study are applicable to geophysical flows. An overview of the long wave motion throughout the basin was provided by wave-gauge mea- surement, whereas Acoustic Doppler Velocimeters (ADVs), were also used to provide external forcing information and data on coupling between the external tsunami forcing and the generated TCS. Inside the TCS generation area, a dense array of measurements is used to quantify the sur- face flow as thoroughly as possible. Surface velocities are inferred from optical tracking of surface tracers using Particle Tracking Velocimetry (PTV) techniques. In the channel formed between the breakwater and the side-wall, surface velocities are inferred in more detail using a close-up stereo camera configuration (3D PTV), which allows for the extraction of additional surface elevation information. From this quantitative data-set, the water surface velocities are used to study the eddy circu- lation growth, kinematic energy decay timescales and infer the two-dimensionality of the flow. The measurements can be used for numerical model validation, while the results of the study offer insights to the complex dynamics of tsunami-induced coherent structures in a well controlled envi- ronment. 11 Chapter 2 Experimental Setup The experiments took place in the tsunami basin of Oregon State University’s O.H. Hinsdale Wave Research Laboratory. The 44 26:5m in–plan–view, basin is one of the largest in the world. A schematic of the experimental set-up is shown in Figure 2.1. A breakwater was built across the width of the basin, at a 27 o angle with respect to the wavemaker, allowing a 3:1m gap between the breakwater tip and the side wall. The basin was thus physically divided to an offshore (on the wavemaker side) and an inshore section. The 26:5 long, 0:6 m wide and 0:8 m high breakwater was built with 12” 8” 16” cinder blocks. The sides were covered with white acrylic Plexiglas sheets to make the breakwater impermeable and create a smooth sidewall surface. Visual data were collected only for the part of the basin plan area denoted by the blue-dashed- dot box in Figure 2.1. In the study area, the uni-struts running through the concrete basin floor and side-walls were covered and the surfaces were painted white. The white background offered good contrast with the black surface tracers used for the PTV experiments. It also served well as a background color for the dye experiments. The experiments were completed in a preliminary stage in October 2013 and a final stage in July-September 2014. During the preliminary experiments, one breakwater design and various wavemaker boundary conditions were tested. The goal of the basin response tests was to generate a stable TCS at the tip of the breakwater. Once this was accomplished, various camera configurations were tested to maximise the surface area covered visually. Particles were introduced in the flow 12 study area TCS formation breakwater offshore basin 63 0 wavemaker 26.5m 44m bridge x y ADV1 ADV2 ADV3 CAM A CAM D CAM B CAM C Figure 2.1: The experimental set-up. Dots denote acoustic Doppler velocimetry (ADV) instrument locations, squares denote the location of the overhead HD cameras used for the 2D PTV and dye experiments, blue dashed box denotes the study area of the basin and the red dashed box denotes the area shown in Figure 2.4b. The local coordinate system y-axis starts from the center-line and is positive upwards. The local coordinate system x-axis starts from the wavemaker mean position and is positive towards the basin. while recording, to ensure they can be traced by conventional PTV algorithms and extract surface velocities. Dye visualization trials were also performed to get visual insights about the 3D structure of the TCS. Finally, the entrainment of scaled-down ships into the eye of the “whirpool” was simulated in the basin (replicating the Oarai harbor observations). All the quantitative data were collected during the experiments of 2014. The final stage experimental-schedule also included complementary dye-visualization trials. Table 2.1 provides a summary of all the experimental trials. The same experiment was repeated many times to exam- ine the repeatability of the experiment and collectively obtain a satisfactory density of data. The following sections provide an overview of the set-up and instrumentation used for the experiments. 13 Table 2.1: Experimental trial short description. Experiment # of trials Description-Remarks 2D particle tracking 22 4 HD camera configuration with small FOV overlap. 4cm parti- cles were introduced in the flow as surface tracers. The ADV data were collected during these trials. 3D particle tracking 62 2 HD camera stereo configuration. 8.5mm particles were intro- duced in the flow as surface tracers. The cameras were moved to 8 different positions around the tip of the breakwater, covering 5 trials per position. Dye tests 14 4 HD camera configuration. Dye was introduced to visualize the 3D structure of the flow. Wave gauges 29 16 wave gauges were mounted on the basin bridge to collect sur- face elevation data. The bridge was moved to 27 different posi- tions. 2.1 Boundary conditions, water level and experiment scaling The boundary conditions in the wave basin consist of solid walls in three of the basin’s sides and a wavemaker on the other (Figure 2.1). The piston-type wavemaker consists of 29 individual boards with 2.0 m maximum stroke and 2.0 m/s maximum velocity. It is capable of creating regular (monochromatic), polychormatic and multi-directional waves, with periods ranging between 0.5 to 10 s (OSU website). The wavemaker motion and water level were optimized to generate a stable TCS at the tip of the breakwater. Numerical modeling provided initial parameters, later fine-tuned during preliminary experiments in October 2013. The piston displacement was kept as simple as posible: a slow push forward and a sudden retreat, taking advantage of the full piston stroke. The result is a single asymmetric N-wave pulse (Figure 2.2) with 42 s period. The water level was set at h= 55cm, resulting in a wave length of 97.6 m (approximately twice the length of the basin). For a “typical” prototype harbor channel depth of 15 m, the Froude scaling is given in Table 2.2, resulting into a 2.66 km prototype wavelength and 3.66 min prototype period, well within nearshore geophysical tsunami scales. Thus, this experiment has the rare property of being realistically scaled with respect to both length and time. The aforementioned boundary conditions and water level were used for all experiments described in this thesis. 14 The Froude number, defined as Fr= U max p gh , where U max is the maximum velocity through the channel, was calculated using U max = 0:89 m/s which was measured at the time of first flow rever- sal in the harbor channel (see Section 4.3). The corresponding depth Reynolds number in the channel is Re h = U max h=n= 43705 O(10 5 ), while the Froude number Fr= 0:38, and therefore the flow is both sub-critical and turbulent. Previous studies (Jirka, 2001; Socolofsky and Jirka, 2004; Bryant et al., 2012) have shown that, for Fr < 0:4 and Re h > 10 3 , the laboratory results are scale independent, and therefore we surmise that the results of this study are applicable to full-scale geophysical flows. The Rossby number for a horizontal scale of 10m O(10) and velocity V 1m=s O(1) cor- responds to Ro= V 0:5 f L 2000, where f = 2Wsinf 0 10 4 s 1 is the Coriolis parameter for the rotation rate of the earth (W= 7:2921 10 5 rad/s) and geographic latitude (f 0 45 o at Corval- lis, Oregon). Since Ro>> 1, experiment scaling suggests that the Coriolis force is negligible and there is no background rotation. The results of this study are therefore restricted to highly-turbulent small-scale geophysical flows. Coherent structure formation in ports and harbors during a tsunami (e.g. Lynett et al. (2012)) meets this criteria, and this is the primary focus of this study. In studies concerning tidal flushing, and starting-jet vortices through inlets (e.g. Bryant et al. (2012)), an important parameter which controls the motion pattern of the generated vortex is the Strouhal number defined as: K w = W UT ; (2.1) where W is the width of the channel, U is the channel-average velocity and T is the wave period. The setup of the parent experiment (T 68 s, W = 3:12 m and U = 0:89 m/s, see Section 4.3) results to K w = 0:05, less than the critical value of K w 0:123, above which the formed dipole is attracted back into the channel sink (Wells and van Heijst, 2003). 2.2 Cameras and light Four Panasonic AW-HE60S (Panasonic, 2015) cameras were used both for the PTV and dye exper- iments. They can record HD video and have pan, tilt and zoom capabilities to fine-tune the field of 15 Table 2.2: Froude scaling Variable Unit Model Prototype Scale factor Remarks Depth L 0.55 m 15 m l = 27:3 scaled to water depth h Wave Length L 97.6 m 2.66 km l = 27:3 scaled to water depth h Velocity L/T 2.32 m/s 12.13 m/s p l = 5:2 scaled to SW p gh Wave Period T 42 s 3.66 min p l = 5:2 scaled to SW L= p gh time (s) 0 50 100 150 200 displacement (m) -1 -0.5 0 0.5 1 (a) wavemaker displacement time (s) 0 50 100 150 200 η (cm) -10 -5 0 5 10 (b) free surface elevation @ x=3m, y=2m Figure 2.2: Wavemaker boundary conditions. (a) The piston displacement time history (0 dis- placement corresponds to x= 0 in basin coordinates), starting from the fully retreated position and reaching maximum stroke; (b) the recorded surface elevation near the face of the wavemaker. view (FOV). For all runs, video was recorded with HD resolution (1920 1080) at a 29:97 frame per second rate. For the 2D PTV and dye experiments, the four cameras were mounted on the basin roof steel beams (see Figure 2.3). Maximizing coverage of the basin study area necessitated minimizing the overlap between the camera FOVs. This camera configuration allowed for extraction of spatial information only in the two horizontal dimensions, thus henceforth termed 2D (see Section 3.3). Optimally, it is desirable to achieve maximum coverage with the highest possible camera resolution (pixels/reference length) at the water surface. Different pan-tilt-zoom configurations were used to study different aspects of the flow (that require different FOVs and resolutions). For the 2D PTV experiments, the chosen camera zoom was such to achieve at least 6pixels/particle diameter (1.5pixels/cm for the s = 4cm diameter tracers used for the 2D PTV experiments) at the water surface, which is an adequate resolution for the tracer recognition algorithm. For the 3D PTV experiments, two Panasonic AW-HE60S cameras were mounted on the basin’s bridge–carriage (Figure 2.4a). The bridge spans over the width of the basin and is sitting on 16 wavemaker CAM A CAM D CAM B CAM C Figure 2.3: Camera fields of view (FOVs) used for the 2D PTV experiments. rails allowing it to be translated along the basin’s long axis. Contrary to the 2D PTV camera positioning, the FOV overlap was maximized for the 3D PTV experiments to achieve a stereo camera configuration; stereo image pairs allow for the extraction of spatial coordinates in all three dimensions. The cameras were positioned 3.27 m above the water surface, achieving a more detailed resolution at the water surface (7 pixel/cm), at the cost of covering a smaller area compared to the camera FOVs of the 2D PTV experiments. Smaller diameter surface tracers were used for the stereo camera configuration experiments (diameter s = 8:5mm) to increase the resolution of the velocity field. The cameras were moved to eight different positions to cover the area around the basin channel (shown in Figure 2.4b), by either translating the bridge (in the x-direction) or by moving the camera mounting (in the y-direction). The camera shutter speed choice was controlled, as always, by the available light. The side windows of the basin were covered with drapes to limit reflections on the water surface. The overhead light panels also created intense glares on the water surface, making it difficult to trace particles going through them. Therefore, rectangular reflective panel (cardboard panels covered with aluminum foil) were hung underneath the lights, to direct the light towards the white ceiling, through which the light was diffused into the room. The reflective panels softened the intensity of the glares but significantly reduced the amount of light reaching the water surface (Figure 2.5). 17 (a) 3D PTV camera setup. pos 1 pos 2 pos 3 pos 4 pos 5 pos 6 pos 7 pos 8 4m (b) Overlapping 3D PTV FOVs. Figure 2.4: 3D PTV experimental setup. The camera shutter speed was adjusted to achieve the brightest possible image, without capturing shadows of the motion of the particles in the still frames. (a) Water surface with intense overhead light glares. (b) Water surface with softened overhead light glares. Figure 2.5: Basin light adjustments. For the 3D PTV experiments, where the tracer velocities were higher (in terms of inter-frame pixels/time), additional light sources were used to increase the camera shutter speed. Two flood lights were used, mounted below the bridge and the side wall (Figure 2.4a). Custom-made light diffusers were used to diffuse the light and soften the light intensity. 18 2.3 ADV’s and wave gauges Complementary quantitative data were collected using Acoustic Doppler Velocimetry (ADV) and wave gauges. The former comprised of three Nortek Vectrino ADV devices, two of them mounted at the basin-floor unistruts and one on a side-wall unistrut near the breakwater tip (see locations in Figure 2.1). The Vectrinos (Nortek, 2015) are high-resolution acoustic velocimeters used to measure 3D water velocity in the laboratory and field. Each acoustic sensor has one transmitter and four receiver transducers, while the sampling volume has a diameter of 6 mm and is located 5 cm away from the sensor so that it doesn’t interfere, as much as possible, with the receiver signal. Velocities are calculated by processing the echo, and finding the Doppler shift as particles pass through the volume, given the speed of sound in the medium. Due to the clarity of the fresh water used in the laboratory, fine seeding material was introduced in the water to provide a surface where the sound signal bounces off. ADV data were collected over four trials at 50 Hz sampling frequency, after which ADV1 (see Figure 2.1) was removed, since its mounting appeared to significantly affect the flow. The sensors of all three ADV’s were positioned mid-depth in the water column (z= 0:275m, see Figure 2.6a). Surface elevation data were collected using 16 p-frame wave gauges mounted on the basin’s instrumentation bridge. The gauges were hanging off the side of the bridge, with only the stainless steel wires being in contact with the water. The bridge was translated in-between 27 trials to cover the whole basin area (Figure 2.7), the two last trial gauge positions (inset in Figure 2.7) adding more density to the surface elevation measurement’s grid. More detailed data were collected in the study area, where both the bridge translation and gauge spacing were the smallest. Each wave gauge has a p-frame supporting stainless steel wires, held apart by an insulated material, which is in turn attached to an aluminum mounting rod (similar to the wave probe shown in Figure 2.6b). Electric current is applied trough the wires, which form a resistor, the load of which depends on the level of immersion in the conducting fluid, in a classic Wheatstone bridge configuration. A simple calibration prior to the experiments provided the relationship between 19 (a) The Vectrino sensor and mounting of ADV 3. (b) Wave probe. Source: www.edesign.co.uk Figure 2.6: ADV and wave gauge used for the experiments. voltage and immersion level. V oltage measurements were collected for 30 min at a 50 Hz sampling frequency. 20 Figure 2.7: Positions of gauges recording water surface elevation. The inset shows the position of the wave gauges during the last two trials (red dots). 21 Chapter 3 Post Processing 3.1 Image pre-processing Flow tracking is possible only through identifying individual particles in images and correlating them with subsequent images. To identify the surface tracers in the camera still frames, first the images have to be pre-processed to enhance their appearance. The pre-processing of the still frames of each trial video file includes the following steps: 1. Remove image distortion by applying the a priori-determined camera calibration coefficients (see Section 3.1.1). 2. Create and apply a mask file that blackens out any features outside the perimeter of the waterline (sidewall, breakwater etc). 3. Convert RGB images to grayscale and enhance the intensity of particles. The third step involves converting the RGB image format to grayscale images. HD camera still frames are integer-valued matrices with dimensions 1080 1920 3, where the third dimension corresponds to the three color intensities (taking values that range between 0 to 255) that make 22 each pixel, i.e. Red(R), Green(G) and Blue(B). Typically, a color-weighed method is applied for the grayscale conversion. Matlab’s rgb2gray routine applies the following algorithm: grayscale intensity= 0:2989R+ 0:587G+ 0:114B; (3.1) but since we want the black particles to appear as white, intensities of Equation 3.1 are inverted (255grayscale intensity). The conversion is applied to all undistorted and masked frames of each trial video file. The mean grayscale image is obtained by averaging all the values each pixel is taking throughout the duration of the video. The grayscale frames are then subtracted from the mean image to obtain an image that contains less noise on the water surface (basin floor features, surface glares and reflections are averaged out). To further enhance the intensi- ties of surface tracers in the image, the intensity values of the final image are scaled up by using intensity= intensity 255=max(intensity) to reach the maximum intensity for the brightest sur- face tracer pixels (provided at least one tracer is captured in the FOV). The steps are illustrated in Figure 3.1. 3.1.1 Camera calibration - intrinsic parametrization Extrinsic parameters are related to the location and the field of view of the camera, whereas the intrinsic relate to the properties of the lens and the camera. The lens of each camera creates a dis- torted image of the FOV . Objects that are linear will appear curved in the image. The pinhole cam- era model, which describes the relation between camera and world coordinates (see Section 3.5), assumes no image distortion. This may be true for certain high end camera lenses, but in general lenses produce non-negligible image distortion. In cameras with zoom lenses (varying focal length and FOV), the distortion is highly zoom-dependent. The internal camera model of Bouguet (2014) is adopted, which is very similar to the one used by Heikkila (2000). It solves for both the radial (6 th order) and tangential distortion coefficients. All internal parameters are obtained by processing calibration images of a flat checkerboard, taken 23 Figure 3.1: Image pre-processing. (a) The original RGB still frame, as extracted from the video; (b) lens distortion is removed; (c) converted to grayscale and applied mask; (d) subtracted from mean video grayscale intensity and enhanced surfaced tracer intensity. at different angles (Figure 3.2). The process was repeated for all zoom levels used for the exper- iments. The calibration toolbox determines the skew and distortion coefficients, the focal length and the principal point of the lens. 3.2 Particle identification The next step is to apply an automated algorithm that detects the center of each surface tracer in the processed image. The methodology of Crocker and Grier (1996), developed to trace colloidal suspensions trough digitized video microscope images, is adopted. The original code of Crocker and Grier (1996) was written in the IDL (Interactive Data Language) programming language, but packages in other environments are available online, including the Matlab version of Kilfoil and Pelletier (2015) which was used here. An overview of the particle detection methodology and an “accuracy” analysis is provided below. 24 Figure 3.2: Camera calibration images used to obtain camera intrinsic parameters. Example shown for one zoom level. Since the particles to be traced have the brightest intensities in the processed grayscale image I, the first step is to find all the image pixels above a user-defined threshold intensity value I min . However, to look for particle center candidates, it is more appropriate to look for local intensity maxima above a threshold value. To avoid picking up noise in the background, more that one detections within a user-defined radius w are rejected, where w is an integer greater than the particle diameter (in pixels) and smaller than the inter particle distance. Crocker and Grier (1996) show that the position of the local maxima and the actual particle center lie within a pixel to each other, for their laboratory conditions. To increase to subpixel accuracy, the center of mass of the integrated intensity is sought within a circle of radius w (the choice of integer w is related to the size of the particles) around the local peak positioned at(x p ;y p ). Integrated intensity (also called the mass or the zeroth moment) is given by (Crocker and Grier, 1996): m 0 = å i 2 + j 2 w 2 I(x p + i;y p + j); (3.2) 25 where i; j are integers representing pixel column and row numbers respectively. The offset of the center of mass is then given by the same as, 0 @ e x e y 1 A = 1 m 0 å i 2 + j 2 w 2 0 @ i j 1 A I(x p + i;y p + j); (3.3) while the updated position of the particle center becomes x= x p +e x ;y= y p +e y . If eitherje x j or je y j are greater than 0.5 pixels, Equation 3.3 is re-applied based on the updated position (Crocker and Grier, 1996). For Equation 3.3 to be valid, the image intensity distribution of the disk around the center of the particle has to be axisymmetric. This is not the case for the unfiltered image of Figure 3.1(d), mostly due to uneven room lighting and also due to nonuniform sensitivity among the camera’s pixels (Crocker and Grier, 1996). To alleviate this, Crocker and Grier (1996) recommend the application of a filter that subtracts the long-wavelength modulation of the background brightness, modeled by a boxcar average over a region of extent 2w+ 1, and suppresses the digitization noise in the CCD camera. The convolution kernel is given by (Crocker and Grier, 1996): K(i; j)= 1 K 0 1 B exp i 2 + j 2 4 1 (2w+ 1) 2 ; (3.4) where B=[å w i=w exp(i 2 =4)] 2 and K 0 = 1=B[å w i=w exp(i 2 =2)] 2 (B=(2w+ 1) 2 ) are normal- ization coefficients. A particle center detection example for two particles of Figure 3.1(d) is given in Figure 3.3, for both the filtered and unfiltered images. It can be seen that the background noise surrounding the particles has disappeared in the filtered image and the intensity distribution of the particles became Gaussian. To establish the accuracy of the particle center detection algorithm, the inter-frame displace- ment of a certain particle is tracked, using both the filtered and unfiltered images as input (see Figure 3.4). For both cases, the inter-frame displacement contains high frequency noise, represent- ing a horizontal fluctuation of the surface tracer certainly not occurring during the experiments. To measure the amplitude of the noise, a low-pass Butterworth filter is applied to the inter-frame displacements time series to obtain a smoother curve (red curve of top two plots in Figure 3.4). The 26 (a) Unfiltered image. (b) Filtered image (using Equation 3.4). Figure 3.3: Particle center detection using a local peak above threshold algorithm (red crosses) and by applying the brightness weighted centroid offset of Equation 3.3 (green crosses) using w= 4. absolute difference between the two curves corresponds to the fluctuation. The mean value of the fluctuation for the example shown in Figure 3.4 is 0:046 and 0:049 pixels, with 0:039 and 0:041 standard deviation for the filtered and unfiltered input images respectively. The same experiment was repeated for 160 other particles, yielding a mean fluctuation amplitude value of 0:043 pixels and 0:045 pixels for the filtered and unfiltered input images respectively. It is concluded that the benefit of applying the boxcar and background noise filter of Equation 3.4 is marginal. Other geometric characteristics of the candidate features can be used to further distinguish sur- face tracers with background noise, or particles forming groups. The Matlab feature-routines of Kilfoil and Pelletier (2015) compute the eccentricity and radius of gyration of candidate features, for which the user can define thresholds and refine the selection process. These parameters are par- ticularly important for particles that appear skewed in the corners of the FOVs and to distinguish single particles from groups of particles. Efforts were made during the experiments to prevent par- ticles from conglomerating, but they were only partially successful. To exclude grouped particles from the analysis, the feature geometry parameters were set to detect only single particles (see Figure 3.5). The final output of this stage is an extracted list of surface tracer image coordinates, for all frames and cameras. 27 frame # 0 500 1000 1500 2000 2500 3000 3500 displacement (pixels) 0 0.5 1 1.5 Inter-frame displacement fluctuation for filtered image frame # 0 500 1000 1500 2000 2500 3000 3500 displacement (pixels) 0 0.5 1 1.5 Inter-frame displacement fluctuation for unfiltered image frame # 0 500 1000 1500 2000 2500 3000 3500 pixels 0 0.1 0.2 0.3 Abs difference in pixels coordinates dx dy Figure 3.4: Inter-frame displacement fluctuation using both the filtered and unfiltered images to detect the center of a random surface tracer (upper two figures): blue curves show the raw inter- frame displacement and red are smoothed using a low-pass Butterworth filter. The bottom figure shows the absolute difference in position (x- and y-coordinates) of the detected particle center between the two input images for each frame number. 3.3 2D Particle Tracking Velocimetry (PTV) Having established a file with the positions r(r;t) of surface tracers, the next task is to identify which particles in time step n correspond to particles in the next step n+1, and hence create particle trajectories. For the 2D PTV tracking, a minimum displacement method is applied to achieve the task. Particle P i;n is expected to lie within a radius L from its previous position. All particles within that radius in the next frame constitute candidates, in what is called the screening process. Finally, matching, or selection, is solely based on the overall minimum inter-frame displacement. The 28 Figure 3.5: The feature detection algorithm identifies single particles (green circles) and rejects groups of particles (red dots). application of this method relies on the fact that the maximum inter-frame displacement is small, compared to the mean inter-particle spacing. The algorithm used is based on the work of Crocker and Grier (1996), written in Matlab script by Kilfoil and Pelletier (2015). It utilizes the theory developed by Einstein (1905) for Brownian motion of suspended particles in a fluid. It can be shown that the probability of a particle undergo- ing diffusion, being displaced by a distance d within the time interval t, is given by (Qian et al., 1991; Crocker and Grier, 1996): P(d;t)= 1 4pDt exp d 2 4Dt ; (3.5) where D is the particle’s self-diffusion coefficient. For the case of surface tracers in this experiment, D can be estimated from the Einstein-Smoluchowski equation D= <jr(t+t)r(t)j 2 > 4t (Crocker and Grier, 1996). The probability for N non-interacting particles (independent motions) is given by the N-product of Equation 3.5 (Crocker and Grier, 1996): P(fd i g;t)= 1 4pDt N exp N å i=1 d 2 i 4Dt ! : (3.6) 29 The code’s objective is to find particle combinations that minimize the sum of the distances covered by the particles within a time step (å N i d 2 i ). Equivalently, we seek to maximizeP(fd i g;t). Pairing all particles with each other would require O(N!) computations, which is unreasonably large. Therefore, the code confines the pairing to a radius L, where L should be somewhat smaller than the typical inter-particle spacinga. Optimally, Crocker and Grier (1996) recommend thatd < L<a=2. This poses a limitation to the video frame rate to achieve convergence of the algorithm without resorting to difficult particle combinatorics. For the current experimental setup, L= 10 pixels yielded optimal results. The algorithm also keeps track of particles that go missing in-between frames, either because they reached the boundaries of the FOV , or because they were not traced/tracked. The last known position of missing particles is kept in memory for a certain number of steps. In order to solve Equation 3.6, they are assigned a displacementd i = L for each time step. If a new particle appears sufficiently close to the predicted location, the trajectory of the missing particle is continued. After the application of the PTV algorithm, a list is created containing particle centers (in image coordinates) with assigned surface tracer and frame numbers. An example of surface tracer paths extracted for camera D is shown in Figure 3.6. The next step is to convert the particle center image coordinates to Cartesian using Equation 3.21 and extract information in physical units. This is the subject of Section 3.5.1. 3.4 3D Particle Tracking Velocimetry (PTV) The camera setup for the 3D PTV experiments is different: only two cameras are used at the same time in stereo configuration and both cameras are located closer to the water surface (see Sec- tion 2.2). The tracers have a smaller diameter (s = 8:5mm) compared to the tracers used for the 2D PTV experiments (s = 4cm). In the captured still images, the apparent particle size (s 6 pix- els) has remained almost the same, meaning the system magnification M s=s has increased by 5 times. System magnification directly influences the inter-frame maximum expected displace- ment. For the same maximum water surface velocity u max = 1:5m=s and the same camera frame 30 Figure 3.6: Selected surface tracer paths projected on the mean grayscale FOV intensity. rate f = 29:97Hz, as in the 2D PTV experiments, the maximum inter-frame particle displacement becomes 5 times greater. A measure of the difficulty to track particles is given by the ratio of inter-particle distance over maximum expected inter-frame displacement as discussed in (Malik et al., 1993): p= D 0 u 0 Dt ; (3.7) whereD 0 is the average inter-particle spacing over the entire frame and u 0 Dt is the distance cov- ered between successive frames. When p>> 1, it becomes relatively easy and accurate to track particles using conventional PTV methods. As p approaches unity (p! 1), it becomes difficult to establish the inter-frame particle correspondence. For the 2D PTV setup, the average inter- particle distance is computed asD 0 100, which leads to p 10. On the other hand, for the 3D PTV setup, the average inter-particle distance is computed as D 0 125, which results to p 3. As a result, the conventional PTV methodology described in Section 3.3 (minimum inter-frame displacement matching) fails to track the surface tracers. Instead, the V orono¨ ı particle tracking method is adopted, which allows for identifying larger inter-frame tracer displacements. 31 3.4.1 Vorono¨ ı particle tracking method Continuing from Section 3.2, given the surface tracer positions in timer(r;t), the task is to identify which particles in time step n correspond to particles in the next step n+ 1. Compared to the minimum displacement method applied in the 2D PTV experiments, for the 3D PTV one more element is added in the selection process, the topology of the neighboring particles of a candidate. An elegant and efficient way to define the topology is through V orono¨ ı cells. Capart et al. (2002) proposed a methodology to utilize properties of V orono¨ ı cells for PTV and it was successfully applied to liquid-granular flows where the inter-frame particle displacement is comparable to the mean inter-particle distance and the flow is rapidly deforming. This methodology has been adopted for the current 3D PTV analysis. The Delaunay triangulation DT(P) of a set of points P is such that the circumcircle of any triangle in DT(P) contains no point of P in its interior (Lee and Schachter, 1980). The triangulation implementation involves increasing the spatial dimension of the set P and computing the convex hull. Since the convex hull is unique for a set of points, so is the Delaunay triangulation. A V orono¨ ı diagram, which is a dual graph of Delaunay triangles, consists of cells defining the area closest to a seed (tracer in this case) (see Figure 3.7a). V orono¨ ı cells have many useful properties that find a wide range of applications. The property exploited in this application is the geometry of the V orono¨ ı 1-star. Given i= 1;2;3:::;N seeds at time frame n, the V orono¨ ı 1-star S i;n is formed by connecting a seed P i;n with the seeds of the neighboring V orono¨ ı cells (see Figure 3.7b). The position of the vertices relative to its center vertex, defines the local topology which can be used as a matching template. In many cases, the Delaunay edges connected to a seed will match the V orono¨ ı 1-star edges. However, Delaunay triangles are more sensitive to changes in the positions of the neighboring seeds and therefore are deemed unsuitable for this application (Capart et al., 2002). Between two successive frames n= 1;2, candidate matches P j;2 for particle P i;1 are sought in an area defined by a radius L around P i;1 , where L is the maximum expected inter-frame displace- ment. For all candidates P j;2 within radius L, a corresponding V orono¨ ı 1-star S j;2 is formed. The 32 P i,n (a) Network of Delaunay triangle edges (- -) connect- ing nearest (particle) neighbors, with vertices P i , and the corresponding V orono¨ ı cells (–). P i,n S i,n (b) The V orono¨ ı 1-star S i (-), for a particle with vertex P i , and its corresponding V orono¨ ı cell (- -) defining the area closest to P i . Figure 3.7: Symbols and definitions for V orono¨ ı polygons. coordinates of the S j;2 vertices are then translated so that the center vertex coincides with P i;1 . The goodness-of-fit of any V orono¨ ı 1-star pair S 1 ;S 2 is evaluated using a distance measure which Capart et al. (2002) define as: dist S (S 1 ;S 2 )= median k 1 =1::m 1 min k 2 =1:::m 2 j(r k1;1 r 0:1 )(r k2;2 r 0;2 )j ; (3.8) where k 1 = 1::m 1 and k 2 = 1::m 2 denote the vertices of each V orono¨ ı 1-star S 1 and S 2 and k= 0 corresponds to the central vertex of the star. The function in the brackets computes all vertex pair distances and identifies the vertex correspondence by seeking to minimize the total distance. The steps this far in the algorithm are shown schematically in Figure 3.8. After calculating the distance measure for all possible pairings, the matching of particles between successive frames is achieved by minimizing the global objective function (Capart et al., 2002): min(n 1 ;n 2 ) å dist S (S i(k);1 ;S j(k);2 ); (3.9) 33 Figure 3.8: The particle screening (left plots) and selection (right plots) processes of the V orono¨ ı PTV method. (a) candidate matches P j;2 () for particle P i;1 (O) are sought in an area defined by a radius L around P i;1 . (b) remaining candidate particles P j;2 () after screening. (c) the V orono¨ ı 1-stars of a candidate S j;2 and its pair S i;1 . (d) V orono¨ ı 1-star S j;2 is translated to match the central vertex of S i;1 and apply Equation 3.8. where n 1 and n 2 are the number of V orono¨ ı 1-stars in the examined successive frames. Both the vertex and particle matching functions are solved using the so called Hungarian cost-optimization algorithm, also known as Munkres’ algorithm (Bourgeois and Lassalle, 1971). If a vertex or parti- cle is missing its pair, the algorithm doesn’t assign a pair to them. Capart et al. (2002) suggest a filtering method to remove erroneous particle matches. It is based on the definition of a star-averaged velocity of particle P i as the weighted sum: v i = m å k=1 l k v k ; (3.10) where l k are the weights assigned to the m neighbouring vectices with velocity v k . The weights are computed using the inverse distance: l k = 1 jr k;i r 0;i j : (3.11) 34 Defining a second–order spatial differential operator: D 2 v i = v i v i å m k=1 jr k;i r 0;i j (3.12) allows the imposition of a user-defined limit to local curvaturejD 2 v i j< tol and filter out spurious matches. Having individually tracked the particles for each camera, the next step is to stereoscopi- cally pair them. This is accomplished with the ray pairing method, which is the subject of the next subsection. 3.4.2 Ray intersection and matching algorithm Assuming the two cameras are perfectly synchronized, the pairing of particles appearing in the FOVs of the two cameras (Left and Right) is obtained using the ray intersection and matching algorithm proposed by Douxchamps et al. (2005). Two candidate features for a stereoscopic match, P (L) i and P (R) j are located along 3D rays with parametric equations (see Equation 3.29): r (L) i (a)= p (L) i +aq (L) i ; r (R) j (b)= p (R) j +bq (R) j ; (3.13) where the free parametersa andb are unknown. Referring to Figure 3.9, we seek to minimize the distance between the two rays. This is achieved by solving fora andb in the system 2 4 q (L) i q (L) i q (L) i q (R) j q (R) j q (L) i q (R) j q (R) j 3 5 2 4 a b 3 5 = 2 4 q (L) i (p (R) j p (L) i ) q (R) j (p (R) j p (L) i ) 3 5 ; (3.14) which dictates that the line segment linking the edges of vectors r (L) i and r (R) j is orthogonal to both rays. The midpoint r i; j and distance d i; j of the line segment are defined as r i; j = 1 2 (r (L) i (a)+ r (R) j (b)); d i; j =jr (R) j (b) r (L) i (a)j: (3.15) The global matching problem is solved by seeking to minimize the objective function min(n 1 ;n 2 ) å d i; j (3.16) 35 by imposing the constraint that any given particle can be assigned one matching. V ogel’s approx- imate matching algorithm is used to solve the objective function: each particle i is assigned a best and second best match (d i; j1 and d i; j2 ) and the matching order is dictated by the difference jd i; j1 d i; j2 j in descending order. A filter similar to Equation 3.12 is applied to filter out spurious matches. The final world coordinates are obtained from Equation 3.26 in Seciton 3.5.2. After obtaining the pairing of particles between the two cameras, and the resulting world coor- dinates, the corresponding 2D image–coordinate paths from Equation 3.9 can be linked together to compute inter-frame displacements in all three directions. If the cameras are not perfectly syn- chronized, the 2D image–coordinate paths can also be used to interpolate the position of a given particle in time. Figure 3.9: Intersection of two rays associated with the left and right camera viewpoints. Source: Douxchamps et al. (2005) 36 3.5 Extrinsic parametrization - coordinate transformation Assuming the camera lens does not introduce distortion to the image plane (pinhole camera model), objects with world coordinates (x;y;z), relative to an arbitrary 3D Cartesian system can be trans- lated to image coordinates (u;v) using the collinearity equations (Holland et al., 1997): u u 0 =C u m 11 (x x c )+ m 12 (y y c )+ m 13 (z z c ) m 31 (x x c )+ m 32 (y y c )+ m 33 (z z c ) ; (3.17) v v 0 =C v m 21 (x x c )+ m 22 (y y c )+ m 23 (z z c ) m 31 (x x c )+ m 32 (y y c )+ m 33 (z z c ) ; (3.18) where C u = f=l u and C v = f=l v are coefficients relating the scale factorsl u andl v , to the effective focal depth f . x c ;y c ;z c are the coordinates of the optical center and m i j are elements of a 3 3 orthogonal matrix related to the azimuth, tilt and roll of the camera (see Holland et al. (1997) for more details). The pinhole camera model assumes that the object, the optical center and the image point are all positioned along a straight line. The collinearity equations can be combined to form the following set of equations (Holland et al., 1997): u= L 1 x+ L 2 y+ L 3 z+ L 4 L 9 x+ L 10 y+ L 11 z+ 1 ; v= L 5 x+ L 6 y+ L 7 z+ L 8 L 9 x+ L 10 y+ L 11 z+ 1 ; (3.19) where L 1 ;L 2 ;:::;L 11 are the Direct Linear Transformation (DLT) coefficients. If the DLT coeffi- cients are known, the transformation from world to image coordinates is direct. The opposite is not true. Inverting the system of equations leads to the following (Holland et al., 1997): 2 4 L 1 L 9 u L 2 L 10 u L 3 L 11 u L 5 L 9 v L 6 L 10 v L 7 L 11 v 3 5 2 6 6 6 4 x y z 3 7 7 7 5 = 2 4 u L 4 v L 8 3 5 ; (3.20) which is underdetermined. In the simplified case, of an area of interest that is relatively flat (e.g. water surface with small amplitude waves), the third dimension can be eliminated (z= 0) and Equation 2.5 becomes: 2 4 x y 3 5 = 2 4 L 1 L 9 u L 2 L 10 u L 5 L 9 v L 6 L 10 v 3 5 1 2 4 u L 4 v L 8 3 5 ; (3.21) 37 which gives a direct correspondence from image to world coordinates, provided DLT parameters L i of Equation 3.21 are known. This simplification makes the coordinate transformation possible with only one camera. Therefore, the FOVs of the cameras do not need to overlap, allowing a greater surface area coverage for a given number of cameras. There are two disadvantages with this simplification: i) loss of accuracy in the coordinate transformation, ii) the elevation information of the surface tracers are not extracted. The second drawback is a direct consequence of omitting the elevation from the coordinate transformation. The first drawback, on the other hand, arises from the fact that the surface is assumed to be flat when it isn’t (see Section 3.6). 3.5.1 2D rectification In the 2D PTV experiments, and in order to determine the DLT coefficients for each of the four cameras, well-distributed ground control points (GCPs) were laid on the dry basin floor, before the start of the experiments (Figure 3.10), i.e., evenly distributed in the FOV and with control points near the corners of the FOV to better constraint the transformation. GCP targets were set on the basin side wall and on both sides of the breakwater, at an elevation equal to the still water level. Extra GCPs were set on vertical poles at three different elevations. The x;y;z (basin) coordinates of the GCPs were measured using a total station. The middle row GCPs corresponds to the water level (z= 0:55cm), from which the DLT coefficients of Equation 3.21 are obtained. These are the necessary coefficients to calculate the basin coordinates of the surface tracers with no overlap between the camera FOVs. 38 Figure 3.10: FOV from camera B showing the GCPs layed out on the (dry) basin prior to the 2D PTV experiments. Re-writing Equation 3.21 in the form Ax= b, where x is the vector with the unknown DLT parameters results into, Ax= b=> 2 4 x i y i 1 0 0 0 u i x i u i y i 0 0 0 x i y i 1 v i x i v i y i 3 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 L 1 L 2 L 4 L 5 L 6 L 8 L 9 L 10 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 4 u i v i 3 5 ; (3.22) which corresponds to ground control point (GCP) i with known image (u i ;v i ) and basin (x i ;y i ) coordinates. Given at least four well distributed (in the FOV) GCPs, unknown vector x of Equa- tion 3.22 can be obtained using a least square method (e.g. the Optimization Toolbox of Matlab r ). 39 The obtained DLT coefficients are applied to transform the image coordinates of the detected sur- face tracers to Cartesian. They may also be applied to rectify an image, an example of which is shown in Figure 3.11. Figure 3.11: Image in Figure 3.10 rectified using the calibration coefficients obtained from Equa- tion 3.22. The image acquires the right coordinates for any given point on the plane z= 0:55m. GCP basin coordinates are marked with the red crosses and the transformed coordinates from Equation 3.21 are shown with blue circles. Once the 2D DLT coefficients of Equation 3.22 are computed, the accuracy of the image to the transformation to world coordinates can be assessed by comparing the predicted world coordinates of the GCPs (using Equation 3.21) to the corresponding world coordinates, as measured with the 40 total station. The error e xy is given by the mean distance between the real (x;y) and predicted (x p ;y p ) coordinates: e xy = 1 N N å i=1 q [x p (i) x(i)] 2 +[y p (i) y(i)] 2 ; (3.23) where N is the number of GCPs used for the extrinsic parametrization. Unfortunately, not all the GCPs set up on the dry basin floor prior at the start of the experiments were used. During the course of the experiments, the cameras inadvertently moved slightly from their original positions. Other than the fixed side-wall and breakwater GCPs, a selected number of moveable GCPs were introduced in the FOVs three times a day, so as to update the extrinsic parametrization. After each parametrization, the resulting DLT coefficients and associated errors were slightly different. For the 2D PTV trials, the most recent camera-calibrated DLT coefficients were used to transform the surface tracer coordinates. The number of GCPs used for the extrinsic parametrization and the accuracy of the coordinate transformation for each camera are given in Table 3.1. Note that Table 3.1 shows the mean values from all extrinsic parametrizations. Table 3.1: Mean and max coordinate transformation error for each camera. Camera No of GCPs (N) e xy (cm) maxe (cm) Standard deviation A 7 0.58 1.30 0.0038 B 15 0.43 1.11 0.0028 C 17 0.87 1.87 0.0045 D 13 0.54 1.02 0.0026 The GCPs at elevations different than the water level were used to determine the remaining DLT coefficients (L 3 ;L 7 ;L 11 ), which are related to the third dimension, i.e. the surface elevation. These coefficients are necessary to obtain the elevation of the surface tracers in the 3D PTV experiments, and also to calculate the optical center coordinates x c ;y c ;z c of each camera, which proved to be 41 useful in the post-processing of the 2D PTV data (see Section 3.6). For all three dimensions, Equation 3.22 is: 2 4 x i y i z i 1 0 0 0 0 u i x i u i y i u i z i 0 0 0 0 x i y i z i 1 v i x i v i y i v i z i 3 5 2 6 6 6 6 6 6 6 6 6 6 6 6 6 4 L 1 L 2 : : : L 11 3 7 7 7 7 7 7 7 7 7 7 7 7 7 5 = 2 4 u i v i 3 5 ; (3.24) which corresponds to GCP i with known image (u i ;v i ) and basin (x i ;y i ;z i ) coordinates. For the 3D case, at least six GCPs are needed to determine the eleven unknown DLT coefficients. The optical center coordinates can then be calculated as: 2 6 6 6 4 x c y c z c 3 7 7 7 5 = 2 6 6 6 4 L 1 L 2 L 3 L 5 L 6 L 7 L 9 L 10 L 11 3 7 7 7 5 12 6 6 6 4 L 4 L 8 1 3 7 7 7 5 : (3.25) 3.5.2 3D rectification For the 3D PTV experiments, a calibration table was used to obtain the camera extrinsic parameters (Figure 3.12); extrinsic parameters are related to the location and the field of view of the camera, whereas the intrinsic relate to the properties of the lens and the camera. The calibration table featured GCPs at two different elevations: the lower elevation GCPs ( 1cm above still water level) were in the form of a square-celled grid painted on the board, and the upper GCPs ( 10cm above still water level) were drawn on 4 4’s which were attached to the calibration board. Ideally, the water level would be in-between the two GCPs elevations. However, this could not be achieved, since the water level of the basin would have to be lowered during every position calibration session, to allow for the lower (board) GCPs to be over the water, and be clearly visible in the FOVs. 42 Every time the cameras were moved to a new position, the calibration table was placed in the basin, centered in the camera FOVs, and its four corner coordinates were measured using a total station. The relative positions of the rest of the GCPs were measured (with a total station) on a dry floor prior to the start of the experiments. The four corner reference coordinates were used to obtain the absolute basin coordinates of all the GCPs, in each camera position, using a rigid body coordinate transformation, i.e. translation and rotation, keeping the relative coordinates between the GCPs intact. Figure 3.12: Calibration table used to obtain camera extrinsic parameters for the 3D PTV experi- ments. The bottom row GCPs correspond to the intersection points of the table background grid, while the top row GCPs are laid on the 4 4’s. To obtain the 3D basin coordinates of surface tracers, the underdetermined system in Equa- tion 3.20 needs to be solved, using at least two cameras (stereo) looking at the same object. With an additional camera, Equation 3.20 becomes (Holland et al., 1997): 2 6 6 6 6 6 6 6 4 [L A 1 L A 2 u A i ] [L A 2 L A 10 u A ] [L A 3 L A 1 u A ] [L A 5 L A 9 v A ] [L A 6 L A 10 v A ] [L A 7 L A 11 v A ] [L B 1 L B 2 u B ] [L B 2 L B 10 u B ] [L B 3 L B 1 u B ] [L B 5 L B 9 v B ] [L B 6 L B 10 v B ] [L B 7 L B 11 v B ] 3 7 7 7 7 7 7 7 5 2 6 6 6 4 x y z 3 7 7 7 5 = 2 6 6 6 6 6 6 6 4 [u A L A 4 ] [v A L A 8 ] [u B L B 4 ] [v B L B 8 ] 3 7 7 7 7 7 7 7 5 ; (3.26) 43 which is an overdetermined system, and therefore makes the inverse problem (converting image- to world-coordinates) possible by utilizing a least square algorithm. Superscripts A and B denote the DLT coefficients from two different cameras, and image coordinates (from two different cameras) of a target with world coordinates x;y;z. However, before applying the coordinate transformation, the same object has to be identified in both cameras to obtain its two pairs of image coordinates. For the camera extrinsic parametrization part, the pairing was done manually, whereas for the pairing of the surface tracers in 3D PTV , the ray method (described in Section 3.5.3) is used. Particle identification and surface tracer matching is the subject of the following sections. The accuracy of the image-to-world-coordinate transformation for the stereo camera configu- ration can be assessed by comparing the predicted world coordinates of the GCPs (using Equa- tion 3.26) to the corresponding world coordinates measured with the total station. For the 3D case, the errore r is given by the mean distance between the real(x;y;z) and predicted(x p ;y p ;z p ) coordinates: e r = 1 N N å i=1 q [x p (i) x(i)] 2 +[y p (i) y(i)] 2 +[z p (i) z(i)] 2 ; (3.27) where N is the number of GCPs used for the extrinsic parametrization. The accuracy of the 3D coordinate transformation is given in Table 3.2, where e xy = 1 N å N i=1 p [x p (i) x(i)] 2 +[y p (i) y(i)] 2 and mean e z = 1 N å N i=1 jz p (i) z(i)j. Note that Table 3.2 shows the mean values for all 8 camera positions. Table 3.2: Accuracy of 3D coordinate transformation. e r (mm) meane xy (mm) maxe xy (mm) meane z (mm) maxe z (mm) 2.01 1.15 3.26 1.47 5.18 3.5.3 3D ray parametric equation This section provides an alternative representation of the image-to-world coordinate DLT transfor- mation given in Section 3.5, that is only used for the ray intersection matching algorithm described in Section 3.4.2. This formulation uses a free parametera which allows to directly relate the image 44 coordinates to rays between the camera pinhole and the tracer. The 3D ray parametric equation relates the image (u;v) to world (x;y;z) coordinates through (Faugeras and Luong, 2001; Doux- champs et al., 2002): a 0 B B B @ u v 1 1 C C C A = A 0 B B B @ x y z 1 C C C A + b; (3.28) where matrix A and vector b are obtained using a least-square algorithm and the known GCP coordinate pairs. Inversely, any object in the image with image coordinates(u;v) belongs to a 3D ray with parametric equation (Douxchamps et al., 2002): r(a)= p+aq; (3.29) wherea is a free parameter and vectors p and q are given by p=(A) 1 b q=(A) 1 0 B B B @ u v 1 1 C C C A : (3.30) 3.6 Surface tracer’s horizontal position correction With reference to Figure 3.13, a camera with known optical center coordinates (x c ;y c ;z c ) is observ- ing a surface tracer that initially rests on the still water level (z= 0:55cm for this experiment). The wavemaker motion results in an instantaneous positive change elevationDz of the water surface. When rectifying the still camera image with the assumption the water level is still at rest, the tracer horizontal coordinates will be displaced by a distance Dx= rDz z c h 0 (3.31) where r is the horizontal distance between the camera optical center and the surface tracer (at rest position) and h 0 is the still water depth. All parameters, except h 0 depend, of course, on time. 45 Figure 3.13: Geometrical depiction of the flat surface simplification: the horizontal coordinates of the surface tracer, that is elevated from the still water level, are erroneously calculated. Knowledge of its vertical displacement can reduce the error. The magnitude of the errorDx given in Equation 3.31 can be evaluated by looking at the range of values the independent parameters are taking. The horizontal distances covered by the camera FOVs in the 2D PTV setting are shown in Figure 3.14(a). With the exception of the bottom right part of camera C’s FOV (which is not of big importance for the experiment), maximum horizon- tal distances correspond to about 12-14m. Also, from the surface elevation information collected during the wave gauges stage of the experiments (see Section 3.13), it is known that surface eleva- tion fluctuation ranged between15cm to+3cm. In this range, the resulting errorDx is depicted in Figure 3.14(b). The maximum anticipated error ( 22cm) occurs when the surface tracer is located at the center of the eddy (sharp local depression), in the first stages of the experiment, and far away from the camera observing it. A more realistic value for the maximum error isDx= 18cm (for r = 11m, Dz = 15cm), and occurs when the clockwise rotating eddy forms in the inshore side (close to the tip of the breakwater), after the leading depression wave arrives. To reduce this error, the surface elevation measurements can be applied to Equation 3.31, so as to estimate the correction needed for the x;y coordinates of the surface tracers, obtained from Equation 3.21. It is important to stress here, that while the error of Equation 3.31 depends on the instantaneous surface elevation h, the error for the instantaneous velocity depends on the rate of change of the surface elevation ¶h=¶t. As long as the rate ¶h=¶t remains small, the surface velocity inferred from PTV is unaffected by the assumption of a flat water surface. However, the velocity vectors 46 horizontal distance from cameras (m) 5 10 15 20 x(m) 0 2 4 6 8 10 12 y(m) 0 5 10 15 horizontal error (cm) 0 5 10 horizontal distance from camera (m) -15 -10 -5 0 dz (cm) -20 -10 0 Figure 3.14: Left: horizontal distances between the overhead cameras and the water surface, cov- ered by the respective FOVs. Right: error in estimation of horizontal coordinates of tracers due to the flat surface assumption (see Equation 3.31). The camera height was set to z c = 10m (a mean value of the elevation of the four cameras) and mean water depth was set to h 0 = 55cm. will be misplaced, if the correction of Equation 3.31 is not applied and the spatial gradientsÑh (strongest near the TCS center) will affect the evaluation of metrics involving gradients, such as the vorticity. The application of Equation 3.31 can be tested for the overlapping FOV regions (see Fig- ure 2.3). The position r i of each tracer P i in the overlapping FOVs is calculated separately for each camera (r (A) i , r (B) i , r (C) i and r (D) i ). The mean distance mean jr (A) i r (B) i j (3.32) calculated for camera pairs with overlapping FOVs (Equation 3.32 camera labels are shown as an example for cameras A&B) provides an indication of how well the coordinate transfor- mation works, with and without the water level correction. The comparison is shown in Fig- ures 3.15&3.16. Evidently, the application of Equation 3.31 significantly reduces the horizontal positioning error, which remains almost constant, and equal to the error in the estimate of the unperturbed water surface. 47 0 20 40 60 80 100 120 140 160 180 200 0 5 10 15 mean error time (sec) error (cm) A−B A−D B−C B−D Figure 3.15: The mean error of tracer’s horizontal position in the overlapping camera FOVs without applying the water level correction. Legends correspond to the overlapping cameras FOV’s shown in Figure 2.3. 0 20 40 60 80 100 120 140 160 180 200 0 1 2 3 4 mean error time (sec) error (cm) A−B A−D B−C B−D Figure 3.16: The mean error of tracer’s horizontal position in the overlapping camera FOVs after applying the water level correction. Legends correspond to the overlapping cameras FOVs shown in Figure 2.3. 3.7 2D surface velocity extraction The data collected from the 2D PTV methodology described in the previous sections are used to extract velocity vectors at the water surface. Velocities in the horizontal directions are computed from successive frames, for each tracer P i , at time t j using the backward finite-difference scheme: u j i = x j i x j1 i t j t j1 v j i = y j i y j1 i t j t j1 ; (3.33) Even though the particle center detection has sub-pixel accuracy, the tracer velocity timeseries experienced a high-frequency fluctuation, subsequently removed using a low-pass filter with a 1.5Hz cutoff. Figure 3.17 shows the extracted velocity vectors over all camera FOVs for a single 48 time step. The velocity vectors from all cameras and for each camera frame were merged together. In the overlapping FOV regions, where particles were traced from more than one camera, the tracer velocity vectors and positions were averaged. Figure 3.17: The four camera FOVs merged together for a single time step. Tracer velocity vectors are overlaid. Each camera is assigned a different vector color. 3.8 3D surface velocity extraction For the 3D PTV experiments, the velocity vectors were also calculated using the backward finite difference scheme (Equation 3.33), but also computing the vertical velocity (w). The tracer vectors from different camera positions were merged together to visualize the flow around the tip of the breakwater (Figure 3.18), as is done for the 2D PTV experiments. However, the tracers from different camera positions correspond to different trials, and if the experiment is not perfectly repeatable, the velocities at the overlapping FOV regions would not necessarily match. For the same reason, the tracer velocity vectors in the overlapping FOV areas were not averaged. 49 2 m/s 2m 191 142 191 142 0 20 40 60 80 time (s) -0.1 0 0.1 η (m) 0 20 40 60 80 time (s) -0.1 0 0.1 η (m) Figure 3.18: Surface velocity vectors as obtained by the close-up stereo PTV analysis, 45 sec after the initiation of the experiment, showing how the flow separates at the tip of the breakwater and the the TCS is generated. Insets display a comparison between the surface elevation measured by two wave gauges (black curves) and the corresponding obtained by the stereo PTV analysis (red curves). The location of the wave gauges is shown with the red crosses. 3.9 Conversion to polar coordinates Studying the turbulent coherent structure (TCS) evolution requires the coordinate system to be formulated relative to the center of the TCS over time. This coordinate transformation also allows one to tie the velocity vectors from individual trials to a common spatial reference, since due to the turbulent nature of the flow, the TCS path is not common between trials (even though the same 50 experiment is repeated). Since a vortical structure is by nature cylindrical, the standard coordi- nate system employed is the cylindrical (or polar, in 2D). The equations to transform Cartesian coordinates(x;y), to polar(r;q) are given by: r= q (x x c ) 2 +(y y c ) 2 ; q = tan 1 y y c x x c ; (3.34) where(x c ;y c ) are the coordinates of the polar coordinate-system center, which in this application correspond to the TCS center. The process of identifying this TCS center is given in detail in Sec- tion 3.9.1.The velocity vectors referenced to the Cartesian coordinate system(u;v) are transformed to polar(u r ;u q ) using: u r = ucosq+ vsinq; u q = vcosq usinq: (3.35) 3.9.1 Identification of the TCS center Typically, vortices and their center are identified via vorticity, streamfunction, orl 2 -operator maps (Jeong and Hussain, 1995; Raffel et al., 2007; Seol and Jirka, 2010). The local extrema of these operators should define the center of flow rotation, provided the velocity field is well represented by the velocity maps. The computation of the operators requires discretization of the velocity field on a regular grid and applying a finite difference scheme. Whereas spatial discretization of the velocity field is inherent in PIV analysis, an interpolation scheme must be employed to regularize the PTV vector data. In this application, the remapping of the velocity field on a regular grid has proven to be a challenging task. In the early stages of the TCS development, the flow around the TCS center is characterized by zones of flow convergence and divergence, thus affecting the density of surface tracers. Even when the flow was well-seeded at the time of flow reversal trough the breakwater gap, the density of tracers was greatly reduced within 30 s after the TCS generation (Figure 3.19). Unless the flow was re-seeded during the evolution of the TCS, the tracers either conglomerated in the TCS center due to flow convergence, or were drawn away from the center due to flow 51 divergence. As a result, a ring formed around the TCS center in which no flow information could be extracted, thus affecting the data sampling for interpolation. In light of the above, the TCS center was identified using two different methods: 1. From vorticity maps, for as long as the density and spatial distribution of tracers allowed to create an accurate representation of the operator. 2. By tracking the center of mass of the conglomerated tracers at the TCS center. The procedure of identifying the TCS center is illustrated in Figure 3.21 and described in detail in the following sections. Figure 3.19: Tracer density around the TCS center at different times. Flow convergence/divergence either draws the tracers towards the TCS center or away from it. Vorticity-extracted TCS center To create votricity maps for each trial, the velocity fields (u;v) were interpolated on a grid with 20 cm cell size, using the natural-neighbor interpolation scheme (see Section 3.10.1). Lloyd et al. (1995) showed that this interpolation scheme was more accurate compared to other interpolation methods to regularize PTV vector data. From the interpolated velocity fields(u;v), the vorticity defined by w = ¶v ¶x ¶u ¶y (3.36) was computed, where the spatial derivatives were evaluated using the central-difference scheme (Equation 3.52a). Instead of defining the TCS center as the location of maximum vorticity w max , 52 the center of mass of the vorticity contour 0:7w max was found to produce more stable results. Seol and Jirka (2010) also applied this technique to infer the vortex center. Tracer conglomerate-extracted TCS center Once surface flow convergence/divergence significantly reduced the number of tracers in the TCS structure, the tracers accumulated at the core were used to determine the TCS center. The edges of the tracer-conglomerate were detected via image processing (detecting discontinuities in the brightness of the image), and the its center of mass defined the TCS center. In later stages of the TCS development, the associated flow field becomes universally weakly-divergent, with the tracer-conglomerate breaking-up into smaller fractions. At this point, instead of tracing the center of mass of a single, coherent tracer-conglomerate, the center of mass of all separated fractions (within a radius from the previous time step) was traced. Examples of this process are shown in Figure 3.20. Figure 3.20: Examples of detecting the edge of the tracer-conglomerate at the TCS center through image processing. At the early stages of TCS development, flow convergence keeps the tracer- conglomerate compact and in one piece. As time progresses, the flow becomes divergent, and the tracer-conglomerate expands and breaks up into smaller fractions. 3.9.2 TCS-center paths An example of a TCS-center path extracted from a single trial using the above methods is shown in Figure 3.22. The path is oscillating, in both the vorticity- and tracer conglomerate-extracted parts. In the case of the vorticity-extracted path, the high-frequency oscillations can partly be explained by the fact the the density and distribution of tracers in the TCS structure were not 53 sufficient to obtain accurate vorticity map representations. On the other hand, the lower-frequency oscillations of the tracer conglomerate-extracted path follow a pattern, which is the result of the tracer-conglomerate center spinning around the TCS-center and at the same time being translated with the speed of the TCS-center. The low-frequency motion pattern can be re-created for a perfectly azimuthal 2D flow field, such as in the case of a shielded-Gaussian vortex flow (see Equation 5.14b in Section 5.2). Two tracers are positioned at different radii (r 1 ;r 2 ) around the center of a shielded-Gaussian vortex with azimuthal velocity u q = 1 2 w 0 r exp((r=R) 2 )] (wherew 0 is the vorticity at the vortex center and R is the vortex core radius), with the vortex center being translated at a constant speed u. The result- ing motion shown in Figure 3.23 evidently follows the same pattern as that of the conglomerate- extracted TCS center path. The period of oscillation around the true path is a function of w 0 , whereas the amplitude of oscillation depends on the off-center radius. To obtain the true TCS-center paths, the raw paths were low-pass filtered using variable cutoff frequencies, appropriately decreasing as the TCS energy and TCS-center velocity decayed. The resulting TCS-center paths of all the individual trials, which provided the basis for the coordinate transformation of the PTV-extracted velocity vectors, are shown in Figure 3.24. 3.9.3 TCS-center propagation velocity The TCS-center velocity provides the means to represent the velocity field in the frame mov- ing with the TCS-center (Fl´ or and Eames, 2002) and is computed from the filtered TCS paths of the individual trials (Figure 3.25). The individual velocity timeseries appear noisy due to the inaccuracies in determining the TCS center through the vorticity maps of sparse velocity vectors. To increase the confidence level of the estimation, the mean TCS-center velocity time history was used instead of that of individual realizations. The standard error of the both the u and vvelocity expected values (defined as s= p n, where s is the standard deviation of the sample and n is the number of realizations) is less than 1 cm/s after 54 s (i.e. 95% confidence interval is2 cm/s of the mean). After 280 s, the TCS-center velocity components are very small and the velocity fields are assumed to be stationary. 54 Figure 3.21: Methodology for the extraction of the TCS center. In the early stages of the TCS development (up to point c1) - when the density of tracers around the TCS center is sufficient to interpolate the velocity field - the center is extracted from vorticity maps. In the next stage of TCS development, the edges of the blob are detected to compute the center of the TCS as the center of mass of the blob perimeter. The vortex path is subsequently filtered using different cutoff frequencies: 0.5s up to point f1, 0.2s up to point f2, and 0.015s until the end of the record. The filtered path is colored according to the corresponding vortex center speed. All images/subplots are to scale except the blob images showing the edge detection, which are scaled by a factor of two. 55 Figure 3.22: TCS center path from one trial. As time progresses, the TCS moves slower and the frequency of the oscillations decreases. The path was filtered using multiple low-pass filters of decreasing cutoff frequencies. 3.10 Lagrangian to Eulerian - gradients, interpolation, ensem- ble averaging, and azimuthal averaging Computing differential operators is more convenient when the velocity field is represented in an Eulerian framework. To re-map the scattered (Lagrangian) velocity vectors on an Eulerian grid, many methods have been suggested in the literature. One family of methods involves spatial inter- polation and the other space-averaging inside an interrogation window, similar to the PIV window averaging. 56 0 0.5 1 1.5 2 x (m) -0.2 -0.1 0 0.1 0.2 off-center fluctuation (m) r=12.9cm r=4.9cm 0 2 4 6 8 10 t (s) -0.2 -0.1 0 0.1 0.2 off-center fluctuation (m) r=12.9cm r=4.9cm Figure 3.23: Paths of off-center tracers in an idealized shielded-Gaussian vortex flow withw 0 = 10 s 1 , R= 1 m, and vortex-center translation velocity u= 0:25 m=s. Figure 3.24: TCS center paths from all trials. The mean path is shown with the bold curve. 3.10.1 Interpolation methods Adaptive Gaussian window The most-used interpolation method for PTV data is the Adaptive Gaussian Window (AGW) method (Ag¨ u´ ı and Jim´ enez, 1987). The Eulerian velocity field u(x;t) is computed by assigning distance-dependent weightsa i to measured (scattered) velocity values v i : u(x)= å i a i v i å i a i ; a i = exp jx x i j 2 H 2 ; (3.37) 57 −0.6 −0.4 −0.2 0.0 0.2 u (m/s) 50 100 150 200 250 time (s) −0.4 −0.2 0.0 0.2 0.4 v (m/s) 50 100 150 200 250 time (s) Figure 3.25: TCS center velocity realizations (u;v) in x- and y-directions. The mean is shown with the bold curves. where H is the width of the Gaussian function. If the mean inter-particle distance between N data points in an area A isd = p A=N, then Ag¨ u´ ı and Jim´ enez (1987) recommend H=d = 1:24, so as to minimize the interpolation error. The smallest scale that can be resolved based on the inter-particle spacing is equal to the Nyquist wave length l= 2d (St¨ uer and Blaser, 2000). The natural-neighbor method The interpolation method used here is the natural-neighbor algorithm which utilizes the V orono¨ ı tessellation of a set of discrete points in space, as described in 3.4.1. The estimate G(x;y) of a function f at point(x;y) in space is given by: G(x;y)= n å i=1 w i f(x i ;y i ); (3.38) where w i are weights applied on neighboring points(x i ;y i ). The weights are computed using two sets of V orono¨ ı tessellations: one with the point(x;y) included in the set and one without it. The weights are a function of the area of the individual polygons formed from the intersection of the V orono¨ ı cells of the neighbors with the cell corresponding to(x;y), when(x;y) is included in the points set. 58 Radial basis functions Although simple to implement, AGW and the natural-neighbor interpolation methods have been claimed as sensitive to noise and as providing erroneous gradient results in patches of scarse data and at the domain edges (Paul, 2012). Instead, the use of more sophisticated (and more computa- tionally demanding) interpolation schemes, such as Kriging, global basis functions, and Gaussian radial basis functions, have all been suggested to improve the predictions (e.g. Spedding and Rig- not (1993); Paul (2012); Casa and Krueger (2013)). Radial basis function (RBF) interpolation for the jth particle is computed using: u j (x)= N å i=1 b ji exp " jxm ji j 2 2H 2 ji # ; (3.39) where N is the number of RBFs, b ji are the RBF weights, m ji are the locations of RBF centers, and H ji are the widths of the Gaussian RBFs controlling the radial extend of the RBF (Casa and Krueger, 2013). The RBF weightsm ji are determined using a least-square optimization algorithm, u(x), and then its derivatives can be computed in user-defined locations in the flow. 3.10.2 Space-averaging Space-averaging of velocity vectors is performed to obtain the instantaneous mean velocity field. The Eulerian grid nodes P are assigned a regular grid spacingDx=Dy and an interrogation window of radius W R is applied at each grid node P i to locate velocity vectors lying inside the window. The sample size in each interrogation window is denoted here as N. In the Reynolds decomposition, the velocity u is decomposed in mean U and fluctuating u 0 components as u= U+ u 0 : (3.40) The mean velocity is derived by time-averaging using an appropriate time interval T as U(x;t)= 1 T Z t+T t u(x;s)ds; (3.41) 59 which is equivalent to the space-averaged mean velocity using a filter length scale L U(x;t)=hui(x;t)= 1 L Z L 0 u(x;t)dx; (3.42) (here hi denotes space-averaging), provided that turbulence is stationary and homogeneous (ergodic hypothesis). Since the scope of this work is not to quantify turbulence, but rather to obtain the mean flow field, it is sufficient to assume that turbulence is stationary and homogeneous, and then obtain the mean flow velocities using spatial-averaging. To obtain the mean value of the local velocity field in the interrogation window at a certain time, the top-hat filter is applied (Discetti et al., 2015). All velocity vectors inside the interrogation window are assigned equal weight to compute the mean velocity components(U(P i );V(P i )), with- out taking into account the spatial variability of the velocity field. The resulting mean velocities and fluctuations are a function of the chosen filter length. To extract the mean velocity field, an ensemble flow field was created from the instantaneous surface velocity vectors of the individual trials. In the initial stage, during TCS generation, the ensemble frame of reference is stationary (Figure 3.26), and the 3D PTV velocity vectors are also incorporated in the ensemble. For the time steps for which the TCS center known through the vorticity maps or conglomerate centers of the individual trials, the ensemble is TCS-centered (Figure 3.27), in which only the 2D PTV data are considered. The former ensemble is used to study the TCS generation and growth at the early stages, whereas the latter is used to study the TCS properties, transition of the flow to Q2D and energy decay due to bottom friction. Stationary ensemble For the stationary ensemble, the chosen interrogation window radius is W R = 0:35 m and the grid spacing isDx=Dy= W R , and results into an overlap of 39% between successive windows. In other words, the overlap area is 39% of each individual interrogation area. This combination was 60 + +...= trial 1 trial 2 ensemble Figure 3.26: Stationary ensemble creation from the surface velocity vectors of the individual trials. This ensemble is used to extract the mean flow velocities during TCS generation. found to provide adequate flow resolution and sample size per interrogation window. The mean spacing of tracers in the stationary ensemble domain with area A is approximated using d = s A N grid ; (3.43) and ranges between 12 and 15 cm, leading to a mean number of particles per grid areaDx 2 ranging 5 9, and per interrogation window ranging 17 27. This sample size satisfies the Nyquist sampling criterion (St¨ uer and Blaser, 2000), which requires that the grid node spacingDx is at least two times the sample spacingd, or Dx 2d: (3.44) 61 Figure 3.27: TCS-centered ensemble creation from the surface velocity vectors of the individual trials. The Nyquist criterion translates to requiring at least N grid 4 samples per grid area and represents a minimum since the tracers are not equally distributed in space. Velocity vectors in areas with sparse tracer distribution (N < 4) are obtained using the natural-neighbor interpolation scheme, discussed in Section 3.10.1. The shortest wavelength of a flow feature that can be resolved using this grid spacing is 2Dx= 0:7m. A sample mean velocity field is shown in Figure 3.28. 16 18 20 22 24 x(m) 6 7 8 9 10 11 12 13 y(m) ensemble velocity vectors 16 18 20 22 24 x(m) 6 7 8 9 10 11 12 13 y(m) mean velocity field Figure 3.28: Ensemble surface velocity vectors and resulting mean surface velocity field near the breakwater tip at t = 27s, using W R = 0:35 m andDx=Dy= W R . Frame shown is after the arrival of the leading depression wave. 62 The accuracy of the estimation of the mean depends on the number of samples available in the interrogation window and on the fluctuations of the velocity field. The standard error SE U in the mean velocity U is given by (Discetti et al., 2015) SE U = s ¯ u 02 +s 2 rand N ; (3.45) where ¯ u 02 is the root mean squared fluctuation ands rand is the random error due to the positioning of the surface vectors. ¯ u 02 can be defined with respect to the 2D turbulence intensity T:I:= s 0:5( ¯ u 02 + ¯ v 02 ) U 2 +V 2 ; (3.46) where( ¯ u 0 ; ¯ v 0 ) and(U;V) are the root-mean square velocity fluctuations and mean velocities in the x and y directions respectively. Another method of estimating the error of the mean is using the Jacknife technique (St¨ uer and Blaser, 2000). The estimate of the standard error in U is given by SE JU = s N 1 N N å i=1 (ˆ u i U) 2 ; (3.47) where N is the sample size and ˆ u i is the mean velocity computed by leaving out the i th observation. For an unbiased sample and ¯ u 02 >>s rand , equations 3.45& 3.47 should produce the same result. The mean and maximum standard errors (using Eqn. 3.47) for the stationary ensemble are shown in Figure 3.29. The metrics are evaluated for the time period during TCS generation (t = 1073s). There is a higher density of samples around the channel due to the integration of the 3D PTV vectors in the ensemble. Most of the tracers were initially placed near the channel to capture the TCS surface velocities, which explains the low density of tracers in the domain boundaries. Mean and maximum velocity standard errors of the order of 5 and 20 cm=s are observed at the separation zone of the breakwater tip and near the top wall respectively, where tracer density is inherently low and turbulence intensity is highest. The mean ratio of the standard error to the velocity magnitude is 10% around the breakwater tip, but remains 6% in most of the domain. 63 6 8 10 12 y (m) mean no. of tracers per window 0 10 20 30 40 50 mean(SE/velocity_magnitude) (%) 0 2 4 6 8 10 6 8 10 12 y (m) mean SE−U (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 mean SE−V (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 6 8 10 12 y (m) 4 6 8 10 12 14 16 18 20 22 24 x (m) max SE−U (m/s) 0.00 0.04 0.08 0.12 0.16 0.20 4 6 8 10 12 14 16 18 20 22 24 x (m) max SE−V (m/s) 0.00 0.04 0.08 0.12 0.16 0.20 Figure 3.29: Number of tracers per interrogation window and standard error of the mean velocities (Eqn. 3.47) for the stationary ensemble. Metrics evaluated in the time interval t=10-73 s using W R = 0:35 m andDx=Dy= W R . TCS-centered ensemble Since the main purpose of the TCS-centered ensemble is to study the TCS evolution, its domain extent is limited in radius to the closest vertical boundary of the basin (d min , see Figure 3.30). This distance represents the maximum possible radius the TCS may attain. At any given time step (i), d i min is defined as the global minimum of the minimum distances (d i j ) between the TCS centers (X c i j ) of the individual trials ( j= 1:::N trials ) and the closest vertical boundary (¶B) as d i min = min j=1::N trials d i j = min j=1::N trials kX c i j ¶Bk min ; (3.48) where the boundary ¶B corresponds to the offshore basin perimeter, defined by the side walls, breakwater and retracted wavemaker. 64 0 500 1000 1500 2000 2500 3000 time (s) 0 1 2 3 4 5 6 7 8 d min from boundaries (m) Figure 3.30: Minimum distance of the TCS-center to a vertical boundary. The mean of all trials is shown with the solid line, whereas dashed lines correspond to the min/max envelope for the individual trials. Similar to the stationary-ensemble, the TCS-centered ensemble domain, with spatial extent x;y2[d i min ;d i min ], was discretized into nodes of regular spacing Dx;Dy. The grid was used to obtain the mean velocity field and evaluate the metrics that require regular grid spacing, such as vorticity and swirl strength. The collected surface velocity information was at times limited by the spatial coverage of the camera FOVs (see Figure 2.3). As a consequence, part of the grid domain was cut-off at the FOV boundaries. The number of contributing trials and mean tracer spacing in the evaluation domain are plotted in Figure 3.31. The mean tracer spacing ranges from 14cm at t = 100s - when the number of contributing trials was high - to 40cm at t= 3000s. Only four trials extended to this point. Taking the above into account, the choice of spacing for the TCS-centered ensemble is Dx=Dy= 0:4m with an interrogation window of radius W R = 0:4m. The uncertainty of the mean velocity field for the TCS-center ensemble is estimated from the velocity fluctuations at the TCS- center. Making the realistic assumption that the center represents the point with the highest turbu- lence intensity, this provides an upper limit of uncertainty. The standard deviation of the sample corresponding to the node at the TCS-center is shown in Figure 3.32a for all times considered in the TCS-centered ensemble. The maximum standard error can be computed in combination with the mean spacingd shown in Figure 3.31b: a mean number of samples per interrogation window 65 0 1000 2000 3000 0 5 10 15 20 25 time (s) # of trials (a) 0 1000 2000 3000 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 time (s) mean tracer spacing (m) (b) Figure 3.31: Statistics of the TCS-centered ensemble for the duration of the experiment. (a) The number of trials. (b) Mean tracer spacing in the domain extending to x;y2[d i min ;d i min ]. N mean is computed using Eqn. 3.43 and the maximum standard error is computed using Eqn. 3.45, assuming ¯ u 02 >> s rand . The resulting maximum standard error is shown in Figure 3.32b. As with the stationary ensemble, velocity vectors in nodes with sparse tracer distribution (N < 4), are obtained using the natural-neighbor interpolation scheme. 0 1000 2000 3000 time (s) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 std (m/s) (a) u r u θ 0 1000 2000 3000 time (s) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 SE (m/s) (b) u r u θ Figure 3.32: Uncertainty of the mean velocity field for the TCS-centered ensemble. (a) Standard deviation of the u r and u q velocity components, evaluated at the TCS-center. (b) Standard error of the mean. 66 Azimuthal-averaging Since the laboratory TCS under study is not axisymmetric, it is practical to examine its properties on an azimuthal-averaged profile. An easy way to define the evaluation points for the profile is by using constant dr and dq spacing (Figure 3.33a). This method will keep the number of evaluation points per radius constant to obtain the mean value, with the drawback that the spacing of the points becomes a function or r (i.e. the larger the radius the further apart the nodes are along the circle). Then the evaluation points cease to contain equally-spaced information on the larger circle perimeters. By contrast, the evaluation points can have constant spacing and a variable number of points per radius r (Figure 3.33b) by setting dq equal to dq(r)= 2arcsin dr 2r : (3.49) If a metric is evaluated on the sample points using space-averaging, the spacing dr can be defined in terms of the interrogation window radius W R . Here, the equally-spaced nodes are not used to obtain an average, but rather to interpolate vorticity values evaluated on a regularly-spaced grid (as described in the previous section). The interpolated vorticity values are then averaged in the azimuthal direction. d θ d r (a) (b) (c) Figure 3.33: Azimuthal-average sampling points for (a) variable spacing and constant dq, and (b) for constant spacing and variable dq; (c) interrogation windows in the form of annuli. The dashed red lines (b) correspond to interrogation windows of radius W R = dr and blue dashed lines (c) correspond to the annuli used as interrogation windows. Another method that can practically be employed to obtain azimuthal-averaged velocity pro- files, is to employ annuli as interrogation windows (Figure 3.33c). This method doesn’t involve 67 interpolation, only averaging of all velocity vectors located inside each annulus. The annuli spac- ing (Dr) is dictated by the average spacing of the vectors in the radial direction (d r ), satisfying the Nyquist criterion (Dr 2d r ), and the annuli are overlapping by 50% to achieve higher resolution. The chosen radial spacing that produced satisfactory results isDr= 0:25m coupled with an annu- lus width of W = 0:5m. An application example using this azimuth-averaging method is shown in Figure 3.34. -5 0 5 distance from center (m) -5 0 5 distance from center (m) 0 2 4 r (m) 0 0.1 0.2 0.3 u θ (m/s) 0 2 4 r (m) -0.1 -0.05 0 0.05 0.1 u r (m/s) Figure 3.34: Azimuthal-averaging of the TCS azimuthal (u q ) and radial (u r ) velocity components - example shown for time t 333s. Black vectors and dots correspond to the scattered velocities and red curves to the azimuthal-averaged profiles. Error bars indicate the standard deviation of the sample velocities in each annulus of width W = 0:5m. 3.10.3 Spatial derivatives Local gradients of any function f are computed using the four point, second-order accurate, least square differential operator (Raffel et al., 2007) given by d f dx i 2 f i+2 + f i+1 f i1 2 f i2 10Dx + O(Dx 2 ): (3.50) This scheme has been designed to reduce the effect of the measurement uncertainty and is recom- mended to be used to produce vorticity maps using PIV data (Raffel et al., 2007). The drawback is that it has the tendency to over-smooth the differential due to the higher weights of the sampling points at 2Dx which has a direct effect on the calculation of maximum vorticity w max at the TCS center. Due to the steep vorticity gradient near the TCS center, w max computed on a grid with 40 68 cm cell size is already expected to be much lower than the actual value. The tracer-configurations provide more detailed vorticity measurements near the center, as discussed in next section. The total vorticity, evaluated over an area around the TCS center (e.g. circulation) using dif- ferent differential operators, was found to be less sensitive on the method used. The above is illustrated for a specific time instance in Figure 3.35. The comparison of the least square operator with the third-order accurate Richardson scheme (Eqn. 3.51, Raffel et al. (2007)) shows a large deviation for the maximum vorticity value, but both operators produce approximately the same result for the circulation over an area around the vortex center. d f dx i f i2 8 f i1 + 8 f i+1 f i+2 12Dx + O(Dx 3 ): (3.51) Figure 3.35: V orticity maps created using the least-squares and Richardson differential schemes at time t = 166s. Maximum vorticity (w max ) and circulation (G, evaluated inside the green circle) values are noted. The mean velocity field is evaluated on a grid withDx=Dy= 0:4m. Around the boundaries, where the stencil is limited, the central, forward and backward differ- ence schemes are used d f dx i f i+1 f i1 2Dx + O(Dx 2 ); (3.52a) d f dx i f i+1 f i Dx + O(Dx); (3.52b) d f dx i f i f i1 Dx + O(Dx): (3.52c) 69 3.11 Vorticity measurements To extract the vorticity near the TCS-center, surface–tracer configurations were used that allowed to compute the local flow rotation. The two configurations that were tested are shown in Fig- ure 3.36. Both consist of four tracers that constitute the vertices of a square. In the first, the sides are interconnected (hereafter called square-tracer), whereas in the second the diagonals were inter- connected to form a cross (hereafter called cross-tracer). The sides of the two configurations are equal to 18 and 7 cm, with corresponding diagonals of 25:5 and 9:9 cm respectively. r=12.7cm (a) dθ i (c) tracer k vertex i, i=1:4 time j time j-1 r=4.9cm (b) Figure 3.36: (a), (b) The two tracer configurations used to infer the vortex–center vorticity. The left configuration uses a 12.7 cm radius, whereas the right configuration uses a 4.9 cm radius; (c) Schematic of measuring the angular velocity of a square-tracer. The square-tracer was used in one trial, with only one such tracer being introduced in the flow. The square-tracer merged to the center of the TCS, with its vertices spinning around it. Multiple cross-tracers were used in two trials, with one tracer merging to the TCS center in each trial. The vertices of the tracer configurations were tracked using the 2D PTV method described in Section 3.3. For the cross-tracer trials, where more than one tracer configurations were used, the cross-tracer perimeter was used as a template to identify the individual configurations. The center of mass of each configuration was used to convert the the x-y basin coordinates to polar(r;q). The angular velocity (W, measured in rad/time; counterclockwise rotation is positive) of each vertex i, at time j was computed using the backward finite difference scheme (see Figure 3.36c): W j i = q j i q j1 i t j t j1 = dq j i t j t j1 : (3.53) 70 The instantaneous vorticity of the tracer configuration k was computed as 2 times the mean angular velocity of the four vertices w j k = 2hW j i=1:4 i; (3.54) since u q = rW, and the vorticity assuming solid-body rotation is given byw z (r)= 1 r (u q + r ¶u q ¶r )= 1 r (rW+ r ¶(rW) ¶r )= 2W. The paths of the tracers which followed the TCS center, and their respective vorticity time- series are shown in Figures 3.37&3.38a. The paths indicate that the center of mass of the tracers was rotating around the TCS-center path and therefore the tracer angular speed measured doesn’t exactly correspond to the TCS center. The fluctuation in the vorticity timeseries is a result of the off-center rotation and the translating motion of the tracers. Since it wasn’t feasible to reference the vorticity measurements to the TCS center, the timeseries were low-pass filtered to obtain a smooth w max decay curve (Figure 3.38b). The resulting decay curves for the two cross-tracer trials match well, justifying the repeatability of the experiment with the tracer configurations. The vorticity offset in the decay curves for the two different configurations is due to the difference in the radial distance to the configuration tracers. 3.12 ADV data The ADV data are provided in raw format. The Vectrino measures four velocity components: u, v, w 1 and w 2 , where w 1 and w 2 are independent and redundant measurement of the vertical velocity (useful for turbulence measurements). The output file of each ADV contains a velocity, signal amplitude and correlation column for each of the four beams. From the background noise and amplitude readings, the signal to noise ratio (SNR) can be computed as SNR= 20log 10 amplitude signal amplitude noise ; (3.55) measured in dB. The correlation is a measure of self-similarity, affected by the strength of the acoustic echo and the transmitter signal interference. The manufacturer recommends SNR 10dB 71 Figure 3.37: The paths of the tracer configurations following the TCS center, overlaid on rectified images of the trial in which a square-tracer was deployed. 0 500 1000 time (s) 0 1 2 3 4 5 6 7 8 ω max (s -1 ) (a) square-tracer cross-tracer cross-tracer 0 500 1000 time (s) 0 1 2 3 4 5 6 7 8 ω max (s -1 ) (b) Figure 3.38: Raw (a) and filtered (b) vorticity decay curves near the TCS center, as measured using the square- and cross-tracers. Colors match the paths shown in Figure 3.37. and correlation> 70%, for a relatively clean signal. Values below the recommended imply that the medium was not properly seeded and that there was not enough scattering. Out of the five 72 trials in which all ADVs were mounted, only the third trial resulted in a relatively clean signal (Figure 3.39, Table 3.3). The signal was de-spiked using the algorithm of Mori (2014), which is based on the work of Mori et al. (2007), and then a low-pass filter was applied. −0.50 −0.25 0.00 0.25 0.50 (m/s) 0 50 100 150 200 250 ADV1 0 50 100 150 200 250 ADV1 u v w1 −0.50 −0.25 0.00 0.25 0.50 (m/s) 0 50 100 150 200 250 ADV2 0 50 100 150 200 250 ADV2 −1.0 −0.5 0.0 0.5 1.0 (m/s) 0 50 100 150 200 250 time (s) ADV3 0 50 100 150 200 250 time (s) ADV3 Figure 3.39: ADV raw (left column) and processed data (right column). Table 3.3: Mean SNR and correlation values for all ADV beams of trial 3. ADV 1 ADV 2 ADV 3 u v w 1 w 2 u v w 1 w 2 u v w 1 w 2 correlation(%) 84 88 82 89 82 85 80 82 85 84 82 86 SNR(dB) 9.0 7.8 6.8 5.1 8.3 6.4 8.3 6.4 10.9 8.4 8.4 6.4 The ADV data was used to validate the 2D PTV methodology. Figure 3.40 compares the hori- zontal velocities from ADVs 1&3 with the surface velocities extracted from the 2D PTV analysis. The 2D PTV curves are the result of interpolating the surface velocity vectors using the natural neighbor interpolation algorithm. The comparison is excellent for both the u and v velocity com- ponents. The small differences observed for ADV3 can be partly attributed to its location, which didn’t allow for accurate interpolation of the scarce PTV velocity vectors around this ADV probe. 73 Discrepancies for ADV1 at later times (> 200 s) can be attributed to the fact that the flow field becomes highly dependent on the TCS path of the particular trial the ADV measurements were taken. −0.50 −0.25 0.00 0.25 0.50 u (m/s) 0 50 100 150 200 250 ADV1 −0.50 −0.25 0.00 0.25 0.50 v (m/s) 0 50 100 150 200 250 ADV1 ADV PTV −1.0 −0.5 0.0 0.5 1.0 u (m/s) 0 50 100 150 200 250 time (s) ADV3 −0.50 −0.25 0.00 0.25 0.50 v (m/s) 0 50 100 150 200 250 time (s) ADV3 Figure 3.40: Comparison between horizontal velocities extracted from ADV 1 (located at x = 4:52m, y= 0:00m), ADV 3 (located at x= 20:69m, y= 13:18m) and 2D PTV . 3.13 Free surface elevation Wave gauges were placed in a grid layout with variable grid spacing (see Figure 2.7). For each wave gauge, the voltage output was converted to surface elevation. The collected dataset is used to obtain a series of continuous surface elevation maps across the basin plan area and visualize the long wave motion forced by the wavemaker. A Delaunay triangulation mesh was created from the gauge nodes (Figure 3.41), from which the surface elevation values at a regular node mesh was interpolated. The breakwater constitutes a hard boundary over which interpolation is not allowed. Due to the constraint imposed by the breakwater, the triangulation mesh excludes wet areas around the breakwater perimeter. For the areas outside the triangulation mesh, an extrapolation algorithm needed to be applied. The interpolated surface elevation maps, apart from being used for visualization, are also used as input into the surface tracer horizontal position correction described in Section 3.6. For this 74 0 5 10 15 20 25 30 35 40 45 −15 −10 −5 0 5 10 15 wavemaker Figure 3.41: Delaunay triangulation mesh (blue lines) to interpolate the wave gauge surface eleva- tion data (red dots). The breakwater perimeter is set as a hard boundary. purpose, a separate interpolation mesh of 10 cm cell size was created, that covers the study area. To avoid leaving out wet areas around the breakwater perimeter, each interpolation mesh point is assigned a separate set of sampling-points. The sampling-point set of each interpolation mesh point corresponds to all gauges the interpolation point can be linearly connected to, without going over a hard boundary. The bi-harmonic spline interpolation/extrapolation scheme of Sandwell (1987) is then applied in Matlab. Although much more computationally demanding, the resulting maps provide continuous surface elevation data for all wet points of the study area. An example of an interpolated surface elevation map for the study area is shown in Figure 3.42. The wave gauge data can also be used to validate the results from the 3D PTV analysis. Fig- ures3.43-3.44 compare the zcoordinate extracted from the surface tracers to the surface eleva- tion timeseries at selected wave gauges around the breakwater tip. The PTV curves represent the interpolated mean surface-elevation of all trials. The numbering of the gauges is provided in Fig- ure 3.45. The comparison shows that, in general, the 3D PTV analysis assigns accurate elevations to the surface tracers. The largest miss-match can be observed during the time periods when the 75 Figure 3.42: Interpolated surface elevation map in the study area 5 s after wavemaker start. The inset shows the time relative to the wameker motion. Wave gauge locations shown with black dots. TCS center was in the vicinity of the the wave gauge. This time period is expressed in the time- series as a local minimum with small-amplitude fluctuation in both the PTV- and gauge-extracted excitations. This mismatch is attributed to the chaotic–turbulent nature of the flow during that time. The center of the inshore vortex followed different trajectories in every trial, and therefore the surface elevation distribution was also varying between trials. 76 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 95; x 22.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 111; x 21.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 143; x 19.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 159; x 18.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 175; x 17.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 191; x 16.0 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 431; x 16.3 (m) y 11.5 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 174; x 17.0 (m) y 10.2 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 190; x 16.0 (m) y 10.2 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 427; x 16.3 (m) y 10.2 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 142; x 19.0 (m) y 10.2 (m) 0 20 40 60 80 η (m) -0.1 0 0.1 wave gauge: 158; x 18.0 (m) y 10.2 (m) time (s) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 94; x 22.0 (m) y 10.2 (m) time (s) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 110; x 21.0 (m) y 10.2 (m) Figure 3.43: Comparison of surface elevation extracted from the surface tracers (3D PTV , red dashed curves) and from the wave gauges (black curves). 77 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 93; x 22.0 (m) y 9.0 (m) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 109; x 21.0 (m) y 9.0 (m) time (s) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 92; x 22.0 (m) y 7.7 (m) time (s) 0 20 40 60 80 η (m) -0.05 0 0.05 wave gauge: 108; x 21.0 (m) y 7.7 (m) Figure 3.44: Comparison of surface elevation extracted from the surface tracers (3D PTV) and from the wave gauges (continued from Figure 3.43). x (m) 15 16 17 18 19 20 21 22 23 y (m) 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 12.5 95 111 143 159 175 191 431 174 190 427 142 158 94 110 93 109 92 108 Figure 3.45: Wave gauge numbering corresponding to Figures 3.43-3.44. 78 Chapter 4 TCS generation and circulation growth This chapter examines the hydrodynamic processes that lead to the generation and short-term evo- lution of the turbulent coherent structure. The analysis on vortex circulation growth is based on the theory developed on vortex ring formation (e.g. (Didden, 1979)), which is also applicable for vortex generation due to topographic forcing in geophysical scales (e.g. (Wells and van Heijst, 2003)). 4.1 Theoretical background on circulation growth Figure 4.1: Schematic of (a) nozzle and (b) orifice type piston-cylinder vortex ring generators. Source: Krueger (2005) 79 For zero-swirl flows, it can be postulated that the vortex circulation growth is a result of the vorticity influx u x w q through the exit plane (Krueger, 2005). The rate of change of the circulation for an axisymmetric flow through a cylinder of radius R 0 can thus be represented as (Krieg and Mohseni, 2013) dG dt = Z R 0 0 u x w q dr= Z R 0 0 u x ¶u r ¶x ¶u x ¶r dr: (4.1) By defining the radius to the edge of the boundary layer forming near the inner cylinder wall as r= y 1 and the corresponding velocity as u x (y 1 )= u 1 , the integral in Eqn. 4.1 can be separated into parts as (Didden, 1979) dG dt = Z y 1 0 u x w q dr+ Z R 0 y 1 u x ¶u r ¶x dr Z R 0 y 1 u x ¶u x ¶r dr= I 1 + I 2 + 1 2 u 2 1 : (4.2) This definition allows one to evaluate the circulation influx without the need of having detailed flow information in the boundary layer, where it is difficult to collect measurements. Further assuming that the radial velocity component is negligible (u r 0 and thus w q ¶u x =¶r), and that the velocity u x is constant for r< y 1 leads to what is commonly referred to as the slug model equation: dG sm dt = 1 2 u 2 1 ; (4.3) where u cl is the velocity at the center-line of the exit cylinder plane. Finally, assuming that the flow in the cylinder is uniform (outside the boundary layer) and equal to the piston velocity u cl u p (t), the circulation time history can be calculated as G sm (t)= 1 2 Z t 0 u p (t 0 ) 2 dt 0 : (4.4) Apart from its application in vortex ring generation, the slug model is also commonly used to approximate the circulation growth of dipoles generated at coastal inlets due to tidal flushing (e.g. Wells and van Heijst (2003)). The methodology of Krueger (2005) to determine the circulation build-up in the domain by a piston-cylinder mechanism is adapted and adjusted to the geometry of this experiment. The rate of change of the circulation on an arbitrary plane (evaluated at certain depth in the water column), 80 with bounding contour C, can be derived from the vorticity transport equation for incompressible flow as ¶G ¶t = I C (w u+nÑw) dc: (4.5) The termw u represents vorticity influx, whereas the termnÑw represents a vorticity source (Krueger, 2005). Eqn 4.5 is evaluated on the planar contour shown in Figure 4.2, which corresponds to the perimeter of the offshore basin. All the terms w;¶=¶z;w x ;w y are ignored assuming the flow is 2D. C 1 2 3 4 5 (x 2 ,y 2 ) (x 1 ,y 1 ) (x 1 ,y 2 ) x y u v wavemaker Figure 4.2: Bounding contour C to determine the circulation influx to the offshore basin during 1 st flow reversal. Evaluating Eqn 4.5 along segment 1 produces ¶G 1 ¶t = Z y 2 y 1 uw z n ¶w z ¶x x=x 1 dy= Z y 2 y 1 u ¶v ¶x n ¶w z ¶x x=x 1 dy; (4.6) where the second integral is the result of R y 2 y 1 u(¶u=¶y)dy= 0, due to no-slip boundary conditions at y= y 1 ;y 2 . Thus, the total rate of change of circulation is given by ¶G ¶t = ¶G 1 ¶t + 5 å i=2 ¶G i ¶t = Z y 2 y 1 u ¶v ¶x x=x 1 dy+n I C (Ñw) dc: (4.7) 81 The dominant first term, corresponding to the vorticity influx through the channel, is only a func- tion of the horizontal gradient of the v velocity. It can be visualized as the curvature of the stream- lines at the channel. The second term, corresponding to the vorticity generation/diffusion in the perimeter of the planar contour, is expected to be small and therefore can be typically omitted (e.g. Wells and van Heijst (2003),Bryant et al. (2012)), leading to ¶G ¶t Z y 2 y 1 u ¶v ¶x x=x 1 dy: (4.8) Comparing Eqn. 4.3 with Eqn. 4.8, it is immediately apparent that the widely used slug model does not apply to this experimental setup. In the vortex ring-generation and tidal-starting jet exper- iments, circulation production is primarily the result of the influx of vorticity generated at the cylinder/orifice boundary layers. In this setup, the resulting circulation is again a result of vorticity influx through the exit plane (channel), but the domain lacks symmetry to make any further sim- plifications. The same methodology for calculating circulation production through vorticity influx through the channel is also applicable to the inshore basin. 4.2 Circulation growth theoretical framework for the experi- mental layout The growth of the TCS circulation is studied for the first cycle, during which coherent structures are generated on either side of the breakwater (see Figure 4.3). The generation of the TCS in the inshore basin is interrupted by the channel flow reversal, forming the main TCS on the off- shore side. It is of interest to analyze both flow phases, since the flow profiles across the channel, generating the coherent structures, have different characteristics. In the circulation-growth analysis that follows, in order to make the results comparable to other studies using different experimental setups, data are non-dimensionalized using scale variables. In their study of tidal starting jet vortices at barotropic inlets, Nicolau del Roure et al. (2009) used the peak channel-flow velocity and inlet width to scale time. In this experiment, these scale variables correspond to peak channel-flow velocity ¯ U max 0:67m=s and channel width W = 2 82 0 10 20 30 40 50 60 70 80 90 100 −0.5 0 0.5 inshore TCS generation offshore TCS generation t=9s t=40.6s t=78s time (s) channel velocity (m/s) channel−average velocity during 1st cycle Figure 4.3: Channel-average velocity evaluated using the offshore arc of length 3:2m shown in Figure 4.4 as U = Q=3:2, where Q is the flow rate per unit depth: Q= R l 0 ru r dq. 3:2= 6:4m - the width is multiplied by a factor of 2 to match the geometry of a symmetric inlet, as in experiments of Nicolau del Roure et al. (2009). Since the channel forcing signal is not symmetrical in this experiment, two different peak channel-flow velocities are considered, one for the inshore (1) and one for the offshore TCS generation phases (2): ¯ U max1 = 0:43m and ¯ U max2 = 0:67m respectively. Hereafter, for either flow phase, the peak channel-flow velocity is noted as U for simplicity. Using this convention, dimensionless time (designated a notation) is defined as t=(tt 0 )U=W; (4.9) where t 0 corresponds to the time at which circulation starts being produced in each flow phase: t 0 9s for the inshore TCS generation, and t 0 40:6s for the offshore TCS generation (see Figure 4.3). An alternative time scaling can be defined through an “effective” piston stroke L, related to the channel velocity u c as L= R t 0 u c (t 0 )dt 0 (Gharib et al., 1998). If the channel flow was idealized as being driven through a piston, L would represent the equivalent piston stroke. The dimensionless timet proposed by Gharib et al. (1998) to scale time for vortex ring formation experimental data, referred to as the “formation time”, is given by t = L=W; (4.10) where again W = 2 3:2= 6:4m is twice the channel width. In the definition of Gharib et al. (1998), W represents the cylinder diameter. 83 4.2.1 Local coordinate system To estimate the vorticity influx in the inshore and offshore domains, the channel velocities and velocity gradients are defined along arcs connecting the side wall with the breakwater tip (see Fig- ure 4.4). The center (X c = 14:39m;Y c = 13:25m) of the circle used to define the arc corresponds to the intersection point of the the top wall and the extension of the breakwater tip face (the circle has the breakwater as tangent and is centered along the top boundary). V orticity is never created or destroyed in the interior of the fluid, but only at the domain boundaries (Batchelor, 1967). Assum- ing that the vorticity generation in the channel stops at the breakwater corner, two separate arcs are defined to study the circulation growth of the inshore and offshore coherent structures. The (r;q;z) components of the cylindrical coordinates system correspond to the along-radius (radial), along the arc (tangential), and out-of-the-page (vertical) directions respectively. The derivatives ¶=¶q are evaluated along the arc r= 0, whereas¶=¶r derivatives are evaluated by defining arc offsets. This coordinate system set-up allows to impose no-slip (u r = 0) and solid wall (u q ;¶u q =¶q = 0) boundary conditions at the channel sides (q = 0;l). ¶u r =¶q is evaluated taking into account that R q 0 ¶u r =¶q 0 dq 0 ; R q l ¶u r =¶q 0 dq 0 = u r (q) (Didden, 1979). This coordinate system also serves for estimating the flow rate Q through the channel per unit depth h (Q= R l 0 ru r dq, m 3 =s), as shown in Figure 4.3, where l is the length of the arc. The vertical component of vorticity is defined in cylindrical coordinates as w z = 1 r ¶(ru q ) ¶r ¶u r ¶q = ¶u q ¶r + 1 r u q 1 r ¶u r ¶q ; (4.11) and therefore the vorticity influx (Eqn. 4.8) is written in polar coordinates as dG dt = Z l 0 ru r w z dq = Z l 0 ru r ¶u q ¶r + 1 r u q 1 r ¶u r ¶q dq = Z l 0 ru r ¶u q ¶r + u r u q dq: (4.12) Following the convention that the component normal to the domain boundary (u r ) has to be positive when the flow is entering the domain, the sign of u r is inverted to compute the vorticity influx into the offshore domain. 84 Integrating the vorticity influx (Eqn. 4.12) over time leads to G(t)= Z t 0 Z l 0 ru r ¶u q ¶r + u r u q dqdt 0 = Z t 0 Z l 0 ru r ¶u q ¶r dqdt 0 + Z t 0 Z l 0 u r u q dqdt 0 = I 1 +I 2 : (4.13) The individual contributions of terms I 1 and I 2 will be examined in the circulation growth sections of this chapter. The velocities along the arcs were evaluated using the spatial averaging technique (described in Section 3.10.2), with 20cm along-arc spacing and W R = 20cm. When the number of samples in the interrogation window were less than four (N sample < 4), velocities were interpolated using the natural-neighbor interpolation scheme. x (m) 5 10 15 20 y (m) 4 8 12 θ r (X c ,Y c ) θ=l inshore arc x (m) 5 10 15 20 y (m) 4 8 12 θ r (X c ,Y c ) θ=l offshore arc Figure 4.4: Polar coordinate system (r;q) used to compute circulation influx through the experi- mental channel. Coordinate z is the vertical (out of the page). The blue polygons define the areas over which the total circulation is evaluated in each side of the basin. 4.3 Vortex identification and quantification 4.3.1 Vortex identification Individual vortices in the flow are identified through the local velocity gradient tensor in a 2D flow given by D 2D = 2 4 ¶u ¶x ¶u ¶y ¶v ¶x ¶v ¶y 3 5 : (4.14) The tensor either has two real eigenvalues or a pair complex conjugate eigenvalues (Adrian et al., 2000). In areas where the local tensor undergoes swirling motion, the swirl strength l ci is equal to the imaginary part of the complex conjugate pair, otherwise l ci = 0. This measure, which 85 represents the local frequency of rotation (the period of one revolution corresponds to 2p=l ci ), has the advantage over vorticity maps in that shear layers are excluded (Adrian et al., 2000). Just as vorticity, swirl strength is a Galilean invariant quantity, and therefore the velocity field frame of reference doesn’t need to be stationary to evaluate the metric. 4.3.2 Circulation quantification The circulationG is defined as the integral of vorticity inside a closed surface S. In the case of a discretized flow field, it is approximated as the discrete sum of the cells inside the curve at each time step G vortex = Z Z S w z dS= dxdy N å i=1 w zi ; 8(x i ;y i ) i=1:N 2 A: (4.15) To measure the TCS circulation, A= A vortex is defined from the l ci = 0 contour enclosing the TCS-center. A vortex also provides an equivalent radius to quantify the TCS size as R vortex = r A vortex p : (4.16) Swirl-strength iso-contours were also used by Nicolau del Roure et al. (2009) and Bryant et al. (2012) to define starting-jet vortex boundaries. An alternative definition of the vortex boundary is given through the vorticity maps, by defining the limiting vorticityw c which defines the boundary as a fraction of the maximum vorticity at the TCS center w max (t): w c (t)=bw max (t), where b is a threshold factor close to zero. Given the resolution of the vorticity maps for this experiment, b = 0:02 gives the most reasonable results. V orticity-defined vortex boundaries have been used for the vortex ring experiments of Gharib et al. (1998); Krueger (2005); Aydemir et al. (2012) amongst others. In terms of the difference between the two vortex boundary definitions, the swirl-strength con- tour l ci = 0 encloses only the core area of the vortex, or in other words, it extends up to the radius for which the gradient¶u q =¶r is positive (at maximum u q ). This is illustrated in Figure 4.5 through an example for a shielded-Gaussian vortex with an azimuthal velocity profile of the type u q = r exp((r=R) 2 ) (see Equation 5.14b in Section 5.2). Swirl-strength is non-zero only in the 86 circular area around the center, where¶u q =¶r> 0. On the other hand, the vorticity contourw z = 0 (where w z (r) for cylindrical coordinates is given in Eqn. 4.11) extends further out in the radial direction due to the added u q =r term (for a non-axisymmetric vortex, the (1=r)¶u r =¶q vorticity term will also play a role). Therefore, the two vortex boundaries l ci = 0 and w c have different definitions and are not equivalent, even for simplified flow conditions (see Kolar (2007) for more details). If the kinetic energy beyond R vortex (radial distance to u q;max ) increases, without change in R vortex , thew c contour will expand but thel ci = 0 will not. 0 1 2 3 4 5 −0.2 0 0.2 0.4 0.6 0.8 (a) r (m) u θ (m/s) −5 0 5 −5 0 5 (b) λ ci −5 0 5 −5 0 5 (c) ω Figure 4.5: Radius vortex definition through an example for an axisymmetric vortex centered at (0;0). (a) Azimuthal velocity profile. (b) Swirl-strength map and velocity vectors. (c) V orticity map. Yellow circles are plotted every 1m in the radial direction. The magenda circle in (b) shows thel ci = 0 contour, and the red circle in (c) shows thew = 0 contour. The radii corresponding to the two contours are plotted in (a) with vertical dashed lines. It should also be noted here that any of the above vortex boundary definitions, in the early stages of the TCS development do not only include the vortex, but rather the starting-jet vortex as a whole. As an example of this statement, Figure 4.6 shows the swirl strength and vorticity maps for an ensemble vector field during TCS generation on the inshore side of the breakwater. It can be observed that the TCS boundary incorporates the trailing jet near the breakwater tip until it detaches from the trailing jet at what is referred to as the separation time (Aydemir et al., 2012) - this instant is different from pinch-off, which corresponds to the time at which the the jet ceases to supply circulation to the vortex (Dabiri, 2009). Therefore R vortex initially defines the starting-jet vortex radius and the vortex radius after pinch-off. 87 6 8 10 12 y (m) 16 18 20 22 24 x (m) t=39.2s (a) 16 18 20 22 24 x (m) ω −1 0 1 (1/s) (b) 16 18 20 22 24 x (m) λ ci −1 0 1 (1/s) (c) Figure 4.6: Swirl strength (b) and vorticity (c) maps corresponding to an ensemble flow field (a) during TCS generation at the inshore side of the breakwater. Swirl strength and vorticity contours are plotted with a step of 0.1 and 0.3 s 1 respectively. The lowest contours correspond tol ci = 0 andjwj= 0:15 s 1 . The total circulationG total in the inshore and offshore basin sides are evaluated using the poly- gon domains shown in Figure 4.4. The inshore circulation-evaluation area is limited by the spatial extent of the PTV data, and therefore does not cover a large part of the inshore basin. For the time window of the experiment considered in this section, the area encloses the TCS boundary at all times. 4.3.3 Upwelling/downwelling quantification This experimental setup doesn’t allow for the extraction of vertical velocities in the water column, since velocity data are only collected at the water surface (with the exception of the ADV loca- tions). The vertical gradient of the vertical velocity component can be estimated through the flow divergence metric for incompressible flow (conservation of mass equation),Ñ u= 0, and thus ¶w ¶z = ¶u ¶x + ¶v ¶y (4.17) Flow divergence maps highlight areas of flow convergence (positive¶w=¶z) and divergence (neg- ative ¶w=¶z), which correspond to flow downwelling and upwelling respectively. The vertical 88 flow rate Q vert can be estimated from the ensemble velocity field of spacing (Dx;Dy) using the approximate equation (Nicolau del Roure et al., 2009) Q vert = HDxDy ¶w ¶z =HDxDy ¶u ¶x + ¶v ¶y ; (4.18) where H is the water depth. In this formulation, upwelling flow rate is negative. 4.4 Inshore TCS generation The time-lapse of the inshore TCS generation is shown in Figures 4.7-4.9 using different flow metrics. The water column in the channel is initially unperturbed and flow separation starts on the offshore side of the breakwater right after the wave front reaches the channel. The growing separation zone, supplied with kinematic momentum from the jet in the channel, creates negative vorticity at the breakwater tip that with time extends to the lee side. V orticity is transported in the flow by advection and viscous diffusion. The vortex starts forming at the separation zone, then its center shifts just off the lee breakwater corner. It remains almost stationary until t 25s, while it is being supplied with circulation. Once the radius growth becomes limited by the breakwater vertical boundary the vortical structure starts propagating inshore. The vortex, coupled with the trailing jet constitute what is termed a starting vortex. When the flow starts reversing at t 38s, a starting-jet forms between the breakwater and the vortex that separates the vortex from the trailing jet. 4.4.1 Circulation growth The flow profiles during the generation of the inshore TCS across the width of the channel are shown in Figure 4.10. The initial profile (t = 10:1s) attains a shape that resembles closely that of an idealized (potential) flow through an orifice, boundary-conditions aside. Flow separation at the inshore tip of the breakwater creates the growing TCS, affecting the flow profile which starts deviating from the corresponding idealized flow profile, particularly close to the breakwater. The 89 6 8 10 12 y (m) t=18.1s t=20.5s 6 8 10 12 y (m) t=22.8s t=25.2s 6 8 10 12 y (m) t=27.5s t=29.8s 6 8 10 12 y (m) t=32.2s t=34.5s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=36.8s 8 10 12 14 16 18 20 22 24 x (m) t=39.2s Figure 4.7: Velocity vector maps at various times during the inshore TCS generation. 90 6 8 10 12 y (m) t=18.1s −1 0 1 (1/s) t=20.5s 6 8 10 12 y (m) t=22.8s t=25.2s 6 8 10 12 y (m) t=27.5s t=29.8s 6 8 10 12 y (m) t=32.2s t=34.5s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=36.8s 8 10 12 14 16 18 20 22 24 x (m) t=39.2s Figure 4.8: V orticity (w) maps at various times during the inshore TCS generation. Green polygons correspond to thel ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. 91 6 8 10 12 y (m) t=18.1s −1 0 1 (1/s) t=20.5s 6 8 10 12 y (m) t=22.8s t=25.2s 6 8 10 12 y (m) t=27.5s t=29.8s 6 8 10 12 y (m) t=32.2s t=34.5s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=36.8s 8 10 12 14 16 18 20 22 24 x (m) t=39.2s Figure 4.9: Swirl strength (l ci sign(w)) maps at various times during the inshore TCS generation. Green polygons correspond to the l ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. 92 deviation becomes more pronounced in later times (t > 20s) and eventually the flow in the channel starts reversing at t 38s. 10.1s 17.5s 24.8s 32.2s 39.5s (a) 0 0.2 0.4 0.6 0.8 1 θ/l -0.6 -0.4 -0.2 0 0.2 0.4 0.6 u r (m/s) θ 1 /l (b) 0 0.2 0.4 0.6 0.8 1 θ/l -0.6 -0.4 -0.2 0 0.2 0.4 0.6 u θ (m/s) (c) Figure 4.10: (a) Flow profiles across the channel for different times during the generation of the inshore TCS. (b,c) The corresponding velocity components in the polar coordinate system. Profiles are color-matched for time, and the aspect ratio of the flow profiles in (a) is skewed for presentation. The width of the separation zone at the breakwater tip becomes wider as the TCS size grows. This is evident from the distance of the maximum u r (at q =q 1 ) to the breakwater (q = 0). The distanceq 1 (drawn in Figure 4.10(b) for t= 24:8s) also defines the edge of the shear layer through which the TCS is supplied with circulation. Most of the area under the ¶u th =¶r curves is con- centrated in the regionq <q 1 , in similar fashion to the vortex ring generation through a cylinder (Didden, 1979). For q >q 1 , u r is more-or-less constant, attaining small ¶=¶q gradients. This is confirmed through the plot of the spatial gradient components (¶u r =¶q, ¶u q =¶u r ) of w z in Figure 4.11. The(1=r)¶u r =¶q component of vorticity is comparable to the¶u th =¶r component, but it does not contribute to the vorticity influx through the channel due to the boundary conditions. As a consequence, the top-wall boundary-layer, which introduces opposite-signed (positive) vorticity to the flow has also little effect on the inshore TCS generation. Circulation growth during the generation of the inshore TCS is shown in Figure 4.12. The total circulation in the inshore domain slowly builds-up until t 20s, it then reaches a steady slope of 1=5 and at t 40s it starts decreasing due to flow reversal. Total circulation computed from the vorticity influx terms of Eqn. 4.13 are shown in Figure 4.12a. V orticity influx evaluated at 93 0 0.2 0.4 0.6 0.8 1 θ/l -3 -2 -1 0 1 2 3 (1/r) ∂ u r /∂ θ (1/s) (a) 0 0.2 0.4 0.6 0.8 1 θ/l -3 -2 -1 0 1 2 3 ∂ u θ /∂ r (1/s) (b) 0 0.2 0.4 0.6 0.8 1 θ/l -3 -2 -1 0 1 2 3 ω z (1/s) (c) Figure 4.11: Velocity gradient and vorticity profiles in the channel for different times during the generation of the inshore TCS. Profiles are color-matched according to the times shown in Fig- ure 4.10. the inshore arc predicts circulation growth accurately until t 30s, at which point it starts being deficient. This growth deviation is attributed to the lack of accuracy/resolution of the flow near the breakwater tip. As it was demonstrated in the data processing chapter, the standard error of the mean flow is the largest in the flow separation region around the breakwater tip (see Figure 3.29). The standard error increases as the turbulence intensity in that region rises, and that affects the vorticity influx estimation. In terms of the weights of the two terms in the vorticity influx equation (Eqn. 4.15), term I 1 accounts for all the circulation growth until t 28s, when term I 2 starts contributing opposite- signed vorticity to the influx rate. Therefore, the importance of term I 2 is limited in this flow phase and can be ignored for the largest part. The circulation growth of the inshore TCS is shown in Figure 4.12b. All the basin circulation is supplied to the TCS up to time t 27:5s, after which point the inshore basin’s circulation grows faster than the TCS-enclosed circulation, even though the trailing jet is still attached to the TCS. Since the vorticity influx from the jet is shared with the trailing jet, TCS circulation growth becomes deficient compared to the total circulation (Afanasyev, 2006). The vortex eventually separates from the jet when the flow starts reversing at t 37:5s - the starting-jet being developed between the breakwater tip and inshore TCS acts as a barrier. This time corresponds to t 2, consistent with the observation of Nicolau del Roure et al. (2009) for tidal starting-jet dipoles 94 10 20 30 40 50 time (s) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -Γ (m 2 /s) (a) Γ total I 1 I 1 +I 2 10 20 30 40 50 time (s) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -Γ (m 2 /s) (b) Γ total Γ TCS 0 1 2 3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t*=(t-t 0 )U/W Figure 4.12: Circulation growth during the inshore TCS generation. (a) Total circulation compared against the vorticity influx parameters of Eqn. 4.13. (b) Total circulation compared with the TCS- enclosed circulation, as defined in Eqn. 4.15. Arrow denotes the pinch-off time of the TCS. forming in symmetric inlets. The TCS subsequently stops being supplied with vorticity from the jet, and the enclosed circulation gets damped by bottom friction as its being advected towards the channel. 4.5 Main (offshore) TCS generation The time-lapse of the offshore TCS generation is shown in Figures 4.16-4.19 using different flow metrics and with flow-visualization pictures in Figures 4.20- 4.21. The flow inshore of the break- water tip at t = 43:5s is reminiscent of a jet dipole created in a symmetric channel opening (e.g. Fig.3 in Afanasyev (2006)). The dipole and jet become skewed as the jet goes through the channel. Circulation magnitude is growing on the offshore side of the dipole as vorticity is being gener- ated at the inshore breakwater tip due to flow separation. Simultaneously the vortex grows in size, while being almost stationary, until its equivalent radius is no longer limited by the distance to the breakwater (at t 50:5s). The vortex then starts propagating offshore leaving a trailing jet of vorticity behind it. At t 55:2s, the vortex shows signs of separating from the trailing jet and at t 57:5s, the first separation occurs. The trailing jet takes the form of a turbulent wake, with a 95 vortex street pattern of single-signed vorticity. From the vortex street, secondary vortices form, which eventually merge with the main vortex in a series of straining events that deform the vortex shape (Mcwilliams, 1990). As will be discussed in the next section, the secondary vortices supply the main TCS with additional circulation, after initial separation from the trailing jet. As the jet is directed through the channel at t 48:2s, the clockwise-rotating inshore TCS is being advected with the flow. It’s initial circular shape (in plan view) becomes elliptical as it gets confined between the counter-clockwise-rotating TCS and the top wall. The interaction and merging of the opposite-signed (negative) vorticity with the main TCS (positive vorticity) can be better visualized through the dye-images in Figures 4.20- 4.21, in which negative vorticity advected near the top wall is colored red. At t = 55:2s the two flows start mixing through a meandering pattern. The jet - advecting negative vorticity - propagates faster than the vortex structure and surrounds the vortex through the front. The turbulent mixing of the two flows distorts the geometry of the TCS, which takes an elliptical shape at t 57:5s. The flow merging also creates distinct areas of upwelling/downwelling at the front of the TCS, where the mixing is more pronounced, as seen in Figure 4.19 for t 59:964:5s. This pattern of “frontal circulation” has also been observed at the front of propagating vortex dipoles (Sous et al., 2005; Akkermans et al., 2008). The frontal small-scale 3D turbulence is being transported to the main Q2D flow (Sous et al., 2004). 4.5.1 Circulation growth The flow acceleration through the channel during the offshore TCS generation is much more rapid than its inshore counterpart’s. The intensity of the return flow is also more pronounced, due to the faster wavemaker retraction and the reflection of the leading elevation wave on the basin’s back wall. The flow velocity and velocity gradient profiles during the generation of the offshore TCS are shown in Figures 4.13& 4.14. Continuing from the last channel-velocity profile shown in Figure 4.10a, the flow rapidly reverses (within 5s) to create radial velocities of the order of 0:8m=s on the offshore breakwater tip. The profile at t= 53:4s becomes more uniform across the channel, with a pronounced separation-zone near the breakwater tip. The separation-zone narrows 96 as the TCS moves away from the breakwater tip, and eventually the channel flow direction starts reversing again at 78s. 45.0s 53.4s 61.7s 70.0s 78.4s (a) 0 0.2 0.4 0.6 0.8 1 θ/l -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 u r (m/s) (b) 0 0.2 0.4 0.6 0.8 1 θ/l -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 u θ (m/s) (c) Figure 4.13: Flow profiles across the channel for different times (a) during the generation of the offshore TCS. The corresponding velocity components in the curvilinear coordinate system are plotted in (b) and (c) with continuous lines. Profiles are color-matched for time, and the aspect ratio of the flow profiles in (a) is skewed for presentation. The inshore TCS, which is being advected through the channel, can be observed passing through the offshore arc at t = 53:4s, as it distorts the uniformity of the u r velocity profile. It can also be observed in the vorticity profiles shown in Figure 4.14c, where a patch of negative vorticity is being advected at t= 53:4s. In all other flow profiles, vorticity is concentrated near the breakwater tip, with peak vorticity reaching values greater than 3s 1 at t= 53:4s. 0 0.2 0.4 0.6 0.8 1 θ/l -4 -3 -2 -1 0 1 2 3 4 (1/r) ∂ u r /∂ θ (1/s) (a) 0 0.2 0.4 0.6 0.8 1 θ/l -4 -3 -2 -1 0 1 2 3 4 ∂ u θ /∂ r (1/s) (b) 0 0.2 0.4 0.6 0.8 1 θ/l -4 -3 -2 -1 0 1 2 3 4 ω z (1/s) (c) Figure 4.14: Velocity gradient and vorticity profiles in the channel for different times during the generation of the offshore TCS. Profiles are color-matched according to the times shown in Fig- ure 4.10. 97 Circulation growth during the generation of the offshore TCS is shown in Figure 4.15. The total circulation in the offshore domain rises with an almost constant rate of dG=dt 0:55m 2 , until t 60s, with the exception of a mild dip between t = 50 56s due to the advection of the inshore TCS through the channel. It then continues rising at a slower pace until peak circulation is reached at t 63s, about 15s before flow reversal in the channel. The time of peak circulation represents the point at which the rate of vorticity influx through the channel becomes less than the rate of energy loss due to bottom friction in the offshore domain. Circulation growth estimated from the vorticity influx terms of Eqn. 4.13 closely follows the offshore total circulation curve until t 50s (Figure 4.15a). The circulation growth deficiency starting at t 50s coincides with the time that the mean-channel flow reaches its peak value and when the negative vorticity of the inshore TCS reaches the offshore arc, advecting swirling flow. The high turbulent intensity of the flow conditions at this time can again be associated with the uncertainty (and resolution) of the mean flow evaluation, and thus with the estimation of vorticity influx through the channel. In terms of the relative weight of the two terms in the vorticity influx equation (Eqn. 4.15), term I 1 accounts for all the circulation growth until t 46s, when term I 2 starts contributing positive-signed vorticity to the influx rate, thus further increasing the circulation growth. The fact that term I 2 adds more circulation during the offshore TCS generation, as opposed to the inshore TCS generation where term I 2 reduced the circulation growth, this is due to the skewed angle of the breakwater, creating favorable u q velocity profiles during the return flow. For a right-angled breakwater, the contribution of term I 2 is expected to be negligible, as it is customarily assumed when applying the slug model (Eqn. 4.4). The circulation growth of the offshore TCS is shown in Figure 4.15b. All the basin circulation is supplied to the TCS up to time t 55s. In fact,G TCS is initially larger thanG total since the area of the trailing jet, which is included in thel ci contour, extends beyond the offshore arc (inshore of the arc). The vortex pinches-off from the jet at t 56s as seen from both Figures 4.18& 4.15b. The TCS is periodically re-supplied with circulation through merging with secondary vortices shed from the vortex street and retains (but doesn’t exceed) the circulation at separation. After each 98 50 60 70 80 time (s) 0 2 4 6 8 10 12 Γ (m 2 /s) dΓ/dt∼ 0.55 (a) Γ total I 1 I 1 +I 2 50 60 70 80 time (s) 0 2 4 6 8 10 12 Γ (m 2 /s) (b) Γ total Γ TCS 0 1 2 3 4 0 2 4 6 8 10 12 t*=(t-t 0 )U/W Figure 4.15: Circulation growth during the offshore TCS generation. (a) Total circulation com- pared against the vorticity influx parameters of Eqn. 4.13. (b) Total circulation compared with the TCS-enclosed circulation, as defined in Eqn. 4.15. Arrows denote the pinch-off time of the TCS and the time of flow reversal. merging episode, the TCS circulation decays due to bottom friction, and flow mixing with the opposite-signed vorticity jet near the top wall. Similar circulation decay in-between mergers was noted by Mcwilliams (1990) for vortices in two-dimensional turbulence, and it was attributed to either casting off of vorticity filaments in the axisymmetrization process or due to the stripping of weak vorticity patches from the periphery of a vortex due to the strain from neighbouring vortices. The time of pinch-off corresponds to t 1:5 which is less than the expected value of t 2. This deviation can be attributed to either the asymmetry of the channel-velocity periodic signal, and/or the asymmetry of the basin geometry. Note that the velocity scaling is larger in the return flow, and thus dimensionless time is disproportional to the equivalent of the inshore TCS generation phase. The dye-visualization images in Figure 4.20 clearly show that the TCS remains attached to the trailing jet at least until t 69:2s. Separation from the jet is estimated to occur at t 76s by inspection of Figure 4.21. 99 6 8 10 12 y (m) t=43.5s t=45.8s 6 8 10 12 y (m) t=48.2s t=50.5s 6 8 10 12 y (m) t=52.9s t=55.2s 6 8 10 12 y (m) t=57.5s t=59.9s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=62.2s 8 10 12 14 16 18 20 22 24 x (m) t=64.5s Figure 4.16: Velocity vector maps at various times during the offshore TCS generation. 100 6 8 10 12 y (m) t=43.5s −1 0 1 (1/s) t=45.8s 6 8 10 12 y (m) t=48.2s t=50.5s 6 8 10 12 y (m) t=52.9s t=55.2s 6 8 10 12 y (m) t=57.5s t=59.9s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=62.2s 8 10 12 14 16 18 20 22 24 x (m) t=64.5s Figure 4.17: V orticity (w) maps at various times during the offshore TCS generation. Green poly- gons correspond to thel ci = 0 contour around the TCS center (red circles) and black circles cor- respond to the smallest-fitting circle from the TCS center to the breakwater. 101 6 8 10 12 y (m) t=43.5s −1 0 1 (1/s) t=45.8s 6 8 10 12 y (m) t=48.2s t=50.5s 6 8 10 12 y (m) t=52.9s t=55.2s 6 8 10 12 y (m) t=57.5s t=59.9s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=62.2s 8 10 12 14 16 18 20 22 24 x (m) t=64.5s Figure 4.18: Swirl strength (l ci sign(w)) maps at various times during the offshore TCS gener- ation. Green polygons correspond to the l ci = 0 contour around the TCS center (red circles) and black circles correspond to the smallest-fitting circle from the TCS center to the breakwater. 102 6 8 10 12 y (m) t=43.5s −25−20−15−10−5 0 5 10 15 20 25 t=45.8s 6 8 10 12 y (m) t=48.2s t=50.5s 6 8 10 12 y (m) t=52.9s t=55.2s 6 8 10 12 y (m) t=57.5s t=59.9s 6 8 10 12 y (m) 8 10 12 14 16 18 20 22 24 x (m) t=62.2s 8 10 12 14 16 18 20 22 24 x (m) t=64.5s Figure 4.19: Upwelling/downwelling maps at various times during the offshore TCS generation. Units shown in the colorbar correspond to 10 3 cm 3 =s. Upwelling has negative value. Red circles indicate the TCS center. 103 t=43.5s t=45.8s t=48.2s t=50.5s t=52.9s t=55.2s t=57.5s t=59.9s t=62.2s t=64.5s t=66.9s t=69.2s Figure 4.20: Images from a dye visualization experiment captured during the offshore TCS gener- ation. Fluorescent green dye is released from the breakwater tip (inside the separation zone) and fluorescent red dye is released just inshore of the tip. The images are taken from an oblique angle and have not been undistorted/rectified. 4.5.2 Other characteristics of the starting-jet vortex Figure 4.22 shows additional metrics to describe the jet dynamics and Figure 4.23 shows the tra- jectory of the starting-jet vortex centers. The offshore vortex x and ydisplacement time-history shows that the vortex initially has the tendency to advect towards the top-wall, as it moves away from the breakwater to make room for its spatial growth. and subsequently diverges away from it. This change of course is related to the vortex size growth and the presence of the jet shear layer 104 t=52.9s t=55.2s t=57.5s t=59.9s t=62.2s t=64.5s t=66.9s t=69.2s t=71.5s t=73.9s t=76.2s t=78.5s Figure 4.21: Figure 4.20 continued. Images taken from another camera and angle. pushing it away from the top wall. As the kinetic energy of the jet decays at t 69s, the vortex starts moving towards the top-wall again, also being expelled by the jet of the weaker return flow. Figure 4.22a shows that the vortex radius grows linearly against the formation time t until pinch-off from the trailing jet. Following Afanasyev (2006), the growth can be related to the jet fluid supply rate U jet W=2 (per unit depth). Note that W is defined here as double the channel width, to be compatible with the definition of formation time t (Eqn. 4.10) given for symmetric inlets. Since the fraction of fluid supply being integrated in the vortex is yet unknown, the supply rate is multiplied by a constant factorg, to be determined. 105 The vortex area A vortex can be assumed - to first order - to increase linearly with time (Afanasyev, 2006): A vortex = A 0 +U jet W 2 gt: (4.19) The jet speed U jet for non-impulsive forcing can be defined as the running mean of an equivalent piston velocity (Gharib et al., 1998) U jet = ¯ U p = 1=t R t 0 ¯ U(t 0 )dt 0 (where ¯ U is the mean channel velocity), and it is related to the formation timet by ¯ U p =tW=t. Substituting for U jet , Eqn. 4.19 becomes A vortex = A 0 +t W 2 2 g; (4.20) and thus the starting-jet vortex radius grows like R vortex p A vortex =p W q g 2p p t. The experi- mental data follow a growth rate of dR vortex =d p t 5=4, which results in g = 3:13p=W 2 1=4. Therefore, to first order, a quarter of the flow rate through the channel is supplied to the vortex core. From the vortex center xdirection displacement evolution, three phases can be identified: the time interval t= 4049s, in which the starting-jet vortex moves away from the breakwa- ter (see vertical displacement curve) while growing in size. The velocity in the xdirection is0:15m=s. the time interval t = 49 60s, in which the starting-jet vortex is mostly being advected by the jet. The velocity in the xdirection is0:33m=s. Plotting the vortex displacement against the equivalent piston stroke L, reveals a linear relationship, with a slope of 0:56. This observation is in agreement with Afanasyev (2006), who found a relationship with slope 0:5 0:02, using both potential flow theory and experimental data for vortex dipoles. the time interval t = 60 73s, in which the vortex slows-down until t = 67s and then accel- erates again reaching a horizontal translation speed of0:25m=s. 106 0 0.5 1 0 0.5 1 1.5 2 2.5 R/ √ τ ∼ 5/4 √ τ = ! L/W R vortex (m) (a) 0 5 10 15 −22 −20 −18 −16 −14 −12 dx/dL∼ 0.56 equivalent piston stroke, L (m) −x (m) (b) 50 60 70 12 14 16 18 20 22 dx/dt∼−1/3 dx/dt∼−0.15 time (s) x (m) (c) 50 60 70 9.6 9.8 10 10.2 10.4 10.6 10.8 time (s) y (m) (d) Figure 4.22: Metrics describing the offshore starting-jet vortex dynamics. (a) V ortex-equivalent- radius growth as a function of formation timet. (b) V ortex horizontal displacement Vs equivalent piston stroke. (c) Horizontal displacement. (d) Vertical displacement. 12 14 16 18 20 22 24 8 9 10 11 12 13 x (m) y (m) offshore TCS inshore TCS Figure 4.23: Trajectories of the centers of the starting-jet vortices for the isnhore and offshore flow phases. 107 Chapter 5 TCS characteristics, evolution and flow tran- sition to Q2D The previous chapter examined the flow properties that lead to the generation of the offshore coher- ent structure, i.e. during the experimental time period t = 0 80s. This chapter is concerned with the evolution of the coherent structure after formation. The velocity field is represented in a TCS- centered coordinate system: the TCS center from the individual trials is matched and the velocity vectors are shifted accordingly. This representation allows to account for the variation in the TCS- center paths (Figure 3.24) due to the effects of the chaotic/turbulent nature of the flow. The path deviation becomes significant after t 65s and the experimental time period examined here corre- sponds to t= 65 3000s, allowing some overlap with the analysis of the previous chapter. The 2D vortex structure is determined from the scattered velocity data and is characterized using appropriate geophysical vortex profiles from the literature. The secondary flow is also char- acterized and quantified from the available observations. The spatial growth of the vortex due to viscous diffusion and kinematic energy decay due to bottom friction are examined later, and at the end of the chapter an attempt is made to determine the flow transition point from three-dimensional (i.e. with significant secondary flow) to quasi-two-dimensional (Q2D). 108 5.1 Governing equations Turbulent shallow water flows with large horizontal to vertical scale ratios L=H >> 1, implying to the hydrostatic approximation, are typically modeled using the depth-averaged shallow-water equations. For a 2D turbulent flow with depth h(x;y) over a rigid bottom layer of elevation z b (x;y) and no background rotation, the governing equations are described by ¶h ¶t + ¶(h ¯ u) ¶x + ¶(h ¯ v) ¶y ; (5.1) ¶ ¯ u ¶t + ¯ u ¶ ¯ u ¶x + ¯ v ¶ ¯ u ¶y =g ¶ ¶x (h+ z b )+ n e f f h ¶ ¶x h ¶ ¯ u ¶x + n e f f h ¶ ¶y h ¶ ¯ u ¶y t bx rh ; (5.2) ¶ ¯ v ¶t + ¯ u ¶ ¯ v ¶x + ¯ v ¶ ¯ v ¶y =g ¶ ¶y (h+ z b )+ n e f f h ¶ ¶x h ¶ ¯ v ¶x + n e f f h ¶ ¶y h ¶ ¯ v ¶y t by rh ; (5.3) where ¯ u; ¯ v are the depth-averaged horizontal velocities (Socolofsky and Jirka, 2004; Seol and Jirka, 2010). n e f f is the effective viscosity, given as the sum of turbulent and molecular kinematic con- tributions: n e f f = n turb +n - adding more dissipation/diffusion to the flow description due to turbulence. The molecular kinematic viscosity for water is taken asn 1:12 10 6 m 2 =s and the added viscosity due to turbulence can be modeled using n turb u H r c f 2 ¯ UH; (5.4) where u = ¯ U p c f =2 is the bottom shear velocity, ¯ U is a reference velocity and c f is the bed friction coefficient (Seol and Jirka, 2010). The n e f f terms of the momentum equations represent lateral diffusion -by entrainment of ambient fluid at the vortex boundaries. Vertical diffusion due to the no-slip boundaries is repre- sented by the bottom shear stress terms t bx ;t by , which are computed using the quadratic friction law t bx =r c f 2 ¯ u p ¯ u 2 + ¯ v 2 ; t by =r c f 2 ¯ v p ¯ u 2 + ¯ v 2 : (5.5) The bed friction coefficient is defined as c f =l=4, andl is the Darcy-Weisbach friction coeffi- cient to relate the pipe roughness to channel flows. c f incorporates the effects of bottom roughness 109 and internal visous friction (Carmer, 2005). For a logarithmic velocity profile and smooth bottom it can be obtained from the standard smooth-wall relation as (Uijttewaal and Booij, 2000) 1 p c f = 1 k (log(Re h p c f )+ 1:0); (5.6) where k = 0:4 is the von Karman constant, and Re H = ¯ UH=n is the depth Reynolds number, evaluated using the depth-averaged velocity ¯ U. Typical values for c f are c f 0:005;0:01 for the field and laboratory, respectively (Socolofsky and Jirka, 2004). In this experiment, ¯ U is in the range of 0:01 1 m/s, resulting in c f 0:003 0:018. For laminar flows, the vertical velocity profile becomes Poiseuille-like (Satijn et al., 2001) and the bottom friction term is modeled by (Seol and Jirka, 2010) t b = 3nr ¯ U H : (5.7) In the 2D dimensionless Navier-Stokes momentum equations ¶u ¶t +(uÑ)u=Ñp+ 1 Re Ñ 2 u 1 Re l u; (5.8) the laminar flow bottom friction term is represented by Re l =w=l, a Reynolds number associated with vertical diffusion (Satijn et al., 2001). w corresponds to a characteristic vorticity scale in the flow, whereas l =p 2 n=4H 2 is the external (or Rayleigh) friction factor. Re= L 2 w=n is the common Reynolds number for a characteristic length scale L. The transition time-scale from turbulent to laminar flow, can be estimated from the vertical vis- cous relaxation time 4H 2 =9p 2 n (J¨ uttner et al., 1997). For this experiment this time-scale translates to t 3:4hrs, which extends beyond the duration of the trials. Therefore, flow decay properties for laminar conditions are not examined here. However, the bulk of the theory for shallow-water monopolar vortices has been devised for laminar, or turbulent conditions with strong shallowness where the flow profile quickly relaxes to a Poiseuille-like profile. 110 5.2 Monopolar vortex theory This section summarizes the governing equations for a purely azimuthal vortex flow, with no back- ground rotation, for which analytical solutions exist. While the secondary (three-dimensional) flow is strong, it is expected that the assumptions of axisymmetry and purely azimuthal flow will be vio- lated. The deviation from the theory will provide a basis to quantify the effect of the secondary flow on the main (2-D) vortex structure. Assuming a purely azimuthal flow (u r ;u z = 0) with no background rotation, the radial and axial components of the cylindrical Navier-Stokes equations are reduced to u 2 q r = 1 r ¶ p ¶r ; (5.9a) 1 r ¶ p ¶z = g; (5.9b) which represent the cyclostrophic and hydrostatic balance respectively. The azimuthal component of the momentum equations is given by ¶u q ¶t =n 1 r ¶ ¶r r ¶u q ¶r u q r 2 + ¶ 2 u q ¶z 2 ; (5.10) and represents the radial diffusion equation. The corresponding horizontal vorticity diffusion equa- tion (neglecting vertical diffusionn¶ 2 w z =¶z 2 ) is given by ¶w z ¶t =nÑ 2 w z =n 1 r ¶ ¶r r ¶w z ¶r : (5.11) Solving the radial diffusion equation for u q using separation of variables, and assuming a Poiseuille velocity profile, yields a radial-temporal part of the form (Satijn et al., 2001) ˆ R(r;t)= Z ¥ 0 a(p)J 1 (pr)exp(n p 2 t)dt; (5.12) wherea(p) is determined from the initial conditions profile G(r) using Hankel integrals as (Beck- ers et al., 2001) a(p)= p Z ¥ 0 G(u)J 1 (pu)udu: (5.13) 111 The choice for the initial condition’s profile depends on the generation mechanism. Typical profiles for geophysical vortices include the Lamb-Oseen and shielded Gaussian (van Heijst and Clercx, 2009) u q (r)= w 0 R 2 2r 1 exp r 2 R 2 ; (5.14a) u q (r)= 1 2 w 0 r exp r R 2 ; (5.14b) where R is a length scale and w 0 is the maximum vorticity at the vortex center (r = 0). The corresponding vorticity profiles (w =(1=r)¶=¶r(ru q )) are given by w(r)=w 0 exp r 2 R 2 ; (5.15a) w(r)=w 0 1 r 2 R 2 exp r R 2 : (5.15b) The Lamb-Oseen vortex is generated in the lab using the “sink technique”. This technique involves sucking fluid away through a vertically-aligned perforated tube and is used to produce cyclonic vortices (van Heijst and Clercx, 2009). The fluid deflection when moving towards the source due to conservation of angular momentum mimics the Coriolis acceleration. The shielded Gaussian (or isolated) vortex is typically generated in the lab using the “stirring technique”. The stirring technique involves confining fluid inside a rotating cylinder (with no-slip boundaries) which is lifted once a purely azimuthal flow is achieved inside the cylinder. The ini- tially still ambient fluid interacts with the rotating fluid to create a zone of vorticity with an opposite sign than the vortex core. It has been shown that any vortex with some level of axisymmetry and zero initial circulation will eventually evolve into this particular profile (Kloosterziel, 1990), and that it’s impossible to generate a monopolar vortex of single-signed vorticity (Satijn et al., 2001). A more general profile for stirring-type vortices, is given by van Heijst and Clercx (2009) as u q (r)= 1 2 w 0 r exp r R a ; and, (5.16a) w(r)=w 0 1 1 2 a r R a exp r R a : (5.16b) 112 where a is a parameter controlling the shape of the profile - the steepness increases with increasing a. For this profile, the radius of u q;max is given by R max = R a 1=a : (5.17) Stability analysis on this family of isolated vortices has shown that the profile becomes unstable for a> 1:85 (Seol and Jirka, 2010). A third vortex generation method, which is very similar to the stirring technique, consists of injecting fluid along the inner wall of an open thin-walled cylinder submerged in the fluid (Flor and Van Heijst, 1996). The process of injecting fluid stops and after allowing for the flow to evolve into regular rotational motion the cylinder is lifted. This technique is called the tangential injection method. The resulting vortex is another type of a shielded and isolated monopolar vortex with a fitting dimensionless azimuthal profile of the form (Flor and Van Heijst, 1996) ˆ u q (ˆ r)= ˆ r exp 1 ˆ r a a ; (5.18) where ˆ u q is non-dimensionalized using the maximum azimuthal velocity ( ˆ u q = u q =u q;max ) and the radial distance by the radial distance R max corresponding to u q;max , as ˆ r = r=R max . For this profile, the radius of u q;max is given by ˆ r= 1, or R max = R. Thus, contrary to the “stirring” profile, for the tangential “injection” profile there is a direct correspondence between R and R max . The corresponding vorticity profile is given by ˆ w(ˆ r)= 2 1 1 2 ˆ r a exp 1 ˆ r a a ; (5.19) where ˆ w =w=w max . a is steepness parameter similar to the “stirring” profile. Carton et al. (1989) found this type of vortex to be unstable for a> 2. All stirring-generated profiles have zero circula- tion due to the flow confinement inside the cylinder. This can be simply confirmed as G= 2p Z ¥ 0 rw z dr= 2p Z ¥ 0 r 1 r ¶ ¶r (ru q )dr= 2p Z ¥ 0 ¶ ¶r (ru q )dr= 2p[ru q ] ¥ 0 ; (5.20) since the profiles decay exponentially at r!¥ and therefore lim r!¥ ru q (r)= 0. 113 Finally, another similarity solution of the radial diffusion equation on the fplane (with con- stant background rotation) is given by (Fl´ or and Eames, 2002) ˆ u q (ˆ r)= 1 ˆ r exp 1 4 ˆ r 2 W1 2 (a1); 1 2 1 2 ˆ r 2 ; (5.21) where ˆ r and ˆ u q are defined as for the tangential “injection” profile and W is the Whittaker function (e.g. Abramowitz and Stegun (1972), pp. 505) given the free parameter a - the steepness of the velocity profile increases monotonically with a. This family of functions produces the shielded Gaussian and Lamb-Oseen vortex profiles for a= 3 and a= 1, respectively. “Stirring” vortices in general are characterized by 1:4< a 3, and vortices generated using the sink technique are described by a values close to 1 (Fl´ or and Eames, 2002). 5.3 2D flow field evolution This section examines the properties of the flow after the TCS formation. The flow description is provided through the interpretation of the vorticity and divergence maps shown in Figures 5.1-5.2 corresponding to the early development stage (in the time period t= 66129s) and Figures 5.3-5.4 showing later time steps. The vorticity map at t = 66s shows that the TCS is still in an early-development stage. Some of the low-level vorticity contours are still connected with the trailing jet, and the TCS is currently being merged with a secondary vortex from the vortex street. The vortex is characterized by a compact core of intense positive vorticity surrounded by a ring of opposite-signed vorticity. This ring forms from the interaction (shearing) of the vortex boundary with the surrounding irrotational fluid. The vortex boundary expands in time via viscous diffusion to entrain more fluid inside its bulk structure. The advected negative vorticity of the inshore TCS can be identified at the top left corner of the domain at t = 66s. This fluid patch rotates around the vortex structure as time progresses, further reinforcing the strength of the opposite-signed ring. The “wrapping” of the opposite-signed vorticity can also be observed in the dye images of Figures 4.20-4.21. The 114 ring provides the support for the development of an “isolated” vortex (i.e. a vortex with zero net circulation) and can still be observed at later time steps shown in Figure 5.3. −6 −4 −2 0 2 4 6 (m) t = 66 s −1 0 1 (1/s) t = 79 s −1 0 1 (1/s) t = 92 s −1 0 1 (1/s) −6 −4 −2 0 2 4 6 (m) −6 −4 −2 0 2 4 6 (m) t = 104 s −1 0 1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 117 s −1 0 1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 129 s −1 0 1 (1/s) Figure 5.1: TCS-centered vorticity maps. Figure 5.2 shows the corresponding flow divergence maps (with inverted sign). In this plotting format, the vertical velocity gradient¶w=¶z and flow downwelling are positive. From the maps, it can be observed that the TCS center is undergoing downwelling. This is expected, since the surface elevation at the center exhibits a strong depression (low pressure), which is evidence of flow con- vergence and downwelling. The TCS core is initially surrounded by a ring of upwelling, indicative of a 3D re-circulation ring about the horizontal axis. This vortical re-circulation affects the hori- zontal momentum balance, since it brings water of low-momentum from the bottom to the surface and is also linked with sediment transport (Nicolau del Roure et al., 2009). At later times shown in Figures 5.2& 5.4, the upwelling flow ring is interrupted by patches of downwelling extending out from the center in various shapes. These patterns indicate that the flow is not axisymmetric and that the radial velocity has a strong azimuthal gradient. 115 −6 −4 −2 0 2 4 6 (m) t = 66 s −0.2−0.1 0.0 0.1 0.2 (1/s) t = 79 s −0.2−0.1 0.0 0.1 0.2 (1/s) t = 92 s −0.1 0.0 0.1 (1/s) −6 −4 −2 0 2 4 6 (m) −6 −4 −2 0 2 4 6 (m) t = 104 s −0.1 0.0 0.1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 117 s −0.1 0.0 0.1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 129 s −0.1 0.0 0.1 (1/s) Figure 5.2: TCS-centered flow divergence (¶w=¶z=(¶u=¶x+¶v=¶y)) maps. Downwelling is positive. 5.4 Azimuthal-averaged properties In this section, the azimuthal-averaged properties of the experimental TCS are presented and ana- lyzed. Evaluating the vortex profiles in an azimuthal-averaged framework allows one to filter-out non-axisymmetric features in the flow 2D fields shown in the previous section. The profiles were obtained using the procedure outlined in Section 3.10.2. The profile radii are limited by the dis- tance to the closest vertical boundary. The TCS azimuthal-averaged velocity components(u r ;u q ) and the corresponding vorticity pro- files are plotted in Figures 5.6-5.11 at 50s time intervals, up to t = 900s. The vorticity profiles - to the extent of radii considered - exhibit a smooth decay with no identifiable small-scale features. The coarse flow resolution and azimuthal-averaging filtered any high wavenumber features out. The profiles also do not exhibit negative values but tend towards zero close to the domain bound- ary. Since the last radial point in the profile is the minimum distance to a boundary out of all the 116 −6 −4 −2 0 2 4 6 (m) t = 142 s −1.0−0.50.0 0.5 1.0 (1/s) t = 154 s −1.0−0.50.0 0.5 1.0 (1/s) t = 167 s −1.0−0.5 0.0 0.5 1.0 (1/s) −6 −4 −2 0 2 4 6 (m) t = 179 s −1.0−0.5 0.0 0.5 1.0 (1/s) t = 192 s −0.8−0.40.0 0.4 0.8 (1/s) t = 204 s −0.8−0.40.0 0.4 0.8 (1/s) −6 −4 −2 0 2 4 6 (m) t = 700 s −0.2−0.1 0.0 0.1 0.2 (1/s) t = 750 s −0.2−0.1 0.0 0.1 0.2 (1/s) t = 800 s −0.2−0.1 0.0 0.1 0.2 (1/s) −6 −4 −2 0 2 4 6 (m) −6 −4 −2 0 2 4 6 (m) t = 850 s −0.2−0.1 0.0 0.1 0.2 (1/s) −6 −4 −2 0 2 4 6 (m) t = 900 s −0.2−0.1 0.0 0.1 0.2 (1/s) −6 −4 −2 0 2 4 6 (m) t = 950 s −0.2−0.1 0.0 0.1 0.2 (1/s) Figure 5.3: TCS-centered vorticity maps. Continued from Figure 5.1. 117 −6 −4 −2 0 2 4 6 (m) t = 142 s −0.1 0.0 0.1 (1/s) t = 154 s −0.1 0.0 0.1 (1/s) t = 167 s −0.1 0.0 0.1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 179 s −0.1 0.0 0.1 (1/s) t = 192 s −0.1 0.0 0.1 (1/s) t = 204 s −0.1 0.0 0.1 (1/s) −6 −4 −2 0 2 4 6 (m) t = 700 s −0.02 −0.010.00 0.01 0.02 (1/s) t = 750 s −0.02 −0.010.00 0.01 0.02 (1/s) t = 800 s −0.02 −0.010.00 0.01 0.02 (1/s) −6 −4 −2 0 2 4 6 (m) −6 −4 −2 0 2 4 6 (m) t = 850 s −0.02 −0.010.00 0.01 0.02 (1/s) −6 −4 −2 0 2 4 6 (m) t = 900 s −0.02 −0.010.00 0.01 0.02 (1/s) −6 −4 −2 0 2 4 6 (m) t = 950 s −0.02 −0.010.00 0.01 0.02 (1/s) Figure 5.4: TCS-centered flow divergence maps. Continued from Figure 5.2. 118 trials considered, for many of the trials contributing velocity vectors, the actual boundary extends further out in radius. The profile does not extend to the opposite-signed vorticity ring described in the previous section. The u q profiles of Figures 5.6-5.7 show a steady decay of the maximum, from u q;max (t = 100) 0:56m=s to u q;max (t= 950) 0:11m=s. As will be shown in Section 5.5.1, the decay is the result of bottom friction. None of the u q profiles approaches zero in the furthest radius considered. This is expected since vorticity is positive across the radial spectrum and thus the circulation of the vortex maintains a positive value. The maxima of the u r velocity profiles are an order of magnitude smaller compared to the u q maxima, and also show a steady decay with time. The standard deviation of the samples in each evaluation radius is of the same order for both velocity components. However, the standard devia- tion as a percentage of the mean (coefficient of variation) is much greater for the radial component. Regardless of the large uncertainty in the u r azimuthal-averaging, the profiles serve well to visu- alize the flow divergence patterns. In the early stages, there is a strong convergence zone near the vortex center, extending up to r 1:7m. Convergence is followed by a zone of flow divergence in the outer vortex region. This flow pattern is in agreement with and explains the observations out- lined in Section 3.9.1. Similar patterns of surface flow convergence/divergence can be identified in aerial pictures of the coherent structures that formed in Oarai harbor during the 2011 Tohoku tsunami (Lynett et al., 2012). Quantifying the radial extent of flow convergence by evaluating the radius where u r crosses the zero-line leads to noisy results. Instead, it is defined here in terms of the ratio of the negative u r values N over the sample population N of each u r profile: R convergence = d min N N ; (5.22) where d min is the distance to the closest vertical boundary (or the maximum radius of the profile) defined in Eqn. 3.48. This definition yields reliable results up to t 380s but the later results become unstable due to the high coefficient of variation in the measurements. For the time period t= 100 380s, the convergence zone extends to R convergence = 1:65 0:27m. 119 Selected azimuthal-averaged u q and w z profile evolutions are plotted together in Figure 5.5 at a 200s time interval. Both profiles are normalized by their maxima in the ordinate and the u q abscissa is normalized by the radius R vmax corresponding to u q;max . Through this plotting format, it can be readily observed that the u q profiles maintain a constant distribution until t 500s. After t 500s the profile relaxes and the steepness at radii beyond R vmax decreases. On the other hand, the vorticity profiles show a steady lateral spreading from the start, which as shown in Section 5.5.2 is a result of viscous diffusion. 012345 r/R vmax 0 0.2 0.4 0.6 0.8 1 1.2 u θ /u θ,max t=100s t=300s t=500s t=700s t=900s t=1100s (a) 012345 r (m) 0 0.2 0.4 0.6 0.8 1 1.2 ω/ω max (b) Figure 5.5: TCS azimuthal-averaged u q and w profiles. (a) Azimuthal velocity (u q ) normalized with u q;max and its corresponding radius R vmax . (b) V orticity normalized withw max . The scattered azimuthal velocity data are fitted with the “stirring” a-profile given in Eqn. 5.16 and the “injection” a-profile given in Eqn. 5.18. The fit and resulting parameters were found to be equal. Since for the “injection” a-profile there is a direct correspondence R= R max it was chosen here for presentation. The resulting fit is shown in Figures 5.12-5.15. It can be observed that at early times the azimuthal velocity profiles are not accurately described using the “stirring” profile. The a-profile cannot model the steepness of the curve around u q;max and does account for the change in slope at radii beyond R vmax . The u q profiles start to converge towards the isolated vortex profile after t 500s. This is evident from the decrease of the root-mean-square-error of the fit, although the scattered data also become more compact, which influences the statistics of the fit. The steepness factor of the fitted profiles ranges between a= 0:33 0:43, and it decreases as the 120 experiment progresses, reflecting the relaxation of the vortex profile. All a steepness values are well below the critical value (a<< 1:85) and therefore the coherent structure is characterized as stable. The scattered data are also fitted with the Whittaker profile given in Eqn. 5.21, as shown in Figures 5.16-5.17. In contrast with the “injection” a-profile, this similarity solution matches the observations better at the early stages, and is unable to model the latter u q scatter distributions. The early Whittaker profiles match the steep peak at u q;max and also model the inflection point at larger radii. The a parameter of the Whittaker function starts increasing, reaching a peak value of a= 0:43 at t= 250s and then decreases down to a value of a= 0:24 at t= 750s. The initial increase of a does not necessarily represent an increase of the profile steepness, as there is no evidence of that in Figure 5.5. It rather represents the adjustment of the Whittaker profile to the progressive radial domain extension. The time when the a-profile matches the observations better than the Whittaker profile corre- sponds to t 450 500s, as evident from the root-mean-square-errors given in the corresponding figures. This time matches with the period when the radial velocity component starts becoming insignificant. As shown in Figure 3.20, the tracer-conglomerate at the TCS center starts expanding between t 320 485s. This period can be interpreted as a transition time for the flow regime of the vortex. The effect of the pronounced radial velocity component in the early life span of the vortex is to deform the azimuthal velocity profile, which is better described by a similarity solu- tion valid for the fplane. In the latter stages, the decay of the radial velocity profile allows the azimuthal velocity distribution to relax towards the a-profile. 5.5 Decay of mean vortex properties Taking the curl of the 2D shallow-water momentum Eqns. 5.2&5.3, using the continuity Eqn. 5.1 for flat bottom (z b = 0), and evaluating for a small amplitude h, leads to the conservation equation for the depth-integrated vorticity (Seol and Jirka, 2010) Dw Dt =n e f f ¶ 2 w ¶x 2 + ¶ 2 w ¶y 2 c f 2h w p ¯ u 2 + ¯ v 2 + ¯ u ¶ ¶y p ¯ u 2 + ¯ v 2 ¯ v ¶ ¶x p ¯ u 2 + ¯ v 2 : (5.23) 121 u θ (m/s) 0 0.2 0.4 0.6 t=100.1s u θ (m/s) 0 0.2 0.4 0.6 t=150.1s u θ (m/s) 0 0.2 0.4 0.6 t=200.1s u θ (m/s) 0 0.2 0.4 0.6 t=250.1s u θ (m/s) 0 0.2 0.4 0.6 t=300.1s u θ (m/s) 0 0.2 0.4 0.6 t=350.2s r (m) 012345 u θ (m/s) 0 0.2 0.4 0.6 t=400.2s r (m) 012345 u θ (m/s) 0 0.2 0.4 0.6 t=450.2s r (m) 012345 u θ (m/s) 0 0.2 0.4 0.6 t=500.2s Figure 5.6: TCS azimuthal-averaged u q velocity profiles for times t = 100 500s. Error bars represent the standard deviation of the sample velocities. u θ (m/s) 0 0.05 0.1 0.15 0.2 t=550.2s u θ (m/s) 0 0.05 0.1 0.15 0.2 t=600.2s u θ (m/s) 0 0.05 0.1 0.15 0.2 t=650.3s u θ (m/s) 0 0.05 0.1 0.15 0.2 t=700.3s u θ (m/s) 0 0.05 0.1 0.15 0.2 t=750.3s u θ (m/s) 0 0.05 0.1 0.15 0.2 t=800.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 0.15 0.2 t=850.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 0.15 0.2 t=900.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 0.15 0.2 t=950.4s Figure 5.7: TCS azimuthal-averaged u q velocity profiles for times t = 550 950s. Error bars represent the standard deviation of the sample velocities. 122 u r (m/s) -0.1 -0.05 0 0.05 0.1 t=100.1s u r (m/s) -0.1 -0.05 0 0.05 0.1 t=150.1s u r (m/s) -0.1 -0.05 0 0.05 0.1 t=200.1s u r (m/s) -0.1 -0.05 0 0.05 0.1 t=250.1s u r (m/s) -0.1 -0.05 0 0.05 0.1 t=300.1s u r (m/s) -0.1 -0.05 0 0.05 0.1 t=350.2s r (m) 012345 u r (m/s) -0.1 -0.05 0 0.05 0.1 t=400.2s r (m) 012345 u r (m/s) -0.1 -0.05 0 0.05 0.1 t=450.2s r (m) 012345 u r (m/s) -0.1 -0.05 0 0.05 0.1 t=500.2s Figure 5.8: TCS azimuthal-averaged u r velocity profiles for times t = 100 500s. Error bars represent the standard deviation of the sample velocities. u r (m/s) -0.04 -0.02 0 0.02 0.04 t=550.2s u r (m/s) -0.04 -0.02 0 0.02 0.04 t=600.2s u r (m/s) -0.04 -0.02 0 0.02 0.04 t=650.3s u r (m/s) -0.04 -0.02 0 0.02 0.04 t=700.3s u r (m/s) -0.04 -0.02 0 0.02 0.04 t=750.3s u r (m/s) -0.04 -0.02 0 0.02 0.04 t=800.3s r (m) 0246 u r (m/s) -0.04 -0.02 0 0.02 0.04 t=850.3s r (m) 0246 u r (m/s) -0.04 -0.02 0 0.02 0.04 t=900.3s r (m) 0246 u r (m/s) -0.04 -0.02 0 0.02 0.04 t=950.4s Figure 5.9: TCS azimuthal-averaged u r velocity profiles for times t = 550 950s. Error bars represent the standard deviation of the sample velocities. 123 0 0.5 1 1.5 ω z (s -1 ) t=100.1s 0 0.5 1 ω z (s -1 ) t=150.1s 0 0.5 1 ω z (s -1 ) t=200.1s 0 0.2 0.4 0.6 0.8 ω z (s -1 ) t=250.1s 0 0.2 0.4 0.6 0.8 ω z (s -1 ) t=300.1s 0 0.2 0.4 0.6 ω z (s -1 ) t=350.2s 024 r (m) 0 0.2 0.4 0.6 ω z (s -1 ) t=400.2s 024 r (m) 0 0.2 0.4 0.6 ω z (s -1 ) t=450.2s 024 r (m) 0 0.2 0.4 ω z (s -1 ) t=500.2s Figure 5.10: TCS azimuthal-averaged vorticity profiles for times t= 100 500s. Error bars repre- sent the standard deviation of the sample velocities. 0 0.1 0.2 0.3 0.4 ω z (s -1 ) t=550.2s 0 0.1 0.2 0.3 0.4 ω z (s -1 ) t=600.2s 0 0.1 0.2 0.3 0.4 ω z (s -1 ) t=650.3s 0 0.1 0.2 0.3 ω z (s -1 ) t=700.3s 0 0.1 0.2 0.3 ω z (s -1 ) t=750.3s 0 0.1 0.2 0.3 ω z (s -1 ) t=800.3s 0246 r (m) 0 0.1 0.2 0.3 ω z (s -1 ) t=850.3s 0246 r (m) 0 0.1 0.2 0.3 ω z (s -1 ) t=900.3s 0246 r (m) 0 0.1 0.2 ω z (s -1 ) t=950.4s Figure 5.11: TCS azimuthal-averaged vorticity profiles for times t= 550 950s. Error bars repre- sent the standard deviation of the sample velocities. 124 0 0.2 0.4 0.6 u θ (m/s) t=100.1s RMSE=0.051 R=0.94m, α=0.36 0 0.2 0.4 u θ (m/s) t=150.1s RMSE=0.035 R=0.84m, α=0.32 0 0.1 0.2 0.3 0.4 u θ (m/s) t=200.1s RMSE=0.026 R=0.92m, α=0.31 0 0.1 0.2 0.3 u θ (m/s) t=250.1s RMSE=0.026 R=0.92m, α=0.34 0 0.1 0.2 0.3 u θ (m/s) t=300.1s RMSE=0.018 R=0.94m, α=0.34 0 0.1 0.2 u θ (m/s) t=350.2s RMSE=0.021 R=0.97m, α=0.35 012345 r (m) 0 0.1 0.2 u θ (m/s) t=400.2s RMSE=0.016 R=1.01m, α=0.35 012345 r (m) 0 0.05 0.1 0.15 0.2 u θ (m/s) t=450.2s RMSE=0.011 R=1.08m, α=0.35 012345 r (m) 0 0.05 0.1 0.15 0.2 u θ (m/s) t=500.2s RMSE=0.010 R=1.15m, α=0.39 Figure 5.12: a-profile fit for times t = 100 500s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. 0 0.05 0.1 0.15 u θ (m/s) t=550.2s RMSE=0.010 R=1.21m, α=0.39 0 0.05 0.1 0.15 0.2 u θ (m/s) t=600.2s RMSE=0.011 R=1.27m, α=0.43 0 0.05 0.1 0.15 u θ (m/s) t=650.3s RMSE=0.009 R=1.35m, α=0.41 0 0.05 0.1 0.15 u θ (m/s) t=700.3s RMSE=0.011 R=1.41m, α=0.40 0 0.05 0.1 u θ (m/s) t=750.3s RMSE=0.008 R=1.48m, α=0.37 0 0.05 0.1 u θ (m/s) t=800.3s RMSE=0.008 R=1.45m, α=0.37 0246 r (m) 0 0.05 0.1 u θ (m/s) t=850.3s RMSE=0.008 R=1.48m, α=0.35 0246 r (m) 0 0.05 0.1 u θ (m/s) t=900.3s RMSE=0.008 R=1.55m, α=0.34 0246 r (m) 0 0.05 0.1 u θ (m/s) t=950.4s RMSE=0.007 R=1.56m, α=0.33 Figure 5.13: a-profile fit for times t = 550 950s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. 125 0 0.05 0.1 u θ (m/s) t=1000.4s RMSE=0.008 R=1.66m, α=0.34 0 0.05 0.1 u θ (m/s) t=1050.4s RMSE=0.006 R=1.69m, α=0.34 0 0.05 0.1 u θ (m/s) t=1100.4s RMSE=0.008 R=1.66m, α=0.34 0 0.05 0.1 u θ (m/s) t=1150.4s RMSE=0.008 R=1.77m, α=0.30 0 0.02 0.04 0.06 0.08 u θ (m/s) t=1200.4s RMSE=0.007 R=1.72m, α=0.34 0 0.05 0.1 u θ (m/s) t=1250.5s RMSE=0.006 R=1.74m, α=0.34 0246 r (m) 0 0.02 0.04 0.06 0.08 u θ (m/s) t=1300.5s RMSE=0.006 R=1.81m, α=0.35 0246 r (m) 0 0.02 0.04 0.06 0.08 u θ (m/s) t=1350.5s RMSE=0.005 R=1.85m, α=0.35 0246 r (m) 0 0.02 0.04 0.06 0.08 u θ (m/s) t=1400.5s RMSE=0.005 R=1.86m, α=0.36 Figure 5.14: a-profile fit for times t = 1000 1400s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. 0 0.02 0.04 0.06 u θ (m/s) t=1450.5s RMSE=0.005 R=1.92m, α=0.36 0 0.02 0.04 0.06 u θ (m/s) t=1500.5s RMSE=0.005 R=2.02m, α=0.35 0 0.02 0.04 0.06 u θ (m/s) t=1550.6s RMSE=0.004 R=2.00m, α=0.36 0 0.02 0.04 0.06 u θ (m/s) t=1600.6s RMSE=0.004 R=2.04m, α=0.34 0 0.02 0.04 0.06 u θ (m/s) t=1650.6s RMSE=0.004 R=2.11m, α=0.34 0 0.02 0.04 0.06 u θ (m/s) t=1700.6s RMSE=0.004 R=2.17m, α=0.34 02468 r (m) 0 0.02 0.04 0.06 u θ (m/s) t=1750.6s RMSE=0.004 R=2.20m, α=0.35 02468 r (m) 0 0.02 0.04 0.06 u θ (m/s) t=1800.6s RMSE=0.004 R=2.34m, α=0.33 02468 r (m) 0 0.02 0.04 0.06 u θ (m/s) t=1850.7s RMSE=0.004 R=2.31m, α=0.34 Figure 5.15: a-profile fit for times t = 1450 1850s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. 126 0 0.2 0.4 0.6 u θ (m/s) t=100.1s RMSE=0.051 R=0.35m, α=0.33 0 0.2 0.4 u θ (m/s) t=150.1s RMSE=0.035 R=0.37m, α=0.41 0 0.1 0.2 0.3 0.4 u θ (m/s) t=200.1s RMSE=0.025 R=0.40m, α=0.42 0 0.1 0.2 0.3 u θ (m/s) t=250.1s RMSE=0.024 R=0.39m, α=0.43 0 0.1 0.2 0.3 u θ (m/s) t=300.1s RMSE=0.017 R=0.36m, α=0.38 0 0.1 0.2 u θ (m/s) t=350.2s RMSE=0.021 R=0.36m, α=0.38 012345 r (m) 0 0.1 0.2 u θ (m/s) t=400.2s RMSE=0.017 R=0.36m, α=0.37 012345 r (m) 0 0.05 0.1 0.15 0.2 u θ (m/s) t=450.2s RMSE=0.011 R=0.35m, α=0.33 012345 r (m) 0 0.05 0.1 0.15 0.2 u θ (m/s) t=500.2s RMSE=0.012 R=0.35m, α=0.33 Figure 5.16: Whittaker-profile fit for times t = 100 500s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. u θ (m/s) 0 0.05 0.1 0.15 RMSE=0.011 R=0.38m, α=0.35 t=550.2s u θ (m/s) 0 0.05 0.1 0.15 0.2 RMSE=0.013 R=0.42m, α=0.40 t=600.2s u θ (m/s) 0 0.05 0.1 0.15 RMSE=0.011 R=0.40m, α=0.33 t=650.3s u θ (m/s) 0 0.05 0.1 0.15 RMSE=0.012 R=0.39m, α=0.29 t=700.3s u θ (m/s) 0 0.05 0.1 RMSE=0.009 R=0.37m, α=0.24 t=750.3s u θ (m/s) 0 0.05 0.1 RMSE=0.009 R=0.39m, α=0.25 t=800.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 RMSE=0.008 R=0.39m, α=0.25 t=850.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 RMSE=0.009 R=0.41m, α=0.25 t=900.3s r (m) 0246 u θ (m/s) 0 0.05 0.1 RMSE=0.008 R=0.44m, α=0.27 t=950.4s Figure 5.17: Whittaker-profile fit for times t = 550 950s. The root-mean-square-error (RMSE), length scale R and steepness parameter a for each instant are noted. 127 The two terms in the RHS correspond to the turbulent vorticity diffusion (left) and the turbulent bottom friction (right), which are responsible for the vorticity decay. Using appropriate hori- zontal (L), vertical (H), vorticity (W) and velocity (UWL) scales, Seol and Jirka (2010) have shown that the order of magnitude of the vorticity diffusion and turbulent bottom friction terms are O 0:1 p c f =2(H=L)W 2 and O (c f =2)(L=H)W 2 , respectively. The friction coefficient ranges around c f 0:01 for laboratory conditions (Socolofsky and Jirka, 2004), and the orders of mag- nitude become O(0:0071(H=L)W 2 ) and O(0:005(L=H)W 2 ) ,respectively. Therefore, for typical shallow flows for which L >> H, the bottom friction term dominates the vortex spin-down and vorticity decay. 5.5.1 Vortex decay model Following the above dimensional analysis, it is expected that the bottom friction will dominate the TCS energy decay. This section provides a simple force-balance model to match the experimental energy decay curves. The force F exerted on the fluid inside a cylinder of radius r and length H, by the fluid sur- rounding it is given by (Batchelor, 1967) as F = 2pr 2 Hm ¶u q ¶r u q r = 2pr 2 Hs; (5.24) wheres(r)=m ¶u q ¶r u q r is the tangential stress on a fluid element on the surface of a cylinder (2pr) at radius r. Equating the rate of change of angular momentum of the fluid in the cylinder shell to the force acting on it gives (Batchelor, 1967) ¶(2prr 2 Hu q ) ¶t = ¶ ¶r 2pmr 2 ¶u q ¶r u q r ; (5.25) which is the equivalent of the radial diffusion equation (Eqn. 5.10). Simplifying Eqn. 5.25 leads to H ¶(r 2 u q ) ¶t = 1 r ¶ ¶r r 2 m ¶u q ¶r u q r = 1 r ¶ ¶r r 2 s : (5.26) 128 Equating the tangential stresss with the bottom shear stress, leads to the force balance equation for a rotating fluid element at a distance r, in a vortex of constant radius R max (Seol and Jirka, 2010) d dt (rV H)= rt b r ; (5.27) where V = u q (R max ) is the azimuthal velocity at r= R max . The bottom shear stresst b for turbulent is computed from Eqn. 5.5 for a purely azimuthal flow ( ¯ u; p ¯ u 2 + ¯ v 2 = V ) as t b =r c f 2 V 2 : (5.28) For laminar conditions, the bottom friction is modeled using the Rayleigh friction parameter given in Eqn. 5.7. Substituting for c f and integrating Equation 5.27, leads to V = 1 1 V 0 + c f 2H t ; (5.29) for the turbulent case and V = V 0 exp 3n H 2 t (5.30) for laminar conditions, where V 0 is the initial azimuthal velocity. Fitting Equation 5.29 to the azimuthal decay data results to a best-fitting c f = 0:01 with 0.99 correlation coefficient and root-mean-square-error of 9 10 4 m=s (Figure 5.18a). The first-order model captures the physical process of TCS energy decay due to bottom friction very well. The reason for the very good fit of the simple model is that the vortex core does not expand until t 390s and therefore R max (the radial distance to u q;max ) remains approximately constant. The decay of the minimum radial velocity u r is also plotted with blue circles in Figure 5.18a up to t= 800s. Although the decay curve shows significant scatter (compared to the u q;max decay) due to the high uncertainty in the estimation of the u r profiles, it follows a clear trend with a much steeper slope than u q;max . 129 Assuming axisymmetry, vorticity about the zaxis can be approximated as w z (r) 1 r ¶ ¶r (ru q )= u q r + ¶u q ¶r : (5.31) Evaluating vorticity at radius R max corresponding to u q;max , leads to w z (R max )= u q;max R max = V R max => V = R max w z (R max ): (5.32) Replacing V for R max w in Eqn. 5.29 gives w = 1 1 w 0 + c f R max 2H t ; (5.33) and therefore the same methodology is applicable for the vorticity decay curve. The fit shown in Figure 5.18b, yields c f R max 0:004m. The fit is not as good as for the azimuthal velocity (RMSE of 0:022s 1 ), since R max varies after t 390s. V ortex radial growth is the subject of the next Section. 10 −3 10 −2 10 −1 10 0 u θmax , −u rmin (m/s) 0 500 1000 1500 2000 2500 3000 3500 time (s) (a) 10 −2 10 −1 10 0 10 1 ω max (1/s) 0 500 1000 1500 2000 2500 3000 3500 time (s) (b) Figure 5.18: Decay of TCS mean properties. (a) Maximum azimuthal (black line) and minimum radial (blue circles) velocity decay. (b) maximum vorticity decay. The 1 st order decay prediction is plotted with the dashed red lines. 130 5.5.2 TCS growth Initially, the TCS undergoes internal oscillations due to high turbulence while later the TCS radius grows with time due to lateral ambient fluid entrainment. Merging of separate structures also leads to the formation larger structures over time (Mcwilliams, 1990; Dracos et al., 1992). The radial growth of the TCS due to ambient fluid entrainment (viscous diffusion) can be found through a particular solution of the vorticity Eqn. 5.11, given by (Mcwilliams, 1990) w z (r;t)=w 0 t 0 t exp r 2 4nt : (5.34) Taking the derivative with respect to time and solving for dw z =dt = 0 leads to a bulk radius r growth of R(t) p 4nt, or R vortex (t)= q R 2 0 + 4Et; (5.35) where R 0 is the vortex radius at t= 0 and E is the diffusivity. The diffusivity E is equal to viscosity n for laminar conditions and for turbulent conditions be approximated by (Seol and Jirka, 2010) En turb ; (5.36) wheren turb is the added viscosity due to turbulence given in Eqn. 5.4. The vortex radius, as discussed in Section 4.3.2, can be defined both in terms of the radial distance to the maximum azimuthal velocity (through the swirl strength maps) but also in terms of the vorticity decay and profile expansion. The former corresponds to the vortex core radius and the latter to the vortex structure boundary, or bulk radius. In this work, the two radii are evaluated through the azimuthal-averaged vorticity and velocity profiles. The vorticity-defined radius correspond to the radius R vortex at which w(r;t)< 0:02w max (t). Two approaches are used to calculate the vortex core growth from the u q profiles: i) One approach is to evaluate R max from the azimuthal-averaged profiles (using the procedure given in Section 3.10.2). The overlap of the interrogation annulus (of width W R = 0:5m) used for azimuthal-averaging has been increased to get a profile resolution ofDr= 0:1m. ii) The other approach is to obtain R max from the fitted u q “injection” a-profiles (e.g. Eqn. 5.18). 131 The evaluated vortex radii using the above methods are plotted in Figure 5.19a. The reference vortex growth growth rates of the forme p t are plotted for different e and the diffusivity E is computed for c f = 0:01 and V 0 = 0:6m=s. 10 2 10 3 time (s) 10 0 10 1 R vortex (m) t∼ 390s (a) √ t √ 4Et √ 2Et √ t/15 0 1000 2000 3000 time (s) 0.5 1 1.5 2 2.5 Re H ×10 6 (b) Figure 5.19: (a) TCS size growth with time. Dashed lines correspond to viscous diffusion rates as functions of p t, red crosses and green squares to the TCS-core growth (from the azimuthal- averaged u q profiles and the u q profile fits), black circles to TCS bulk radius growth, and the cyan curve to the mean-minimum distance to the closest vertical boundary. (b) Local depth-based Reynolds number decay (Eqn. 5.38). The two methods of calculating the vortex core radius show good agreement. The azimuthal- averaged u q;max is more scattered due to the uncertainty in the evaluation sample. On the other hand, the radial growth using the fitted curves is more stable and fits through the azimuthal- averaged data-points. Both growth curves show that the vortex core radius remains constant until t 390s, then follows the p t growth rate. The transition time between the two regimes also corre- sponds to the time when the secondary flow diminishes and the vortex profile becomes Gaussian. The vortex bulk radius on the other hand grows p t, since the start of the measurements at t= 100s. More specifically, the growth rate follows p 4Et in the time period t 100230s. After t 230s, vortex growth is dictated by the distance of the TCS-center to the vertical boundaries and more or less follows that curve. The TCS is constantly adjusting its position to allow it to preserve its momentum and grow in size by entraining ambient fluid. As the TCS size grows, it 132 eventually positions itself at the point that fits the maximum TCS size. An algorithm that seeks to maximize the minimum distance to the boundaries (keeping in mind that the wave maker is fully retracted, i.e. at x=1 m), gives the maximum circle that fits the polygon. It has a radius equal to r max = 8:5 m and center coordinates (x c ;y c )=(7:46;4:80) (Figure 5.20). The centers of the coherent structures of the trials extending to t 3000s almost reached their final position in the offshore basin. The maximum radius the coherent structure can theoretically sustain, without the spatial growth constraint of the vertical boundaries, can be computed using R max 2H c f = 110m: (5.37) It corresponds to the radius at which the vortex will loose all its rotational energy before making one full rotation (Jirka, 2001). Figure 5.20: TCS center paths from all trials. The mean path is shown with the bold curve. The dashed line corresponds to the maximum circle fitting the polygon defined by the solid boundaries. 133 Defining a local depth-based Reynolds number using the vortex size (Seol and Jirka, 2010) Re H (t)= w max (t)R vortex (t)H n ; (5.38) allows to characterize the state of the flow- an approximate Reynolds number for laminar condi- tions corresponds to Re H 500 (Seol and Jirka, 2010). Figure 5.19b shows the turbulence decay with time and confirms that the flow remains fully-turbulent throughout the duration of the experi- ment. 5.6 Flow transition to Q2D 5.6.1 Theoretical background 2D turbulence is a well studied phenomenon first predicted by Kraichnan (1967), and is related to the spectral turbulent kinetic energy transfer. Whereas in 3D turbulence the energy is transferred in smaller scales and eventually the kinetic energy is absorbed through molecular viscosity, in what is called the enstrophy cascade, in 2D turbulence the coupled constraints of energy and enstrophy conservation lead to the inverse cascade from smaller to larger scales (Boffetta et al., 2000). When forced at a fixed scale, 2D turbulence should exhibit a k 5=3 energy spectrum extending from the source to larger scales, where k is the wavenumber (Paret and Tabeling, 1997). Experimental evidence shows that the inverse energy transfer to larger scales is a result of clustering of small- scale, equal-signed vortices (Boffetta et al., 2000). In enstrophy cascade on the other hand, the flux of rotational momentum is constant among all eddy scales leading to a k 3 distribution of turbulent kinetic energy (Jirka, 2001). Thus, the mechanism of vortex stretching is not applicable for 2D turbulent flows. As a result of the inverse energy cascade, when shallow geophysical flows with large length to depth ratios (L>> H) are subjected to localized or distributed disturbances, the internal oscil- lations grow into large-scale coherent structures (Jirka, 2001). From experimental observations on shallow jets, Dracos et al. (1992) identified three stages of development with respect to the 134 distance x from the source: i) the “near-field” (x=H 0:1), in which the flow contains highly 3D- scale turbulence, ii) the “middle-field”, where a mean 2D flow has developed which interacts with 3D turbulent flow, and iii) the “far-field” region (x=H 10), where the structure of the flow has grown into scales greater than the depth and the flow has a 2D character (Jirka, 2001). There- fore, the inverse energy cascade becomes effective soon after the initiation of the disturbance, but convective distance is needed for the flow to re-organize and become 2D. In two-dimensional (2D) flows, the motion is confined on a plane. This idealized flow is rarely found in nature either due to 3D instabilities or because their two-dimensionality is dependent on external factors that affect their stability, most commonly in the presence of solid boundaries (Dolzhanskii et al., 1992). The no-slip boundary conditions and the resulting boundary layer give rise to transverse velocities that break down the two-dimensionality of the flow. The flow com- ponents which are not part to the main flow structure constitute the secondary flow. The effect of bottom friction can be minimized in experiments by the use of a two-layer stratification (e.g. Paret and Tabeling (1997)) to inhibit vertical motions reaching the upper fluid layer (Akkermans et al., 2008). Under sufficient vertical confinement, apart from the bottom layer in the water column, the flow still behaves in a two-dimensional fashion, and the term quasi-two-dimensional (Q2D) is used to describe it. Dolzhanskii et al. (1992) proposed that the two-dimensionality of the flow depends on two dimensionless parameters: the traditional Reynolds number Re n = UL=n and the Reynolds num- ber in terms of the external friction Re l = U=lL, where U and L and typical velocity and length scales andl = 2n=H 2 is the external friction parameter (equivalent to the Rayleigh friction param- eter defined in Section 5.1). It follows that for Re l << Re n there is no dependence on Re n and the lateral diffusion term Ñ 2 u in Eqn. 5.8 can be ignored. Dolzhanskii et al. (1992) further argue that the relative importance of the two terms Re l << Re n is a necessary condition for a quasi-two-dimensional flow state (taking the ratio Re n =Re l , the inequality condition translates to H= p 2 << L). In other words, any flow (and resulting vertical diffusion) that can be modeled by parameterizing the bottom shear stress as an external force, i.e. a flow dominated by bottom friction effects, can be characterized as Q2D. 135 Satijn et al. (2001) studied analytically and numerically the two-dimensional structure of monopolar vortices in shallow fluid layers. Following Dolzhanskii et al. (1992), they parametrize the flow conditions with two Reynolds numbers: the traditional Re n = UL=n and the external friction Re l = U=lL, where the external friction is defined as l =p 2 n=4H 2 . By varying these two parameters it is possible to examine both the inertial effects and the role of shallowness in the two-dimensionality of the flow. In order to quantitatively characterize the flow as Q2D, they defined two conditions. The fist is related to the ratio of the kinetic energy of the secondary flow to the kinetic energy of the primary flow. The kinetic energy (in each direction i=(r;q;z), in polar coordinates) contained in a circular vortex of radius R and local depth H is defined as E k;i = 2p Z H 0 Z R 0 1 2 r(z)u 2 i (r;z)rdrdz: (5.39) Satijn et al. (2001) defined the base flow by an axisymmetric vortex with an isolated Gaussian azimuthal velocity profile. The radial u r and vertical u z velocities that develop due to the boundary conditions constitute the secondary flow components. Satijn et al. (2001) found that the flow behaves Q2D if the kinetic energy in the radial and vertical velocity components is insignificant compared to the primary flow. This is quantitatively represented through two inequalities: q r (t)= E k;r (t) E k;q (t) 0:01; q z (t)= E k;z (t) E k;q (t) 0:01: (5.40) The second condition is related to the shape of the vortex vorticity profile. It was found by Satijn et al. (2001) that flows do not behave in a Q2D fashion if the secondary flow deforms the vorticity profile from its base state. This is quantitatively formulated as Q= R R 0 [w 0 z (r 0 ;0)w 0 z (r 0 ;t)] 2 rdr R R 0 w 02 z (r 0 ;0)rdr 0:1; (5.41) wherew 0 z (r 0 ;t) is the re-scaled vorticity profile with respect to its radius and amplitude. Whereas the condition given in Eqn. 5.40 requires integration (and consequently information) over the water column, the second requirement can be applied on information available only on the water surface. 136 5.6.2 Experimental observations and secondary flow quantification From the assimilated data and the analysis presented in this chapter, the main observations can be summarized as: The primary flow corresponds to a stirring-type profile of the form given in Eqn. 5.18, with the steepness factor a ranging between a= 0:33 0:43, best described with a= 0:34 after t= 1000s. In the period t = 100 400s, the azimuthal velocity profiles deviate from the base flow profile and are best-matched using a similarity profile valid for the fplane (Eqn. 5.21). The radial velocity profile is characterized by a zone of flow convergence near the TCS center with radius R convergence 1:65m, followed by a zone of flow divergence. The flow convergence zone is associated with downwelling at the TCS center, where pressure also reaches a minimum. The low pressure point is coupled with a strong surface elevation depression. Around the TCS center, zones of upwelling can be identified from 2D flow divergence maps, providing evidence that the flow undergoes vertical re-circulation about the horizontal axis. These 3D flow features have also been observed in propagating shallow water dipoles (e.g. Akkermans et al. (2008)). The decay of the mean properties due to bottom friction, i.e. azimuthal velocity and vorticity, starts before t = 100s and the decay rate is well described using a first-order decay model. The minimum value of the radial velocity component decays faster than the azimuthal, becoming two orders of magnitude smaller by t 800s. The vortex bulk radius expands due to viscous diffusion before t= 100s, whereas the vortex core radius is not growing until t 390s. This timing coincides with the core transition time from convergent to weakly divergent. This finding is confirmed by the visual observa- tion that the TCS-center tracer-conglomerate remains compact until t 320s and then starts expanding and breaking-up into fractions (Figure 3.20). 137 Through the above observations, it can be postulated that there is a direct link between the decay of the radial velocity component and the vortex core radius. When the secondary flow becomes insignificant, the vortex core growth eventually allows the azimuthal velocity profile of the TCS to adjust to the a-profile. An attempt is made here to quantify the relative magnitude of the secondary flow following the methodology proposed by Satijn et al. (2001). Since experimental velocity data over the water column are not available, Eqn. 5.40 is applied only on the water surface, so that the kinetic energy is computed as E k;i =pr Z R 0 u 2 i (r)rdr; (5.42) per unit depth (in units of Nm=m). The integral is evaluated using the azimuthal-averaged profiles described in Section 5.4 and the radius of integration R is kept constant at R= 4m to exclude the effect of the vortex radial growth in this analysis. The resulting kinetic energy decay, evaluated every 10 frames (10=29:97s), is plotted in Fig- ure 5.21. The azimuthal velocity kinetic energy E q exhibits a smooth decay, consistent with the decay of the maximum azimuthal velocity due to bottom friction (Section 5.5.1). The radial veloc- ity kinetic energy E r , as expected from the observations summary of this section, exhibits a rapid decay at t 390s, and remains insignificant thereafter. Both results further reinforce the afore- mentioned observations. The E r decay curve experiences strong fluctuations that require further investigation. Large scattering was expected since the radial velocity profiles are subject to large uncertainty. However, the signal is not random but rather appears to be periodic. The period of the fluctuations is inferred by filtering the decay curve using a moving average scheme (with a window of 20s) and subtracting the raw signal from the mean (Figure 5.21a). Performing a Fast Fourier Transformation on the residual reveals the dominant periods of the fluctuations. From the spectral energy curve shown in Figure 5.21b, three distinct peaks can be identified. In the absence of compelling evidence about the origin of these periodic signals, answers were sought in the basin response to the wave forcing. The resonant periods and modes of the basin were computed from the surface elevation timeseries following the methodology outlined in the work by Maravelakis et al. (2014) (Figures 5.22& 5.23). 138 The period of the first and second peaks in the spectral energy plot (T = 8:05s and T = 16:7s) match with two identified resonant periods, and the period of the third peaks is closely matched with the resonant period of T = 23:1s. However, at this stage of the analysis, until all other sources of uncertainty are ruled out, the origin of the E r decay fluctuation remains speculative. Figure 5.21d shows the decay of the ratio E r =E q using the smoothed E r decay curve. The ratio at all times remains below the E r =E q = 0:01 transition threshold suggested by Satijn et al. (2001). However, since the profiles are not integrated over the depth of the water column, the suggested threshold is not directly applicable. Moreover, whether the trend of the filtered decay curve is physical remains inconclusive until further investigation. It is nevertheless safe to conclude that since the radial velocity associated with the three-dimensional motions becomes negligible at t 390s, it can be considered as the transition point of the flow to quasi-two-dimensional. 139 0 200 400 600 800 1000 time (s) 0 10 20 30 40 50 E r (N*m/m) (a) 0 200 400 600 800 1000 time (s) 0 1000 2000 3000 4000 5000 6000 7000 E θ (N*m/m) (c) 0 20 40 60 T (s) 0 50 100 150 200 SE (N 2 s) T=8.05s T=16.7s T=22.5s (b) 0 200 400 600 800 1000 time (s) 0 1 2 3 4 5 6 E r /E θ ×10 -3 (d) 100 200 300 400 time (s) -20 0 20 residual Figure 5.21: Quantification of the secondary flow magnitude. (a) Raw (black) and smoothed (dashed red) kinetic energy decay of the radial velocity. (b) Spectral energy plot of the radial velocity fluctuations. (c) Kinetic energy decay of the azimuthal velocity. (d) Kinetic energy ratio E r =E q decay. 140 T=78.8s T=27.5s T=23.1s T=20.6s T=16.7s T=14.2s T=12.7s T=11.6s Figure 5.22: Basin resonant periods and modes. 141 T=10.9s T=9.9s T=9.4s T=8.0s T=7.4s T=6.7s T=6.3s T=6.0s Figure 5.23: Basin resonant periods and modes - continued from Figure 5.22. 142 Chapter 6 Field measurements of vortical motions in Ventura harbor during the 2015 Chilean tsunami Tsunami–induced coastal currents are spectacular manifestations of nonlinear and chaotic phenom- ena, even during marginal or otherwise unremarkable events. Due to their long–periods, tsunamis transport significant energy near shore, and shear and turbulent features appear through their inter- action with the ubiquitous irregularity of coastal bathymetry. The oscillatory character of tsunami wave trains lead to numerous flow reversals, spawning persistent turbulent coherent structures (e.g. large “whirlpools”) which can dominate damage and transport potential. However, no quantitative measurements exist which might provide physical insight into this turbulent variability, and not a single video recording to help elucidate how these vortical structures evolve and terminate. We report our measurements of currents in Ventura Harbour, California, generated by the 2015 Chilean M8.3 earthquake. We measured Lagrangian surface–velocities across a large spatial area, as the event unfolded. From this map of the tsunami flow field, we find that a tsunami with a near-shore amplitude of 30 cm at 6 m depth produced unexpectedly large currents of up to 1.5 m/s. Coherent turbulent features across a wide range of scales appear throughout the event, often generating the greatest local currents. 143 On the 16th of September, 2015 a M8.3 (USGS) earthquake struck the central coast of Chile, with its epicentre between Coquimbo and Valparaiso, 490 km north of the 2010 M8.8 Maule earth- quake. The thrust-faulting earthquake ruptured along the subduction zone interface between the Nazca and South America plates at a relatively shallow focal depth. The shaking lasted for three minutes, a harbinger warning for evacuation in near-field coastal areas, where the waves arrived within minutes. The local emergency management was effective in disseminating information, thus no casualties from the tsunami were reported, despite the extensive flooding in hard-hit coastal areas like Coquimbo and Tongoy. The more energetic 2010 Maule tsunami, with run-up over 29 m in Constitucion, claimed 124 lives (Fritz et al., 2011). In the far-field, the West Coast and Alaska Tsunami Warning Center issued an advisory for the SW California coastline, which suggested that tsunami-induced flooding was not expected to be significant, but potentially hazardous currents were possible. Harbour closures were left at the discretion of local harbour masters, with just one of the five coastal counties under the Advisory opting for a temporary closure of all its beaches and harbours. Ventura Harbour remained opera- tional, yet experienced the highest amplitude waves in California, with the tide gauge registering a peak of 32 cm above tide. Nano-tsunamis of this size do not attract attention, because their effects along open coastlines are minimal. Motivated by visual observations and numerical simulations of complex currents in harbours (e.g. Okal et al. (2006a); Lynett et al. (2012)), the arrival of this low-amplitude tsunami in California was an opportunity to collect quantitative flow measurements with minimal risk, so as to help elucidate complex near-shore currents. Due to the focus of most existing tsunami studies on overland flooding depths, and combined with the inherent difficulty of measuring tsunami currents, current magnitudes in ports have relied on eyewitness accounts (Borrero et al., 2015). A small but increasing number of current measurements have been reported from acoustic current meters (Lacy et al., 2012; Borrero et al., 2013; Admire et al., 2014) and by tracking floating objects in eyewitness videos (Admire et al., 2014; Fritz et al., 2012). In contrast, here, we provide a comprehensive description of the flow field, based on a quantitative data of tsunami current effects and large scale coherent flow structures unfolding during a real event. 144 In coordination with Ventura Harbour Patrol, we set up our instrumentation at 5:00 AM local time, 1/2 hr before the arrival of the leading wave. An HD camera was positioned to continuously observe the harbour entrance channel (Figure 6.1) and did so for five hours following the arrival of the first wave. To optically track objects in the flow, we deployed surface tracers, motivated by results from shallow flows in rivers, known to be dominated by two-dimensional structures. We thus infer that surface velocities are reasonable for describing many aspects of the flow (Weitbrecht et al., 2002). Our tracers included 60 cm diameter rubber balls and custom-made 40x40x3 cm styrofoam floaters. Both were large-enough to be distinguishable with distances up to 600 m from the camera. All styrofoam floats carried flashing lights to locate them in the dark, while one carried an InReach Explorer (DEX) GPS unit. This GPS device bears a non-differential receiver that allows for position logging every 2 sec, adequate to resolve the motions of interest. After extracting the tracer’s position time-history (see Methods), the coordinates and velocities were referenced to one of the two local coordinate systems shown in Figure 6.1b. The choice of local coordinate system depended on which branch each tracer would follow at the channel bifurcation. The velocity pairs referenced to the local coordinate system correspond to the channel- parallel (u y ) and channel-normal velocities (u x ). The compiled flow time-histories allowed us also to evaluate the accuracy of a standard numer- ical model used for tsunami simulations. Given the complexity of tsunami near–shore flows, one might surmise that high–resolution and turbulence–resolving models are required to adequately describe the flow. However, uncertainties and imprecision in field-scale simulation can often obscure any such increase in modelled physics. This posed the interesting and fundamental ques- tion as to whether a relatively simple hydrodynamic model might yield useful predictions. We thus used MOST (Method Of Splitting Tsunamis), which solves the nonlinear shallow water equations using a finite-difference scheme(Titov and Synolakis, 1997). MOST has been benchmarked with laboratory results and extensively validated operationally for tsunami events, but only to predict maximum inland penetration. In California, MOST has been used to produce the state’s tsunami inundation maps (Barberopoulou et al., 2011b,a) and to predict currents (Ayca et al., 2012). In a recent study (Lynett et al., 2014), current speeds in Crescent City, CA, during the 2011 Japan 145 Figure 6.1: (a) The camera’s FOV and the ground control points used to obtain the camera’s extrin- sic parameters are shown in the inlet; (b) The Ventura entrance channel, camera position and field of view (FOV), and local coordinate systems (labelled as CS 1 and CS 2). The coordinate system y axes are positioned in the centre-line of the entrance channel and its two branches, whereas the x axes define the distance normal to the centre-line. Three zones are defined by distance x: zone 1 withjxj<= 25m, zone 2 with 50m<jxj<= 25m, and zone 3 withjxj> 75m. 146 tsunami were simulated using MOST and the comparison to estimated flow speeds at the inner harbour entrance (Admire et al., 2014) was shown to be satisfactory. For our purpose, we referenced the output velocities of MOST, in the longitudinal and latitu- dinal directions (see Methods for more information on the model setup), to the local coordinate systems defined in Figure 6.1b, keeping only the along-channel velocity component u y . Instead of comparing the model results to the measurements in a Lagrangian framework, we opted to repre- sent the model results through a mean channel flow speed by averaging u y over numerical points in the harbour channel with 185m< y< 385m. This representation allows for small timing inac- curacies and averages out localised flow speed variations in the model. Flow speed variations are represented through the extrema across the channel length, forming a flow–speed envelope. We summarise our quantitative results in Figure 6.2 which depicts the measured and predicted flow time histories across the harbour entrance channel (see Figure 6.1 for a top view of the tracer paths). The measurements show that the channel-parallel water flow is fluctuating with u y veloc- ities ranging between1:3 m/s. The velocity magnitude reached a maximum of 1.5 m/s at 6 m depth. These values are about four times greater than what simple scaling calculations (long-wave theory) would suggest for a 30 cm wave amplitude. The MOST channel-mean speed prediction fluctuates between0:5 m/s, and the min/max largely envelope the measurements. In general, the phasing matches the measurements well, with the exception of the fourth to sixth waves (6:50 to 7:40 AM local time). This disagreement is also evident in the tide gauge comparison and can be attributed to a range of factors, including inaccuracies in the offshore forcing (Lynett et al., 2014). The flow time history in Figure 6.2 has been separated in three zones, defined by the distance (x) normal to the local coordinate system centre–line. Zone 1, which corresponds to the middle section of the harbour channel, shows the smoothest flow pattern with smaller extrema. In contrast, flow in zones 2 & 3 experiences rapid acceleration/deceleration, as evident from the steep velocity gradients in the measurements. Steep velocity gradients appear when the flow in the channel starts to change direction, and vorticity is introduced from the boundaries. As examples, two tracers (labelled as ex 1 and ex 2) experiencing rapid deceleration are identified in Figure 6.2 and Figure 6.4. In ex 1, the tracer decelerates rapidly during the flow reversal and is caught in a 147 Figure 6.2: (a)-(d) Time history of channel-parallel velocity in Ventura Harbour, measured (points) versus predicted (curves). Crosses (+) represent GPS measurements, whereas circles (o) and squares () represent optical measurements referenced to coordinate systems CS 1 and CS 2 respectively. All GPS measurements are referenced to CS 1. The flow field has been separated in three zones, defined by the distance (x) normal to the channel center-line. The mean flow chan- nel speed u y predicted by MOST is overlaid in all subplots as continuous curves, whereas the dashed curves correspond to the minimum and maximum MOST-predicted values in each channel zone. All results are shown together in (d). Along-channel distance (y) shown with colormap. (e) Ventura’s tide gauge record compared to the surface elevation time history extracted from MOST. The tide gauge record has been high-pass filtered using a cutoff period of 100 min. The location of the tide gauge is shown in Figure 6.1. 148 Figure 6.3: Rectified camera frames captured at (a) 8:29 and (b) 8:43 AM local time, overlaid over ortho-rectified USGS aerial imagery (times are shown in Figure 6.2 with vertical dashed lines). (a) transitional-jet and TCS forming when flow reversed towards the harbour; (b) two transitional-jets expelled from the channel bifurcation corners. Vectors correspond to tracer flow speeds. TCS streamlines and the fronts of transitional-jets are approximately shown with dashed curves. rotational structure at the entrance of the harbour channel. The second example shows a particle accelerating towards the harbour; once it turns east at the bifurcation, it enters a zone of strong turbulence adjacent to the separation region on the lee side and decelerates rapidly. Interesting flow patters are observed after rectifying the image sequences (see Supplementary video of Kalligeris et al. (2016) which shows the evolution of the flow structures). Two frames captured during maximum flow speeds are shown in Figure 6.3. When the flow reversed towards the harbour (positive u y ), an anti–clockwise rotating turbulent coherent structure (TCS) can be identified at the bifurcation, followed by a transitional-jet dipole forming at the sharp seawall corner adjacent to the beach (Figure 6.3a). When the flow reversed towards offshore (negative u y ), transitional–jets were expelled from the bifurcation corners, carrying three tracers that experienced surface speeds up to 1 m/s (Figure 6.3b). This pattern appears similar to larger tidal flows through inlets, where starting-jet and expelled boundary-layer vortices form (Bryant et al., 2012) at the beginning of ebbing. 149 Our measurements of surface flows in Ventura Harbour’s entrance channel during the 2015 Chilean M8.3 earthquake show that the surface current speeds reached 1.5 m/s in a highly turbulent, but sub-critical flow. The vortical structures triggered by this nano tsunami persisted for over four hours, and surface amplitudes decayed slowly. The harbour channel boundaries fed the flow with vorticity, creating transitional–jets that grow in size and form coherent vortical structures (Jirka, 2001). Our findings suggest that forecasts of vortical motions should be considered in operational tsunami warnings, in addition to the now standard forecasts of maximal elevations. Had this tsunami been more energetic, such features would pose a substantial hazard for safe navigation in the channel and for birthed vessels. As it was, the tsunami effects in Ventura Harbour existed in a type of scientific “Goldilocks zone”, where the currents were strong enough to be clearly measured and of scientific relevance, but not strong enough to make navigation and field work in the water unsafe. 6.1 Methods 6.1.1 Optical particle tracking. We recorded HD video at 23.97 Hz when dark, then switched to time–lapse mode, capturing raw frames of 4928x3264 resolution every 10 sec. Extracting velocity vectors from the surface tracers identified in the image sequence required a series of steps. We calibrated the camera in the lab to obtain the intrinsic parameters (Bouguet, 2014) and removed lens distortion. We then applied the pinhole model, which relates image (p;q) to world (x;y;z) coordinates, in the form of the direct linear transformation equations 3.19. Rewriting equation 3.19 as a system of linear equations of the form A x= b, allowed us to obtain the least square solution of DLT coefficients using the ten non-coplanar ground control points (GCPs) for which world coordinates were (Holland et al., 1997). The world coordinates of the GCPs were obtained by a posterior survey (Figure 6.1a), and all elevations were referenced to the mean lowest low (MLLW) datum. 150 Once the DLT coefficients were known, the image coordinates of the surface tracers could be transformed to world coordinates by rewriting equation 3.19 as: 2 4 x y 3 5 = 2 4 [L 1 L 9 p] [L 2 L 10 p] [L 5 L 9 q] [L 6 L 10 q] 3 5 1 2 4 p L 4 + z t (L 11 p L 3 ) q L 8 + z t (L 11 q L 7 ) 3 5 (6.1) where z t is the tracer elevation, which needs to be independently estimated. z t was set equal to the tide level above MLLW without accounting for tsunami and higher frequency wave amplitudes. Due to the camera low view angle and the distances covered, surface elevations not accounted for had a significant effect on the coordinate transformation. A 20 cm wave displaces the horizontal position by 18 m (along the line of sight) at the channel entrance. To alleviate this effect, the world coordinate pairs were transformed to polar (r;q) with respect to the camera, and a low-pass filter with 50 s cutoff period was applied to the radial coordinate time-series to filter-out the high- frequency wind- and boat-generated waves. The filtered polar coordinates of each tracer were converted to UTM (Universal Transverse Mercator), then the velocities in the east (E) and north (N) directions were computed using the central difference scheme. Inaccuracies in the coordinate transformation and tracer centre identification necessitated the use of another low-pass filter in the velocity time-series, in each direction. 6.1.2 GPS. The spherical coordinate pairs extracted form the DEX GPS unit were first converted to UTM projected coordinates, and the velocities in the east and north directions were computed also using a central difference scheme. We note that although the sampling frequency was set to 2 Hz, at low flow speeds the sampling frequency was automatically reduced, since the displacement within the sampling period was not sufficient to register a new point, given the accuracy of the receiver. Sudden jumps in velocity appearing when the sampling frequency changed were removed. GPS receivers inherently suffer from positioning inaccuracies due to satellite- and receiver- dependent errors, ionospheric and tropospheric delays and multi-path errors (Uren and Price, 1994). The absolute position error for non differential receivers can reach up to 3.6m in each 151 horizontal direction (Johnson and Pattiaratchi, 2004). However, errors in velocity computation are rather affected by the relative position accuracy, namely the temporal changes of position accuracy. The total velocity error, which is also a function of the finite differencing error, limits the motions that can be resolved. It has been shown (Johnson and Pattiaratchi, 2004) that the signal-to-noise ratio increases for lower frequency motions and becomes acceptable to resolve frequencies below 0.05 Hz, for 1 Hz sampling. By analogy, we low-pass filtered the velocity time-series using a 0.025 Hz cutoff. The DEX floater was also traced optically, when it appeared in the camera FOV . The velocities calculated from the GPS data at those times were used to validate the optical particle tracking methodology (Figure 6.5). 6.1.3 Numerical Modelling. Our initial conditions are derived from NOAA’s (National Oceanographic and Atmospheric Administration) tsunami-inverted source (Utku Kanoglu, personal communication), obtained using standard procedures (Percival et al., 2011). The tsunami was propagated across the Pacific Ocean using a 4 arc-min grid, and through four telescopic grids of increasing resolution in the near field, while a 1/3 arc-sec resolution (Caldwell et al., 2011) was used in Ventura Harbour and, referenced to the mean low water datum. The overall accuracy of the simulations was evaluated by comparing the calculated surface elevation time series at Ventura’s tide gauge station with the de-tided record- ing (Figure 6.3e). The arrival time, phase and amplitude of the first three waves compare well with the recording. We also used MOST to estimate the background tidal flow by forcing the model boundaries with the predicted tidal amplitude. The simulated maximum currents in the harbour entrance chan- nel were below 20 cm/s over three tidal cycles. During the field measurements (flood tide), they were less than 10 cm/s. 152 Figure 6.4: Projected tracer paths in the harbor entrance channel. Continuous (-) and dash-dot (-.) curves correspond to paths extracted from optical tracking, referenced to coordinate systems CS 1 and CS 2 respectively. Dashed (- -) curves correspond to GPS tracks, all referenced to CS 1. Figure 6.5: Scatter plot of velocities measured with the DEX GPS unit versus velocities measured optically. Squares correspond to velocity magnitude, crosses to the velocity in the eastern direction, and stars to the velocity in the north direction. The R 2 coefficients of determination are noted. 153 Chapter 7 Conclusions The objective of this thesis is to study the generation and evolution of tsunami-induced coherent structures in a well-controlled environment using realistic scaling. A physical configuration in the image of a port entrance was created in a large-scale wave basin. The flow from the generated small-amplitude, long period wave was accelerated through the scaled harbor channel. The long wave momentum and energy was transferred to the starting-jet vortex at the channel. The vortex detached from the trailing jet resulting in a self-propagating coherent structure. The experimental layout was optimized to allow a large coherent structure to form in the offshore basin and data collection extended long-enough to capture the growth and evolution of the TCS until its radial growth became limited by the vertical boundaries of the offshore basin. Surface velocities were successfully extracted using both 2D and 3D PTV techniques. The 2D PTV methodology was applied to the broader study area of the basin, whereas 3D PTV was only applied around the harbor channel. The errors from the inherent 2D PTV simplifications were dis- cussed. The application of the methodology was validated by comparing the extracted data to the measured surface elevation and ADV velocity timeseries. The comparison shows good agreement. V orticity measurements at the TCS core were also collected using two tracer-configurations. The generation and growth of the starting-jet vortex were studied through instantaneous mean velocity fields evaluated on a stationary-ensemble. V orticity maps created from the velocity fields were used to calculate the total circulation build-up in the basin. The vorticity transport equation 154 was simplified to a obtain an equation for the vorticity influx rate through the channel. The compar- ison between the total circulation growth and the vorticity influx rate shows good agreement until the turbulence intensity introduces high estimation uncertainties to the channel vorticity profiles. The vortex boundary was defined through thel ci = 0 contour to infer the vortex circulation growth rate. The dimensionless pinch-off time of the inshore TCS compares well with the observations of Nicolau del Roure et al. (2009), whereas the pinch-off of the offshore TCS occurred somewhat earlier than predicted. The growth of the offshore starting-jet vortex provided an estimate of the flow rate supplied to the vortex core. Finally, the vortex advection velocity (resulting from the vertical boundaries reactionary forces to the vortex flow) was found to be in agreement with the observations of Afanasyev (2006) for vortex dipoles. The later development of the TCS was studied through a TCS-centered ensemble. The scattered velocity vectors from the individual trials were all referenced to the TCS center. The vorticity maps revealed a compact core of intense positive vorticity surrounded by a ring of opposite-signed vorticity which forms from the interaction of the vortex boundary with the surrounding irrotational fluid. The flow divergence maps show that the TCS center exhibits downwelling (coupled with a strong depression) and the surrounding flow experiences upwelling in a re-circulating 3D pattern. The structure of the monopolar vortex was assessed through its velocity components (u r ;u q ) in cylindrical coordinates, so as to characterize its profile through known types of geophysical vor- tices. In the first stages of TCS development, while the radial velocity creates zones of surface flow convergence/divergence, the u q profiles better match a vortex profile valid on the fplane. The u q profile starts to converge towards the Gaussian vortex a-profile at t 390s, which corresponds to the time of rapid decay of the radial velocity component. The vortex bulk radius grows from the initial stages in a rate compatible with viscous diffusion, while the vortex core radius remains stationary until the radial velocity decay. Once the vortex core starts growing, the vortex profile relaxes towards the a-profile, with a 0:34. The decay of the maximum azimuthal velocity due to bottom friction matches the prediction of a first-order force-balance model. The bottom shear stress corresponding to a turbulent boundary layer was found to explain the kinetic decay until the end of data collection and thus laminar flow 155 conditions were not reached. The transition time-scale from turbulent to laminar flow, can be estimated from the vertical viscous relaxation time to be in the order of t 3:4 hrs, which extends beyond the duration of the trials. The kinetic energy of the main and secondary flow components in a fixed radius from the TCS center were evaluated to characterize the two-dimensionality of the shallow flow conditions. It was demonstrated that the radial kinetic energy experiences a rapid decay around t 390s and it can be considered insignificant compared to the azimuthal kinetic energy thereafter. Thus, it can be considered as the transition point of the flow to quasi-two-dimensional. Preliminary analysis identified basin sloshing as the source of the fluctuations in the radial kinetic energy, but at this point it remains inconclusive. The last chapter of this thesis presented observations and surface flow measurements in Ven- tura harbor during the Chilean 2015 tsunami. Rectified camera frames revealed the formation of starting-jets that grow in size and form coherent structure in the harbor entrance channel. Flow measurements through surface tracers showed that the surface current speeds reached 1.5 m/s in a highly turbulent but sub-critical flow. The measurements compared well with the NSWE numerical model’s mean channel-flow envelope prediction. Suggestions for Further Research The TCS circulation growth analysis presented in Chapter 4 requires further work to be com- plete. The starting-jet boundary, defined from the swirl-strength l ci = 0 contour possibly results in very different timescales of detachment from the trailing-jet (pinch-off) compared to defining the boundary using vorticity iso-contours. The computed TCS circulation before pinch-off is also expected to be different depending on the TCS boundary definition. Thus, more work is needed to examine the implications of the above in the circulation growth analysis and provide alternative pinch-off definitions. Chapter 5 provided data on the TCS structure evolution and revealed a strong radial velocity component until t 390s. The effect of the radial velocity component on the TCS core radius can ideally be illustrated analytically (to first-order). This analysis framework should also reveal the 156 physical mechanism behind the deformation of the the azimuthal velocity profile until the radial velocity becomes insignificant. The observations regarding the magnitude of the secondary flow can be extended to earlier times (during the TCS generation) through the stationary-ensemble. This would allow to examine the build-up timescales of the radial velocity kinetic energy E r . Furthermore, the pattern of E r decay should be further studied to determine whether it’s physical (a re-circulation in the TCS structure) and provide an overall estimate of the evaluation uncertainty. The results regarding the transition of the flow to quasi-two-dimensional are not expected to be affected since individual observations collectively match and point at the same time-scale. 157 Reference List Abramowitz, M. and Stegun, I. A. Handbook of mathematical functions with formulas, graphs, and mathematical tables. U.S. Department of Commerce, National Bureau of Standards; Dover publications, applied mathematic series, tenth printing, with corrections edition, 1972. ISBN 0-486-61272-4. Admire, A. R., Dengler, L. A., Crawford, G. B., Uslu, B. U., Borrero, J. C., Greer, D., and Wilson, R. I. Observed and modeled currents from the tohoku-oki, japan and other recent tsunamis in northern california. Pure and Applied Geophysics, 171(12):3385–3403, 2014. doi: 10.1007/ s00024-014-0797-8. Adrian, R. J., Christensen, K. T., and Liu, Z.-C. Analysis and interpretation of instantaneous turbu- lent velocity fields. Experiments in Fluids, 29(3):275–290, 2000. doi: 10.1007/s003489900087. Afanasyev, Y . D. Formation of vortex dipoles. Physics of Fluids, 18(3):037103, 2006. doi: 10.1063/1.2182006. Ag¨ u´ ı, J. C. and Jim´ enez, J. On the performance of particle tracking. Journal of Fluid Mechanics, 185:447468, 1987. doi: 10.1017/S0022112087003252. Akkermans, R. A. D., Cieslik, A. R., Kamp, L. P. J., Trieling, R. R., Clercx, H. J. H., and van Heijst, G. J. F. The three-dimensional structure of an electromagnetically generated dipolar vortex in a shallow fluid layer. Physics of Fluids, 20(11):116601, 2008. doi: 10.1063/1.3005452. Ayca, A. and Lynett, P. J. Effect of tides and source location on nearshore tsunami-induced cur- rents. Journal of Geophysical Research: Oceans, 2016. doi: 10.1002/2016JC012435. Ayca, A., Lynett, P. J., Borrero, J. C., Miller, K., and Willson, R. Numerical and physical model- ing of localized tsunami-induced currents in harbors. Coastal Engineering Proceedings, 1(34): currents.6, 2012. doi: 10.9753/icce.v34.currents.6. Aydemir, E., Worth, N. A., and Dawson, J. R. The formation of vortex rings in a strongly forced round jet. Experiments in Fluids, 52(3):729–742, 2012. doi: 10.1007/s00348-011-1110-6. Barberopoulou, A., Borrero, J. C., Uslu, B., Legg, M. R., and Synolakis, C. E. A Second Genera- tion of Tsunami Inundation Maps for the State of California. Pure and Applied Geophysics, 168 (11):2133–2146, 2011a. doi: 10.1007/s00024-011-0293-3. Barberopoulou, A., Legg, M. R., Uslu, B., and Synolakis, C. E. Reassessing the tsunami risk in major ports and harbors of california i: San diego. Natural Hazards, 58(1):479–496, 2011b. doi: 10.1007/s11069-010-9681-8. Batchelor, G. K. An introduction to fluid dynamics. Cambridge University Press, Cambridge, UK, 1967. Beckers, M., Verzicco, R., Clerx, H., and van Heijst, G. Dynamics of pancake-like vortices in a stratified fluid: experiments, model and numerical simulations. Journal of Fluid Mechanics, 433:127, 2001. doi: 10.1017/S0022112001003482. 158 Boffetta, G., Celani, A., and Vergassola, M. Inverse energy cascade in two-dimensional turbulence: Deviations from gaussian behavior. Physical Review E, 61:R29, 2000. Borrero, J. C., Bell, R., Csato, C., DeLange, W., Goring, D., Dougal Greer, S., Pickett, V ., and Power, W. Observations, effects and real time assessment of the march 11, 2011 tohoku-oki tsunami in new zealand. Pure and Applied Geophysics, 170(6-8):1229–1248, 2013. doi: 10. 1007/s00024-012-0492-6. Borrero, J. C., Lynett, P. J., and Kalligeris, N. Tsunami currents in ports. Philosophical Transac- tions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 373 (2053), 2015. doi: 10.1098/rsta.2014.0372. Bouguet, J.-Y . Camera Calibration Toolbox for Matlab, 2014. URL http://www.vision. caltech.edu/bouguetj/calib_doc/. Bourgeois, F. and Lassalle, J.-C. An extension of the munkres algorithm for the assignment prob- lem to rectangular matrices. Commun. ACM, 14(12):802–804, 1971. doi: 10.1145/362919. 362945. Bruins, H. J., MacGillivray, J. A., Synolakis, C. E., Benjamini, C., Keller, J., Kisch, H. J., Kl¨ ugel, A., and van der Plicht, J. Geoarchaeological tsunami deposits at palaikastro (crete) and the late minoanfIAg eruption of santorini. Journal of Archaeological Science, 35(1):191 – 212, 2008. doi: http://dx.doi.org/10.1016/j.jas.2007.08.017. Bryant, D. B., Seol, D.-G., and Socolofsky, S. A. Quantification of turbulence properties in bubble plumes using vortex identification methods. Physics of Fluids, 21(7):075101, 2009. doi: http: //dx.doi.org/10.1063/1.3176464. Bryant, D. B., Whilden, K. A., Socolofsky, S. A., and Chang, K.-A. Formation of tidal starting-jet vortices through idealized barotropic inlets with finite length. Environmental Fluid Mechanics, 12(4):301–319, 2012. doi: 10.1007/s10652-012-9237-4. Caldwell, R., Taylor, L., Eakins, B., Carignan, K., Grothe, P., Lim, E., and Friday, D. Digital eleva- tion models of santa monica, california: Procedures, data sources and analysis. NOAA Technical Memorandum NESDIS NGDC-46, page 38 pp, 2011. doi: 10.1007/s10652-012-9237-4. Capart, H., Young, D. L., and Zech, Y . V oronoi imaging methods for the measurement of granular flows. Experiments in Fluids, 32:121–135, 2002. doi: 10.1007/s003480100351. Carmer, C. F. Shallow turbulent wake flows: momentum and mass transfer due to large-scale coherent vortical structures. PhD thesis, Universit¨ at Karlsruhe, 2005. Carnevale, G. F. and Kloosterziel, R. C. Emergence and evolution of triangular vortices. Journal of Fluid Mechanics, 259:305331, 1994. doi: 10.1017/S0022112094000157. Carton, X. J., Flierl, G. R., and Polvani, L. M. The generation of tripoles from unstable axisym- metric isolated vortex structures. EPL (Europhysics Letters), 9(4):339, 1989. 159 Casa, L. D. C. and Krueger, P. S. Radial basis function interpolation of unstructured, three- dimensional, volumetric particle tracking velocimetry data. Measurement Science and Tech- nology, 24(6), 2013. doi: 101088/0957-0233/24/6/065304. Choi, B. H., Pelinovsky, E., Kim, D. C., Kim, K. O., and Kim, K. H. Three-dimensional simulation of the 1983 central east (japan) sea earthquake tsunami at the imwon port (korea). Ocean Engi- neering, 35(1415):1545 – 1559, 2008. doi: http://dx.doi.org/10.1016/j.oceaneng.2008.06.002. Crocker, J. C. and Grier, D. G. Methods of digital video microscopy for colloidal studies. Journal of colloid and interface science, 179(1):298–310, 1996. doi: 10.1006/jcis.1996.0217. Dabiri, J. O. Optimal vortex formation as a unifying principle in biological propulsion. Annual Review of Fluid Mechanics, 41(1):17–33, 2009. doi: 10.1146/annurev.fluid.010908.165232. Didden, N. On the formation of vortex rings: Rolling-up and production of circulation. Zeitschrift f¨ ur angewandte Mathematik und Physik ZAMP, 30(1):101–116, 1979. doi: 10.1007/ BF01597484. Discetti, S., Ag¨ uera, N., Cafiero, G., and Astarita, T. Ensemble 3d ptv for high resolution turbulent statistics. 11th international symposium on particle image velocimetry-PIV15, Santa Barbara, CA, September 14-16, 2015. Dolzhanskii, F. V ., Krymov, V . A., and Manin, D. Y . An advanced experimental investigation of quasi-two-dimensional shear flows. Journal of Fluid Mechanics, 241:705–722, 1992. doi: 10.1017/S0022112092002209. Douxchamps, D., Devriendt, D., Capart, H., Craeye, C., Macq, B., and Zech, Y . Stereoscopic and velocimetric reconstructions of the free surface topography of antidune flows. Experiments in Fluids, 39(3):535–553, 2005. doi: 10.1007/s00348-005-0983-7. Douxchamps, D., Spinewine, B., Capart, H., Zech, Y ., and Macq, B. Particle-based imaging methods for the characterisation of complex fluid flows. In Proc of the IEEE Oceans Conf, 2002. Dracos, T., Giger, M., and Jirka, G. H. Plane turbulent jets in a bounded fluid layer. Journal of Fluid Mechanics, 241:587614, Aug 1992. doi: 10.1017/S0022112092002167. Driessen, J. and MacDonald, C. F. The troubled island. minoan crete before and after the san- torini eruption. Technical report, In: Aegaeum 17. Universit´ e de Li´ ege/University of Texas, Li´ ege/Austin., 1997. Duran-Matute, M., Kamp, L. P. J., Trieling, R. R., and van Heijst, G. J. F. Scaling of decaying shallow axisymmetric swirl flows. Journal of Fluid Mechanics, 648:471484, 2010. doi: 10. 1017/S0022112009994034. Einstein, A. ¨ Uber die von der molekularkinetischen Theorie der W¨ arme geforderte Bewegung von in ruhenden Fl¨ ussigkeiten suspendierten Teilchen. Annalen der Physik, 322(8):549–560, 1905. Faugeras, O. and Luong, Q.-T. The geometry of multiple images. MIT Press, 2001. 160 Fl´ or, J. B. and Van Heijst, G. J. F. An experimental study of dipolar vortex structures in a stratified fluid. Journal of Fluid Mechanics, 279:101133, 1994. doi: 10.1017/S0022112094003836. Flor, J. B. and Van Heijst, G. J. F. Stable and unstable monopolar vortices in a stratified fluid. Journal of Fluid Mechanics, 311:257287, 1996. doi: 10.1017/S0022112096002595. Fl´ or, J.-M. and Eames, I. Dynamics of monopolar vortices on a topographic beta-plane. Journal of Fluid Mechanics, 456:353–376, 2002. doi: 10.1017/S0022112001007728. Fritz, H. M., Phillips, D. A., Okayasu, A., Shimozono, T., Liu, H., Mohammed, F., Skanavis, V ., Synolakis, C. E., and Takahashi, T. The 2011 Japan tsunami current velocity measurements from survivor videos at Kesennuma Bay using LiDAR: JAPAN TSUNAMI VIDEO LIDAR. Geophysical Research Letters, 39(7):n/a–n/a, 2012. doi: 10.1029/2011GL050686. Fritz, H., Petroff, C., Cataln, P., Cienfuegos, R., Winckler, P., Kalligeris, N., Weiss, R., Barrientos, S., Meneses, G., Valderas-Bermejo, C., Ebeling, C., Papadopoulos, A., Contreras, M., Almar, R., Dominguez, J., and Synolakis, C. Field survey of the 27 february 2010 chile tsunami. Pure and Applied Geophysics, 168(11):1989–2010, 2011. doi: 10.1007/s00024-011-0283-5. Gharib, M., Rambod, E., and Shariff, K. A universal time scale for vortex ring formation. Journal of Fluid Mechanics, 360:121140, 1998. doi: 10.1017/S0022112097008410. Hayashi, S. and Koshimura, S. The 2011 tohoku tsunami flow velocity estimation by the aerial video analysis and numerical modeling. Journal of Disaster Research, 8(4):561572, 2013. Heikkila, J. Geometric camera calibration using circular control points. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(10):1066–1077, 2000. doi: 10.1109/34.879788. Holland, T. K., Holman, R., Lippmann, T. C., Stanley, J., Plant, N., and others. Practical use of video imagery in nearshore oceanographic field studies. IEEE Journal of Oceanic Engineering, 22(1):81–92, 1997. Hussain, A. K. M. F. Coherent structures-reality and myth. Physics of Fluids, 26(10):2816–2850, 1983. doi: 10.1063/1.864048. Jeong, J. and Hussain, F. On the identification of a vortex. Journal of Fluid Mechanics, 285:69–94, 1995. Jirka, G. H. and Uijttewaal, W. S. J. Shallow flows: a definition. Proceedings of the international conference on shallow flows, pages 3–11, 2004. Jirka, G. H. Large scale flow structures and mixing processes in shallow flows. Journal of Hydraulic Research, 39(6):567–573, 2001. doi: 10.1080/00221686.2001.9628285. Johnson, D. and Pattiaratchi, C. Application, modelling and validation of surfzone drifters. Coastal Engineering, 51(5-6):455–471, 2004. doi: 10.1016/j.coastaleng.2004.05.005. J¨ uttner, B., Marteau, D., Tabeling, P., and Thess, A. Numerical simulations of experiments on quasi-two-dimensional turbulence. Physical Review E, 55:54795488, 1997. doi: 10.1103/ PhysRevE.55.5479. 161 Kalligeris, N., Skanavis, V ., Tavakkol, S., Ayca, A., Safty, H. E., Lynett, P., and Synolakis, C. Lagrangian flow measurements and observations of the 2015 chilean tsunami in ventura, ca. Geophysical Research Letters, 43(10):5217–5224, 2016. doi: 10.1002/2016GL068796. 2016GL068796. Kashiwai, M. Tidal residual circulation produced by a tidal vortex. Journal of the Oceanographical Society of Japan, 40(4):279–294, 1984. doi: 10.1007/BF02302521. Kilfoil, M. and Pelletier, V . Mechanics of the Cell (Software), http://people.umass.edu/kilfoil/downloads.html, Last accessed on 2015-08-28., 2015. URL http://people.umass.edu/kilfoil/downloads.html. Kim, D.-H. and Lynett, P. J. Turbulent mixing and passive scalar transport in shallow flows. Physics of Fluids, 23(1):016603, 2011. doi: 10.1063/1.3531716. Kim, K. O., Pelinovsky, E., and Jung, K. T. Three-dimensional simulation of 2011 east japan- off pacific coast earthquake tsunami induced vortex flows in the oarai port. Journal of Coastal Research, 1:284–289, 2013. doi: http://dx.doi.org/10.2112/SI65-049.1. Kloosterziel, R. C. On the large-time asymptotics of the diffusion equation on infinite domains. Journal of Engineering Mathematics, 24(3):213–236, 1990. doi: 10.1007/BF00058467. Kolar, V . V ortex identification: New requirements and limitations. International Journal of Heat and Fluid Flow, 28(4):638–652, 2007. doi: http://dx.doi.org/10.1016/j.ijheatfluidflow.2007.03. 004. Kraichnan, R. H. Inertial ranges in two-dimensional turbulence. Physics of Fluids, 10(7):1417– 1423, 1967. doi: http://dx.doi.org/10.1063/1.1762301. Krieg, M. and Mohseni, K. Modelling circulation, impulse and kinetic energy of starting jets with non-zero radial velocity. Journal of Fluid Mechanics, 719:488–526, Mar 2013. doi: 10.1017/ jfm.2013.9. Krueger, P. S. An over-pressure correction to the slug model for vortex ring circulation. Journal of Fluid Mechanics, 545:427–443, 2005. doi: 10.1017/S0022112005006853. Lacy, J. R., Rubin, D. M., and Buscombe, D. Currents, drag, and sediment transport induced by a tsunami. Journal of Geophysical Research: Oceans, 117(C9):n/a–n/a, 2012. doi: 10.1029/ 2012JC007954. C09028. Lee, D. T. and Schachter, B. J. Two algorithms for constructing a delaunay triangulation. Inter- national Journal of Computer & Information Sciences, 9(3):219–242, 1980. doi: 10.1007/ BF00977785. Lloyd, P. M., Stansby, P. K., and Ball, D. J. Unsteady surface-velocity field measurement using particle tracking velocimetry. Journal of Hydraulic Research, 33(4):519–534, 1995. doi: 10. 1080/00221689509498658. 162 Lynett, P., Gately, K., Wilson, R., Montoya, L., and others, . Inter-model analysis of tsunami- induced coastal currents. in preparation, 2017. Lynett, P. J., Borrero, J. C., Weiss, R., Son, S., Greer, D., and Renteria, W. Observations and modeling of tsunami-induced currents in ports and harbors. Earth and Planetary Science Letters, 327-328:68–74, 2012. doi: 10.1016/j.epsl.2012.02.002. Lynett, P. J., Borrero, J., Son, S., Wilson, R., and Miller, K. Assessment of the tsunami- induced current hazard. Geophysical Research Letters, 41(6):2048–2055, 2014. doi: 10.1002/ 2013GL058680. 2013GL058680. Malik, N. A., Dracos, T., and Papantoniou, D. A. Particle tracking velocimetry in three- dimensional flows. Experiments in Fluids, 15(4-5):279–294, 1993. Maravelakis, N., Kalligeris, N., Lynett, P., Skanavis, V ., and Synolakis, C. A study of wave ampli- fication in the venetian harbor in chania. Coastal Engineering Proceedings, 1(34):waves.59, 2014. doi: 10.9753/icce.v34.waves.59. Mcwilliams, J. C. The vortices of two-dimensional turbulence. Journal of Fluid Mechanics, 219: 361385, 1990. doi: 10.1017/S0022112090002981. Montoya, L., Lynett, P., Thio, H.-T., and Li, W. Spatial statistics of tsunami overland flow proper- ties. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2016. doi: 10.1061/(ASCE) WW.1943-5460.0000363. Mori, N. Mace softwares. http://www.oceanwave.jp/softwares/mace/index. php?MACE%20Softwares, 2014. [Online; accessed 30-Dec.-2015]. Mori, N. and Tomoyuki, T. Nationwide post event survey and analysis of the 2011 tohoku earthquake tsunami. Coastal Engineering Journal, 54(01):1250001, 2012. doi: 10.1142/ S0578563412500015. Mori, N., Suzuki, T., and Kakuno, S. Noise of acoustic doppler velocimeter data in bubbly flows. Journal of Engineering Mechanics, 133(1):122–125, 2007. doi: 10.1061/(ASCE) 0733-9399(2007)133:1(122). Nicolau del Roure, F., Socolofsky, S. A., and Chang, K.-A. Structure and evolution of tidal starting jet vortices at idealized barotropic inlets. Journal of Geophysical Research: Oceans, 114:n/a– n/a, 2009. doi: 10.1029/2008JC004997. C05024. Nortek. Datasheet vectrino lab. www.nortek-as.com/lib/data-sheets/ datasheet-vectrino-lab/view, 2015. [Online; accessed 30-Dec.-2015]. Okal, E. A. The quest for wisdom: lessons from 17 tsunamis, 2004–2014. Philosophical Transac- tions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 373 (2053), 2015. doi: 10.1098/rsta.2014.0370. Okal, E. A. and Synolakis, C. E. Sequencing of tsunami waves: why the first wave is not always the largest. Geophysical Journal International, 204(2):719–735, 2015. doi: 10.1093/gji/ggv457. 163 Okal, E. A., Fritz, H. M., Raad, P. E., Synolakis, C., Al-Shijbi, Y ., and Al-Saifi, M. Oman field survey after the december 2004 indian ocean tsunami. Earthquake Spectra, 22(S3):203–218, 2006a. doi: 10.1193/1.2202647. Okal, E. A., Fritz, H. M., Raveloson, R., Joelson, G., Panokov, P., and Rambolamanana, G. Mada- gascar field survey after the december 2004 indian ocean tsunami. Earthquake Spectra, 22(S3): 263–283, 2006b. doi: 10.1193/1.2202646. Okal, E. A., Sladen, A., and Okal, E. A.-S. Rodrigues, mauritius, and reunion islands field survey after the december 2004 indian ocean tsunami. Earthquake Spectra, 22(S3):241–261, 2006c. doi: 10.1193/1.2209190. Panasonic. Panasonic AW-HE60s Full-HD Integrated Pan-Tilt Camera. http://business. panasonic.co.uk/professional-camera/integrated-pt-cameras/ aw-he60s, 2015. [Online; accessed 13-July-2015]. Paret, J. and Tabeling, P. Experimental observation of the two-dimensional inverse energy cascade. Phys. Rev. Lett., 79:4162–4165, 1997. doi: 10.1103/PhysRevLett.79.4162. Paul, M. P. Development of a vision-based particle tracking velocimetry method and post-processing of scattered velocity data. Master’s thesis, University of Washing- ton, 2012. URL https://digital.lib.washington.edu/researchworks/ handle/1773/20280?show=full. Percival, D., Denbo, D., Ebl, M., Gica, E., Mofjeld, H., Spillane, M., Tang, L., and Titov, V . Extraction of tsunami source coefficients via inversion of dart buoy data. Natural Hazards, 58 (1):567–590, 2011. doi: 10.1007/s11069-010-9688-1. Qian, H., Sheetz, M. P., and Elson, E. L. Single particle tracking. Analysis of diffusion and flow in two-dimensional systems. Biophysical journal, 60(4):910, 1991. Raffel, M., Willert, C., Wereley, S., and Kompenhans, J. Particle Image Velocimetry: A Practical Guide. Springer-Verlag Berlin Heidelberg, Berlin; Heidelberg; New York, 2nd ed edition, 2007. ISBN 978-3-540-72307-3. Sandwell, D. T. Biharmonic spline interpolation of geos-3 and seasat altimeter data. Geophysical Research Letters, 14(2):139–142, 1987. doi: 10.1029/GL014i002p00139. Satijn, M. P., Cense, A. W., Verzicco, R., Clercx, H. J. H., and van Heijst, G. J. F. Three- dimensional structure and decay properties of vortices in shallow fluid layers. Physics of Fluids, 13(7):1932–1945, 2001. doi: http://dx.doi.org/10.1063/1.1374936. Seol, D.-G. and Jirka, G. H. Quasi-two-dimensional properties of a single shallow-water vortex with high initial Reynolds numbers. Journal of Fluid Mechanics, 665:274–299, 2010. doi: 10.1017/S0022112010003915. Shariff, K. and Leonard, A. V ortex rings. Annual Review of Fluid Mechanics, 24(1):235–279, 1992. doi: 10.1146/annurev.fl.24.010192.001315. 164 Socolofsky, S. A. and Jirka, G. H. Large-scale flow structures and stability in shallow flows. Journal of Environmental Engineering and Science, 3(5):451–462, 2004. doi: 10.1139/S04-032. Son, S., Lynett, P. J., and Kim, D.-H. Nested and multi-physics modeling of tsunami evolution from generation to inundation. Ocean Modelling, 38(12):96 – 113, 2011. doi: http://dx.doi.org/ 10.1016/j.ocemod.2011.02.007. Sous, D., Bonneton, N., and Sommeria, J. Turbulent vortex dipoles in a shallow water layer. Physics of Fluids, 16(8):2886–2898, 2004. doi: 10.1063/1.1762912. Sous, D., Bonneton, N., and Sommeria, J. Transition from deep to shallow water layer: formation of vortex dipoles. European Journal of Mechanics - B/Fluids, 24(1):19 – 32, 2005. doi: http: //dx.doi.org/10.1016/j.euromechflu.2004.06.002. Spedding, G. R. and Rignot, E. J. M. Performance analysis and application of grid inter- polation techniques for fluid flows. Experiments in Fluids, 15(6):417–430, 1993. doi: 10.1007/BF00191784. St¨ uer, H. and Blaser, S. Interpolation of scattered 3d ptv data to a regular grid. Flow, Turbulence and Combustion, 64(3):215–232, 2000. doi: 10.1023/A:1009904013148. Suppasri, A., Mas, E., Koshimura, S., Imai, K., Harada, K., and Imamura, F. Develop- ing tsunami fragility curves from the surveyd data of the 2011 great east japan tsunami in sendai and ishinomaki plains. Coastal Engineering Journal, 54(01):1250008, 2012. doi: 10.1142/S0578563412500088. Synolakis, C. E., Bernard, E. N., Titov, V . V ., Kˆ ano˘ glu, U., and Gonz´ alez, F. I. Validation and Verification of Tsunami Numerical Models. Pure and Applied Geophysics, 165(11-12):2197– 2228, 2008. doi: 10.1007/s00024-004-0427-y. Tavakkol, S. and Lynett, P. Celeris: A GPU-accelerated open source software with a Boussinesq- type wave solver for real-time, interactive simulation and visualization. ArXiv e-prints, 2016. Titov, V . V . and Synolakis, C. E. Extreme inundation flows during the hokkaido-nansei-oki tsunami. Geophysical Research Letters, 24(11):1315–1318, 1997. doi: 10.1029/97GL01128. Titov, V ., Gonzalez, F., Bernard, E., Eble, M., Mofjeld, H., Newman, J., and Venturato, A. Real- time tsunami forecasting: Challenges and solutions. Natural Hazards, 35(1):35–41, 2005. doi: 10.1007/s11069-004-2403-3. Tsai, V . C., Ampuero, J.-P., Kanamori, H., and Stevenson, D. J. Estimating the effect of earth elasticity and variable water density on tsunami speeds. Geophysical Research Letters, 40(3): 492–496, 2013. doi: 10.1002/grl.50147. Uijttewaal, W. S. J. and Booij, R. Effects of shallowness on the development of free-surface mixing layers. Physics of Fluids, 12(2):392–402, 2000. doi: http://dx.doi.org/10.1063/1.870317. Uren, J. and Price, W. Surveying for engineers (3rd ed). Macmillan, Basingstoke, 1994. 165 van Heijst, G. and Clercx, H. Laboratory modeling of geophysical vortices. Annual Review of Fluid Mechanics, 41(1):143–164, 2009. doi: 10.1146/annurev.fluid.010908.165207. Weitbrecht, V ., Khn, G., and Jirka, G. Large scale piv-measurements at the surface of shallow water flows. Flow Measurement and Instrumentation, 13:237 – 245, 2002. doi: http://dx.doi. org/10.1016/S0955-5986(02)00059-6. Wells, M. G. and van Heijst, G.-J. F. A model of tidal flushing of an estuary by dipole formation. Dynamics of Atmospheres and Oceans, 37:223–244, 2003. doi: 10.1016/j.dynatmoce.2003.08. 002. Whilden, K. A., Socolofsky, S. A., Chang, K.-A., and Irish, J. L. Using surface drifter observations to measure tidal vortices and relative diffusion at aransas pass, texas. Environmental Fluid Mechanics, 14(5):1147–1172, 2014. doi: 10.1007/s10652-014-9361-4. Wilson, R. I., Admire, A. R., Borrero, J. C., Dengler, L. A., Legg, M. R., Lynett, P., McCrink, T. P., Miller, K. M., Ritchie, A., Sterling, K., and Whitmore, P. M. Observations and impacts from the 2010 chilean and 2011 japanese tsunamis in california (usa). Pure and Applied Geophysics, 170 (6):1127–1147, 2013. doi: 10.1007/s00024-012-0527-z. Wilson, R. I., Lynett, P., Miller, K. M., Admire, A. R., Ayca, A., Curtis, E., Dengler, L. A., Hornick, M., Nicolini, T., and Peterson, D. Maritime tsunami response playbooks: Background information and guidance for response and hazard mitigation use. Technical report, California Geological Survey Special Report 241, 2016. 166

Abstract (if available)

Abstract The aim of this thesis is to study the generation and evolution of tsunami-induced coherent structures in a well-controlled environment using realistic scaling. The results of this study are intended to provide a step forward towards unraveling the complexities involved in how tsunami-induced coherent structures affect our coastal infrastructure. ❧ A physical configuration was created in the tsunami wave basin of Oregon State University, in the image of a port entrance. A small-amplitude, long period wave is generated, and the flow is accelerated through the harbor channel. A separated region forms, which when coupled with the transient flow, leads to the formation of a turbulent coherent structure (TCS). The wave forcing translates to realistic tsunami prototype scales. The fully turbulent and subcritical flow conditions suggest that the laboratory results are scale independent and applicable to geophysical flows. ❧ Surface velocities are inferred from surface tracers using 2D and 3D PTV techniques. Surface elevation data was collected throughout the basin to provide an overview of the long wave motion. Acoustic Doppler Velocimeter data were also collected to provide external forcing information and data on coupling between the external tsunami forcing and the generated TCS. ❧ The experimental layout was optimized to allow a large coherent structure to form in the offshore basin. The experimental data collection extended long-enough to capture the growth and evolution of the TCS until the TCS radial growth became limited by the vertical boundaries. The circulation growth of the TCS is compared with studies on tidal flushing in near-shore inlets. The structure of the resulting monopolar coherent structure is assessed through its velocity components in cylindrical coordinates, so as to characterize its profile through other known types of geophysical vortices. Kinematic energy decay due to bottom friction and the growth of the vortex structure are compared with analytical models. The three-dimensional motions are examined through the secondary flow components to characterize the two-dimensionality of the shallow flow.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Wave induced hydrodynamic complexity and transport in the nearshore

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Open channel flow instabilities: modeling the spatial evolution of roll waves

PDF

Re-assessing local structures of turbulent flames via vortex-flame interaction

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Understanding the transport potential of nearshore tsunami currents

PDF

Pattern generation in stratified wakes

PDF

The development of a hydraulic-control wave-maker (HCW) for the study of non-linear surface waves

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Nonlinear long wave amplification in the shadow zone of offshore islands

PDF

Investigations of fuel effects on turbulent premixed jet flames

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

High-resolution imaging and monitoring of fault zones and shallow structures: case studies in southern California and on Mars

PDF

Boundary layer and separation control on wings at low Reynolds numbers

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Flow and thermal transport at porous interfaces

PDF

Contributions to structural and functional retinal imaging via Fourier domain optical coherence tomography

Asset Metadata

Creator Kalligeris, Nikos (author)

Core Title Tsunami-induced turbulent coherent structures

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 02/15/2017

Defense Date 01/18/2017

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

Tag OAI-PMH Harvest,ports and harbors,shallow water flows,tsunami currents,tsunami hazard,tsunami-induced eddy,turbulent coherent structures

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Lynett, Patrick (committee chair), Spedding, Geoff (committee member), Synolakis, Costas (committee member)

Creator Email kalliger@usc.edu,nkalligeris@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-337745

Unique identifier UC11258927

Identifier etd-Kalligeris-5061.pdf (filename),usctheses-c40-337745 (legacy record id)

Legacy Identifier etd-Kalligeris-5061.pdf

Dmrecord 337745

Document Type Dissertation

Rights Kalligeris, Nikos

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