University of Southern California Dissertations and Theses

Understanding the transport potential of nearshore tsunami currents

(USC Thesis Other)

Understanding the transport potential of nearshore tsunami currents

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

UNDERSTANDING THE TRANSPORT POTENTIAL OF NEARSHORE TSUNAMI
CURRENTS
by
Aykut Ayca

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)

December 2018

ii
Acknowledgements
This dissertation is the final product of my 26 years long education life, and it has become
possible with the invaluable support of a lot of people. First and foremost, I have to express my
sincere gratitude to my advisor Prof. Dr. Patrick Lynett for sharing, generously, his guidance,
advice, criticism, and insight throughout my Ph.D. studies. He was a great role model and mentor.
This degree wouldn’t be possible without his support and motivation.
I would also like to thank my dissertation committee members Prof. Dr. J.J. Lee and Prof. Dr.
Julian Domaradzki, whom also served in my qualification exam committee along with Dr. Felipe
de Barros, for their valuable suggestions and guidance which definitely made this work better.
During my time in USC, I had the privilege to collaborate with Prof. Dr. Costas Synolakis, who
was also in my qualification exam committee. His enthusiasm, wisdom, and advice helped me
enormously in expanding my vision.
My lab mates (in order of appearance) Sangyoung Son, Haengsik Ko, Hoda El-Safty, Sasan
Tavakkol, Vassilios Skanavis, Luis Montoya, Kostis Douligeris and Ezgi Cinar, thank you for all
the fun, discussions and collaborations that we had together. I need to thank Adam Keen for turning
our long road trips and tedious field surveys into unforgettable memories with the help of loads of
beer and crappy food of course. And a special thanks to Nikos Kalligeris, I was so lucky to come
to LA with you in the first place. Thank you for being an amazing roommate during the most
difficult times away from home.
I am also grateful to my former advisor Prof. Dr. Ahmet Cevdet Yalciner who hired me as his
research assistant in METU and introduced me to tsunamis. My journey since then and this work
is the outcome of his support and guidance in the beginning. Since I mentioned METU, I am

iii
thankful to my former labmate and colleague Gokhan Guler for guiding me in setting up the
OpenFOAM simulations.
I am grateful for the received funding which helped with my financial survival during my Ph.D.
life in California. This research was supported mostly by funding provided by the California
Geological Survey and National Science Foundation.
I can’t express enough my gratitude to all of my family and family-in-law members in Turkey
who not only emotionally but also financially supported me during this journey. Especially my
mother and father encouraged and supported me to pursue my degree abroad, but even before that
they worked hard and sacrificed a lot to provide me with the best opportunities throughout my life.
And my little sister İpek; although I couldn’t be there for you as you grow, I was always proud of
you and felt your constant love and support here with me all the time. Every success I acquired is
also your achievement.
Last but not least, I would like to thank my wife Seda and my daughter Defne. I can’t imagine
reaching the finishing line without you. Seda, your presence made everything less burdensome,
more amusing and bearable together with your endless support, patience, understanding, and love
not only throughout this study but also in our life in a foreign country. And, my little sweetheart
Defne, my main motivation is to be a successful role model for you. My last year in the program
was extra challenging due to running between a Ph.D. and parenthood, but it is the most rewarding
of all. I am so blessed and lucky to have you.
My entire education in Turkey was free, or to be more precise, funded by the people of modern
Turkey. A country that is modern, democratic and independent built over the ruins of an empire
by one of the greatest revolutionists in the history, Mustafa Kemal Atatürk. Without him, his ideas

iv
and vision, and the country he founded, I couldn’t have come this far. Therefore, as a humble
payback for everything you had done for us; this work is dedicated to you.

v
Table of Contents
Table of Contents ................................................................................................................................................... v
List of Tables ....................................................................................................................................................... vii
List of Figures ..................................................................................................................................................... viii
Abstract .......................................................................................................................................................... xii
1. . Introduction ...................................................................................................................................................... 1
1.1 Aim of this Study ............................................................................................................................ 6
1.2 Numerical Model............................................................................................................................. 8
2. . Effects of Tides and Tsunami Source Directionality on Tsunami-Induced Currents .......................................... 11
2.1 Overview ...................................................................................................................................... 11
2.2 Tidal Time Series .......................................................................................................................... 16
2.3 Wave Directionality Analysis ........................................................................................................ 19
2.4 Model Results and Discussion ....................................................................................................... 21
2.4.1 Effect of Tides on Tsunami Currents ..................................................................................... 21
2.4.2 Wave Directionality Simulations: Maximum Water Surface Elevations and Current Speeds .. 25
2.5 Conclusions ................................................................................................................................... 31
3. . Sediment Transport ......................................................................................................................................... 33
3.1 Scalar Sediment Transport Model .................................................................................................. 35
3.2 Erosion and Deposition Rates ........................................................................................................ 36
3.3 Numerical Methods ....................................................................................................................... 38
3.4 Validation of the Model ................................................................................................................. 42
3.4.1 Laboratory Experiment on 2D Dam Break Flow ................................................................... 42
3.4.2 Large Scale Tests in Crescent City and Santa Cruz Harbors................................................... 44
3.5 Conclusions ................................................................................................................................... 50
4. . Modeling the Transport of the Large Vessels with the Effect of Tsunami Currents ........................................... 52
4.1 Inclusion of Ships as Pressure Disturbance ..................................................................................... 53
4.1.1 Waves Generated by Moving Pressure Disturbances ............................................................. 55
4.1.2 Benchmark Test in OpenFOAM ........................................................................................... 63
4.2 Equations of Motion ...................................................................................................................... 72
4.3 Collision Model ............................................................................................................................. 78
4.3.1 Collision Detection ............................................................................................................... 79
4.3.2 Collision Solver.................................................................................................................... 80
4.4 Large Scale Numerical Experiments .............................................................................................. 89
4.4.1 Hypothetical Tsunami Scenario in Port of Long Beach .......................................................... 91
4.4.2 Incident in Ishinomaki Harbor During 2011 Tohoku Tsunami ............................................... 95

vi
4.5 Conclusions ................................................................................................................................. 101
5. . Conclusions and Future Work Ideas .............................................................................................................. 104
5.1 Conclusions ................................................................................................................................. 104
5.2 Suggested Future Works .............................................................................................................. 105
References ........................................................................................................................................................ 107
Appendix ........................................................................................................................................................ 113

vii
List of Tables
TABLE 2-1: GRID SETUP FOR EACH OF THE STUDY SITES .......................................................................................... 15
TABLE 2-2: TIDE STATIONS USED IN THIS STUDY..................................................................................................... 16
TABLE 2-3: AMPLITUDE TO DEPTH RATIOS FOR TIDE AND TSUNAMIS IN ALL THREE STUDY SITES............................... 18
TABLE 3-1: NESTED GRID SYSTEM FOR EACH STUDY SITE. ....................................................................................... 45
TABLE 4-1: THE MAIN PARAMETERS OF THE VESSELS USED IN PORT OF LONG BEACH SIMULATIONS ......................... 93
TABLE 4-2: THE MAIN PARAMETERS OF C.S VICTORY............................................................................................. 96

viii
List of Figures
FIGURE 1-1 COMPARISON OF TIME HISTORIES OF (A) MEASURED (COLORED DOTS) VS. MODELED (CURVES) VELOCITIES
ALONG THE ENTRANCE CHANNEL OF VENTURA HARBOR, B) TIDE GAUGE RECORDED DURING THE TSUNAMI VS.
FREE SURFACE ELEVATION OBTAINED BY MOST (KALLIGERIS ET AL. 2016). .................................................. 10
FIGURE 2-1: LAYOUTS OF THE STUDY HARBORS A) CRESCENT CITY HARBOR, B) HALF MOON BAY AND C) SAN DIEGO
BAY. ........................................................................................................................................................... 15
FIGURE 2-2: A) INPUT TIME SERIES OF TIDES WITH RESPECT TO MEAN LOW WATER (MLW) LEVEL FOR EACH LOCATION
EXTRACTED FROM TIDE GAUGE STATION MEASUREMENTS; B) TSUNAMI TIME SERIES FROM TOHOKU 2011
TSUNAMI, PREDICTED BY MOST AT EACH LOCATION. THE STUDY SITES ARE SHOWN IN FIGURE 2-1. ............... 18
FIGURE 2-3: TIME SERIES FROM ALL THREE HYPOTHETICAL SOURCES PREDICTED BY MOST OFFSHORE PILLAR POINT
HARBOR SHOWN IN FIGURE 2-1. AS CAN BE SEEN, THE MAXIMUM WAVE AMPLITUDE IS 1 M FOR ALL. .............. 20
FIGURE 2-4: MAXIMUM CURRENT SPEEDS PREDICTED BY MOST FOR A) TSUNAMI ONLY, B) WHEN TSUNAMI ARRIVES
AT HIGH TIDE, C) WHEN TSUNAMI ARRIVES AT LOW TIDE, AND D) WHEN TSUNAMI ARRIVES AT MID-HIGH TIDE.
THE VECTORS SHOW THE DIRECTION OF THE MAXIMUM DETIDED CURRENTS. .................................................. 21
FIGURE 2-5: MAXIMUM CURRENT SPEEDS PREDICTED BY MOST FOR A) TIDE ONLY, B) WHEN TSUNAMI ARRIVES AT
HIGH TIDE, C) WHEN TSUNAMI ARRIVES AT LOW TIDE, AND D) WHEN TSUNAMI ARRIVES AT MID-HIGH TIDE. THE
VECTORS SHOW THE DIRECTION OF THE TIDAL CURRENTS AT THE TIME OF THE MAXIMUM CURRENTS OCCUR. .. 22
FIGURE 2-6: NORMALIZED MAXIMUM MEAN CURRENTS PLOTTED AGAINST TIDE LEVEL AT A) SAN DIEGO BAY, B)
HALF MOON BAY, AND C) CRESCENT CITY HARBOR. THE TIDE ELEVATIONS ARE ALSO SHOWN. ..................... 23
FIGURE 2-7: MAXIMUM COMPUTED TSUNAMI HEIGHTS ACROSS THE PACIFIC (LEFT PANEL), AND THE WAVES
APPROACHING THE CALIFORNIA (RIGHT PANEL) FROM A) ALASKA B) CHILE AND C) MARIANA SCENARIOS.
BLACK DOTS MARK THE LOCATION OF THE HALF MOON BAY. ....................................................................... 26
FIGURE 2-8: PROBABILITY DISTRIBUTIONS FROM ALL THREE SOURCES OF A) MAXIMUM CURRENT SPEEDS AND B) FREE
SURFACE ELEVATIONS. ................................................................................................................................. 27
FIGURE 2-9: SCATTER PLOT OF MAXIMUM SIMULATED CURRENT SPEEDS AS A FUNCTION OF WATER DEPTH OF ALL
THREE SCENARIOS. ....................................................................................................................................... 28

ix
FIGURE 2-10: CORRELATION OF THE MAXIMUM CURRENTS SPEEDS OF MARIANA SCENARIO VERSUS MAXIMUM
CURRENT SPEEDS FROM A) ALASKA SCENARIO AND B) CHILE SCENARIO PREDICTED AT EVERY GRID POINT IN THE
INNERMOST GRID. THE SOLID LINE INDICATES THE 1:1 RELATION AND THE TWO DASHED LINES INDICATE 1 M/S
ERROR BOUNDS. ........................................................................................................................................... 29
FIGURE 2-11: A) MAXIMUM CURRENT SPEEDS ESTIMATED BY THE MODEL IN PILLAR POINT HARBOR, B) ABSOLUTE
DIFFERENCE IN MAXIMUM CURRENT SPEEDS FROM ALASKA AND CHILE CASES WITH RESPECT TO MARIANA
CASE, AND C) RELATIVE (%) DIFFERENCE IN MAXIMUM CURRENT SPEEDS FROM ALASKA AND CHILE CASES WITH
RESPECT TO MARIANA CASE. ........................................................................................................................ 30
FIGURE 3-1: EXPERIMENTAL SETUP OF 2D DAM-BREAK FLOW THROUGH A PARTIAL BREACH. ................................... 42
FIGURE 3-2: COMPARISONS OF BOTTOM PROFILE CALCULATED BY MODEL TO EXPERIMENTAL DATA AT A) CROSS-
SECTION 1, B) CROSS-SECTION 2. .................................................................................................................. 43
FIGURE 3-3 COMPARISON OF BATHYMETRIC CHANGES IN CRESCENT CITY HARBOR OBTAINED BY A) MODEL WITH D50
= 0.15 MM, B) MODEL WITH D50 = 0.20 MM, C) MODEL WITH D50 = 0.25 MM, D) MODEL WITH D50 = 0.35 MM TO E)
FIELD DATA PROVIDED BY CGS. ................................................................................................................... 46
FIGURE 3-4: A) SCATTER OF MEASURED DATA VS. MODEL RESULTS IN CRESCENT CITY. RED LINES ARE THE ± 1
STANDARD DEVIATION FROM THE 1:1 REGRESSION LINE SHOWN WITH BLACK. YELLOW LINE INDICATES THE
MEAN OF THE DATA CALCULATED WITHIN 0.5-METER BINS. THE COLORBAR INDICATES THE DENSITY OF THE
POINT CLOUD; B) CHANGE IN MEAN OF THE DATA CALCULATED WITHIN 0.5 METER BINS WITH THE EMPIRICAL
PARAMETER 𝜑. ............................................................................................................................................ 47
FIGURE 3-5: COMPARISON OF BATHYMETRIC CHANGES IN THE ENTRANCE OF SANTA CRUZ HARBOR OBTAINED BY A)
MODEL WITH D50 = 0.15 MM, B) MODEL WITH D50 = 0.20 MM, C) MODEL WITH D50 = 0.25 MM, D) MODEL WITH
D50 = 0.35 MM TO E) FIELD MEASUREMENTS MADE BY CGS. ......................................................................... 49
FIGURE 4-1: THE SNAPSHOTS OF THE FREE SURFACE AT VARIOUS INSTANTS FOR 𝐹𝑟 = 1.0 AND 𝑝𝑚𝑎𝑥 =0.3. LEFT
PANEL IS TAKEN FROM ERTEKIN ET AL., (1986); RIGHT PANEL SHOWS THE RESULTS OBTAINED BY MOST. ...... 57
FIGURE 4-2: THE SNAPSHOTS OF THE FREE SURFACE AT VARIOUS INSTANTS FOR 𝐹𝑟 = 1.2 AND 𝑝𝑚𝑎𝑥 =0.1. LEFT
PANEL IS TAKEN FROM ERTEKIN ET AL., (1986); RIGHT PANEL SHOWS THE RESULTS OBTAINED BY MOST. ...... 58

x
FIGURE 4-3: COMPARISON OF FREE SURFACE ELEVATIONS CALCULATED IN MOST WITH ERTEKIN’S MODEL AT A)
65ℎ𝑜 AND B) 130ℎ𝑜 AHEAD OF THE INITIAL POSITION OF THE PRESSURE DISTURBANCE FOR 𝑝= 0.3 AND 𝐹𝑟=
1.0. ............................................................................................................................................................. 60
FIGURE 4-4: TIME SERIES COMPARISON OF FREE SURFACE ELEVATIONS AT DIFFERENT INSTANCES DURING THE
SIMULATION FOR 𝐹𝑟 =1.05 & 1.2 IN MOST, AND 𝐹𝑟 =1.2 IN ERTEKIN, FOR 𝑝= 0.1. .................................. 62
FIGURE 4-5: SCHEMATIC VIEW OF THE BENCHMARK SET-UP. A) TOP-VIEW, SHOWING THE BLOCK AND THE NUMERICAL
GAUGES WHERE THE TIME SERIES ARE EXTRACTED; B) SIDE VIEW OF THE SET-UP (VERTICALLY EXAGGERATED).
................................................................................................................................................................... 64
FIGURE 4-6: VELOCITY VECTOR MAPS AT VARIOUS TIMES DURING VORTEX GENERATION AND EVOLUTION. LEFT PANEL
MOST; RIGHT PANEL OPENFOAM RESULTS. ................................................................................................ 67
FIGURE 4-7: VORTICITY (𝜔) MAPS AT VARIOUS TIMES DURING VORTEX GENERATION AND EVOLUTION. THE RED LINES
REPRESENT THE 𝜔= 0 CONTOUR. LEFT PANEL MOST; RIGHT PANEL OPENFOAM RESULTS........................... 68
FIGURE 4-8: SWIRL STRENGTH (𝜆𝑐𝑖) MAPS AT VARIOUS TIMES DURING VORTEX GENERATION AND EVOLUTION. THE
RED LINES REPRESENT THE 𝜆𝑐𝑖 = 0 CONTOUR. LEFT PANEL MOST; RIGHT PANEL OPENFOAM RESULTS. ....... 69
FIGURE 4-9: COMPARISON OF THE TIME SERIES OF STREAMWISE VELOCITIES, 𝑢 (𝑚/𝑠), ESTIMATED IN MOST (BLUE
CURVES) AND OPENFOAM (RED CURVES), FOR GAUGES 1-8 (SEE FIGURE 4-5 FOR THE EXACT LOCATIONS OF THE
NUMERICAL GAUGES). .................................................................................................................................. 70
FIGURE 4-10: COMPARISON OF THE TIME SERIES OF STREAMWISE VELOCITIES, 𝑢 (𝑚/𝑠), ESTIMATED IN MOST (BLUE
CURVES) AND OPENFOAM (RED CURVES), FOR GAUGES 9-18 (SEE FIGURE 4-5 FOR THE EXACT LOCATIONS OF
THE NUMERICAL GAUGES). ........................................................................................................................... 71
FIGURE 4-11: NAME OF THE TRANSLATORY AND THE ANGULAR DISPLACEMENTS OF A SHIP, AS PRESENTED IN SNAME
(1950). ........................................................................................................................................................ 73
FIGURE 4-12: A) ELLIPSES 𝐸1 AND 𝐸2, AND THE MINIMUM DISTANCE, 𝑑𝑚𝑖𝑛 BETWEEN THEM AT TIME 𝑡0; B) TWO
BODIES COLLIDE AND INTERPENETRATION OCCURS AT TIME 𝑡0+Δ𝑡; C) THE POSITION OF THE BODIES AT EXACT
TIME OF COLLISION, 𝑡𝑐, WITHIN GIVEN DISTANCE TOLERANCE OF 𝜖. ................................................................ 80
FIGURE 4-13: A) BODY A AND BODY B APPROACHING EACH OTHER WITH PROSPECTIVE CONTACT POINTS 𝑝𝑎(𝑡) AND
𝑝𝑎𝑡; B) TWO BODIES ARE IN CONTACT AT 𝑡 =𝑡0, 𝑝𝑎𝑡0= 𝑝𝑏𝑡0. ................................................................... 81

xi
FIGURE 4-14: A) THE RELATIVE VELOCITY VECTOR AT THE TIME OF THE CONTACT; B) THE IMPULSE DUE TO COLLISION
ACTING ON EACH BODY; C) POST-IMPULSE RELATIVE VELOCITY VECTOR. ....................................................... 83
FIGURE 4-15: TIME LAPSE OF THE ECCENTRIC COLLISION TAKES PLACE BETWEEN TWO BODIES WITH EQUAL MASSES. 86
FIGURE 4-16: TIME LAPSE OF THE COLLISION TAKES PLACE BETWEEN TWO BODIES WITH DIFFERENT MASSES. ........... 87
FIGURE 4-17: TIME-LAPSE OF THE THREE-BODY COLLISION WITH EQUAL MASSES, AND THEIR CENTER OF GRAVITIES
ARE ON THE SAME LINE. THE BODY ON THE LEFT MOVES WITH THE INITIAL VELOCITY OF −0.1 𝑚/𝑠, WHILE THE
OTHER TWO ARE INITIALLY AT REST. ............................................................................................................. 88
FIGURE 4-18: TIME-LAPSE OF A THREE-BODY COLLISION WITH EQUAL MASSES. THE BODY ON THE LEFT MOVES WITH
THE INITIAL VELOCITY OF −0.1 𝑚/𝑠, WHILE THE OTHER TWO ARE INITIALLY AT REST. ................................... 88
FIGURE 4-19: AN EXAMPLE OF THE ELLIPTIC GAUSSIAN FUNCTION EMPLOYED, REPRESENTING A VESSEL WITH THE
LENGTH OF 250 METERS AND WIDTH OF 30 METERS. A) SHOWS THE PERSPECTIVE-VIEW IN 3D, B) SIDE-VIEW OF
THE PRESSURE DISTRIBUTION ....................................................................................................................... 91
FIGURE 4-20: SOUTHEAST BASIN OF THE PORT OF LONG BEACH. A) INITIAL LOCATION OF THE ONLY SHIP IN TEST CASE
ONE; B) INITIAL LOCATIONS OF ALL EIGHT IDENTICAL SHIPS ARE SHOWN. ....................................................... 93
FIGURE 4-21: CURRENT SPEEDS EXTRACTED 100 MINS AFTER THE ARRIVAL OF THE FIRST WAVE. RIGHT PANEL SHOWS
THE CURRENT FIELD FROM THE BASELINE SIMULATION; MIDDLE PANEL ARE THE CURRENT SPEEDS WHEN
VESSELS ARE INCLUDED; RIGHT PANEL ABSOLUTE DIFFERENCE OF CURRENT SPEEDS PREDICTED IN BOTH CASES.
THE FIGURES ON TOP ARE FROM ONE-SHIP SIMULATION, WHILE THE BOTTOM PLOTS ARE FROM THE MULTI-SHIP
CASE. .......................................................................................................................................................... 94
FIGURE 4-22: AERIAL PHOTO OF THE ISHINOMAKI PORT. THE RED CROSS SHOWS THE INITIAL LOCATION OF THE C.S.
VICTORY. THE RED LINE INDICATES THE APPROXIMATE PATH OF THE SHIP AFTER IT BROKE FREE. ................... 95
FIGURE 4-23: SNAPSHOTS OF THE CURRENT FIELD AND THE VESSEL’S LOCATION AT DIFFERENT INSTANTS. ............... 97
FIGURE 4-24: THE COMPARISON OF THE PATH OF THE VESSEL PREDICTED BY THE MODEL FOR THE DRAFT OF 7.5
METERS AND RELEASE TIME OF 66 MINS, WITH THE DATA PROVIDED BY DR. ASAI. .......................................... 98
FIGURE 4-25: THE PATHS OF THE VESSEL ESTIMATED BY THE MODEL FOR DIFFERENT VALUES OF A) RELEASE TIME; B)
DRAFT; C) INITIAL LOCATION; D) DRAG COEFFICIENT FOR LATERAL LOADING .............................................. 100
FIGURE 4-26: THE POINTS SHOW WHERE THE VESSEL IS ENDED UP AFTER 5 HOURS LONG SIMULATIONS. DIFFERENT
COLORS REPRESENT PARAMETERS THAT THE SENSITIVITY ANALYSIS IS CARRIED OVER.................................. 101

xii
Abstract
The purpose of this research is to study the transport processes induced by the tsunami-
generated currents in the nearshore. By this research, it is intended to significantly advance our
ability to provide quantitative guidance to ports and harbors in the development of resiliency and
recovery plans.
First, the effects of dynamic tides and wave directivity on maxima of tsunami-induced currents
were investigated. The goal of this part was to construct a basis for understanding the dynamic
effects of tides and wave directivity on current-based tsunami hazards in a coastal zone through
the application of numerical simulation tools for hazard mapping purposes. The effects of both
were quantified through sensitivity analysis. It was concluded that the tide induced variability in
maximum tsunami currents can be in the order of 25%, while wave directivity has negligible
impact.
Then, an integrated numerical tool was developed to predict the changes in the bottom
topography inside the harbor basins due to sediment transport driven by the tsunami-induced
currents. The recent tsunamis that were affected California also led to led to areas of extreme scour
as well as areas of large deposition caused by the strong currents. Thus, the main goal here was to
be able to provide quantitative guidance to the authorities and/or to stakeholders to develop counter
measures against this phenomenon through numerical modeling.
Finally, the motion of large vessels in ports under the effect of the tsunami-induced currents
was investigated in this thesis. The recent observations revealed that the maneuvering capabilities
of vessels can become very limited when they are captured by strong currents, and within eddies
in particular. The vessels that break free and drift-off in and around the port during tsunamis

xiii
generally exhibits random behaviors, and thus poses a great risk of hazard to the port and its
vicinity. The numerical tool was designed to allow for the interaction between flow and the vessel.
Since the draft of the vessel in shallow waters reduces to the effective flow depth considerably,
this interaction gets intensified when ships are in restricted waters like ports. It was shown that the
developed tool is be able to solve the interaction between the flow and the vessel reasonably The
integrated collision model can evaluate the colliding contact between multiple vessels or a vessel
and port structures in a physically correct manner as long as the vessels have elliptical shapes. In
general, the developed model can provide reasonable predictions on the ship’s path, yet a
sensitivity analysis conducted in Ishinomaki Port revealed that these predictions are extremely
sensitive to the initial conditions. This finding indicates that the process in fact is chaotic and
requires further assessmen

1
1. Introduction
Until recently, the focus of tsunami hazard studies has mostly been overland flooding,
inundation, and/or damage to coastal infrastructure due to the waves. However, the latest
transoceanic tsunamis have shown that even when there is no or little inundation, the currents
generated by tsunami surges can potentially cause significant damage to maritime facilities (Lynett
et al., 2012; 2014), raised new questions regarding the adverse effects of tsunamis. These unwanted
and often un-avoidable effects can range from small damage to vessels and infrastructure to the
complete destruction of harbors (Wilson et al. 2013; Dengler et al. 2008). The hydrodynamic
energy of these currents are also capable of entraining substantial amount of sediments from the
bottom, including meter-long boulders, at distances of several hundreds of meters to as far as a
few kilometers inland (Sugawara, et al. 2014). Significant changes in bottom morphology can
happen in relatively short time scales, which lead to decrease or interruptions in business activity
as well as added dredging costs (Tomita et al. 2006; Goto et al. 2011; Wilson et al. 2012).
Over the past few years, unfavorable nearshore effects of tsunami-induced currents from far-
field sources have been reported from many locations around the world. Several incidents had
occurred during the 2004 Indian Ocean tsunami, and were documented in detail by Okal et al.
(2006a; 2006b; 2006c) and Geist et al (2006). The most astounding of these events took place in
Port of Salalah, Oman, where an approximately 300-m long container ship broke all it’s twelve
mooring lines due to strong currents generated near the loading berth. The vessel drifted by the
currents, spinning uncontrollably for hours both in and around the harbor until it had been beached
over a sandbar by the flow. As reported in Okal et al. (2006a), the maximum tsunami amplitude in
Port of Salalah during the event was about 1.5-m and didn’t cause any overland flooding or

2
inundation. However, the currents generated locally by relatively small tsunami waves were strong
enough to pull a large container ship off the dock. The currents near the vessel were estimated to
be around 4-5 m/s by a numerical simulations of the event (Lynett et al., 2012). In the wake of the
Indian Ocean tsunami, a similar instance was also reported in Leport, Réunion Island. A 196-m
container ship, the MSC Uruguay, similarly broke all its mooring lines and drifted inside the
harbor, striking and damaging the cranes on a nearby dock (Okal et al. 2006b). Likewise, in Port
of Toamasina, Madagascar, a 50-m freighter wandered through the harbor for nearly 3 hours after
breaking its moorings by the strength of the currents (Okal et al., 2006c). All these incidents
occurred well after the arrival of the tsunami that was described as having maximum amplitude
and neither flooded the port facilities. This is important to consider, because practitioners without
specific experience in tsunami hydrodynamics may stop computations shortly after the maximum
runup of the first or second wave is reached. Eddies can persist for days, and this makes
computations more unreliable, as numerical errors accumulate.
The Tohoku tsunami that was triggered by the great Tohoku-oki earthquake (Mw = 9.0) was
documented quite well through comprehensive observations and instrumental measurements.
Although current speed measurements weren’t widespread along the Pacific Rim, there are still a
handful of data sets available. In the extreme near field, Fritz et al., (2012) analyzed the videos
recorded by two eyewitnesses in Kesennuma, Japan to determine the tsunami current speeds. Their
analysis revealed that the current speeds in Kesennuma Bay reached up to 11 m/s making
navigation impossible. However, in the far field, harbors and bays in Hawaii Islands (Cheung et
al., 2013), several harbors in New Zealand (Borrero et al., 2013; Lynett et al., 2012), and some
harbors in the Galapagos Islands, including the largest port of the country, Puerto Ayora (Lynett

3
et al., 2013), also experienced strong currents and damage related to these currents, even though
the incoming tsunami heights were relatively small.
Maritime communities along the U.S. West Coast also suffered due to currents generated by
the 2011 Tohoku tsunami. It caused adverse effects on almost every maritime facility along the
U.S. West Coast, ranging from interruptions of harbor operations to wide-spread damage to harbor
infrastructure. The strongest effects and most severe damage occurred in the Crescent City and
Santa Cruz Harbors (Wilson et al., 2012, 2013). Tsunami surges created very strong currents in
Crescent City’s inner harbor; video analysis indicated the current speeds of 4.5 m/s in the vicinity
of the entrance (Admire et. al., 2014). All docks and boats present at the time of the tsunami were
heavily damaged or completely destroyed. The damage in the harbor cost $20 M more on top of
the damage due to Kuril Islands tsunami in 2006. In Santa Cruz, fast-moving flows appeared
approximately 3 h after the initial arrival of the tsunami. Dozens of docks were damaged in
addition to 14 boats sunk and several others damaged. In addition to these two harbors, moderate
to minor damage to docks and boats happened in Ventura, Moss Landing, Morro Bay, San Diego
Bay, and some Bay Area harbors. During the 2011 Tohoku tsunami, the overall cost of damage
along the California coast, in two dozen of harbors, was over $50 M.
California harbors have suffered tsunami-induced currents from relatively smaller transoceanic
events as well (Dengler et al., 2008; Wilson et al. 2012, 2013). In 2006, an earthquake (Mw = 8.3)
occurred near the Central Kuril Islands and it was followed by a tsunami. The peak height of the
tsunami that arrived at Crescent City, California was about 1.8 m trough to peak, but it did not
cause any flooding since the largest waves coincided with the low tide. Nonetheless, strong
currents began with the arrival of the first waves and became more energetic over time, causing
severe damage. Docks located closest to the entrance of the inner harbor, where the currents were

4
the strongest (Dengler et al. 2008), recorded the highest damage levels, with the tsunami causing
$20 million in losses (Dengler & Uslu, 2011).
Another transoceanic event, the Maule, Chile tsunami of February 2010, generated strong and
damaging currents in California. Some docks were damaged in San Diego Bay near the entrance
of Shelter Island, 20 docks were damaged in the Ventura Keys in Ventura Harbor, and in Santa
Cruz Harbor two boats broke free from their moorings and caused minor damage in collisions with
other boats and harbor infrastructure. In addition, large charter boats left the harbor prior the arrival
of the tsunami and waited 6 h for reentry as strong currents made it difficult to navigate around the
entrance (Wilson et al. 2012).
In addition to the impacts of tsunami-induced currents discussed above, the sediment transport
induced by tsunamis within the harbors has long-term impacts on the recovery and resiliency of
the affected maritime communities. This was highlighted after the 2011 Tohoku tsunami, by the
observations in the near field (Weiss & Bourgeois, 2012), and also in the far field; particularly
along the California coast in harbors like Crescent City and Santa Cruz (Wilson et al. 2012).
Earlier studies have reported issues related to the tsunami-induced sediment transport
worldwide. During the 1960 Chile Tsunami, a jetty in Kesennuma Bay Japan collapsed due to
scour depths of 10 meters near its footing (Takahashi et al, 2001). While, in the aftermath of the
2004 Indian Ocean tsunami, in Kirinda Harbor Sri Lanka, Goto et al. (2011) evaluated the
bathymetric changes where the maximum thickness of sediment deposits was on the order of 4 m
within the harbor, the sediment deposits of 0.2 meters height were found 250 meters inland during
the field survey in American Samoa following the 2009 South Pacific Tsunami (Apostos et al,
2011).

5
California harbors had also suffered multiple times due to sediment transport caused by the
transpacific tsunamis. After the 1960 Chilean and 1964 Alaskan tsunamis, scour and sediment
deposition in some California harbors was reported by Lander et al., (1993). Approximately 4 m
of sediment deposition was reported in the outer harbor of Crescent City during the 1960 tsunami.
In 1964, a number of harbors experienced scour around pilings and sedimentation within access
channels (NGDC, 2011). More recently, the Chilean tsunami of February 27, 2010, caused enough
erosion in the channel entrance to Ventura Harbor that the harbormaster reported $100,000 of
savings in dredging costs (Wilson et al., 2010).
Following the 2011 Tohoku tsunami, the survey organized by California Geological Survey
(CGS) showed that significant channel scour and sediment deposition occurred in Crescent City
and Santa Cruz harbors (Wilson et al., 2012). In the entire Crescent City harbor, about 289,400 m
3

of sediment was removed in an area of 0.67 km
2
, while the minimum deposition volume of 154,600
m
3
was estimated with the sediment covering 55% of the portion of the harbor contained in the
surveys. At the entrance of Santa Cruz harbor, measured erosion volumes ranged from 2550 to
14,800 m
3
, accompanied with fill estimates up to 8750 m
3
, depending upon the surveys used to
describe post-tsunami conditions. It took about a year to remove the significant amount of sediment
deposited in Crescent City and Santa Cruz harbors which extended the recovery duration and
caused interruptions in harbor operations.
In general, the extensive research on sediment deposits starting from earlier tsunamis to today
revealed that the tsunami-induced sediment transport can lead to significant changes in the
geomorphology in coastal areas (Dawson & Shaozhong, 2000; Scheffers & Kelletat, 2003).

6
1.1 Aim of this Study
The issues related to tsunami-induced currents in harbors as summarized above motivated this
proposed research. The first goal of this work is to investigate how the maxima of tsunami-induced
currents vary due to dynamic effects of tides and wave directivity. The hope is to construct a basis
for understanding the dynamic effects of tides and wave directivity on current-based tsunami
hazards in a coastal zone through the application of numerical simulation tools for hazard mapping
purposes. The consideration of these aspects is crucial and yet challenging in the modeling of
tsunami currents. A sensitivity analysis was conducted in three harbors by coupling the tsunami
with the tide signal at twelve different tide levels and allowing the initial tsunami wave to arrive
at various tidal phases. Then, to evaluate the effect of wave directionality on maximum currents,
three earthquakes with different magnitudes were devised along the Pacific, which were also tuned
to create the same maximum near-harbor amplitude. The motivation for testing this hypothesis is
the expectation that decisions about advisories or warnings from a distant source are based
primarily on the local height predictions from National Oceanic and Atmospheric Administration
(NOAA) and not yet on maximum nearshore current speeds.
The next aim is to develop an integrated numerical model to estimate how the bottom
morphology inside the harbor basins changes by the tsunami-induced sediment transport. As
described above, there are numerous instances when strong tsunami currents in harbors led to areas
of extreme scour as well as areas of large deposition. Numerical models are useful in providing
quantitative guidance to the authorities and/or to stakeholders to develop resiliency options. With
the insight gained from numerical hindcast studies, areas of high scour and deposition can be
predicted and plans to over or under-dredge areas during regular maintenance may be
implemented, which ultimately can be used to mitigate practical concerns like scouring of

7
structures or deposition along the waterways of ports and harbors. To achieve these goals, an
advection-diffusion type sediment transport equation has been solved using a finite volume
scheme, which greatly enhances the numerical stability, and closed with semi-empirical formulae
calculating the erosion and deposition rates. Here, only the transport of the suspended particles
was considered in this research, since it is the governing mode of sediment transport in an energetic
tsunami. However, one can easily include the other modes of transport, such as involving the bed
load transport, in the future if desired. The developed model was tested in the laboratory scale
based on an experiment investigated the scour in the downstream due to dam break, and in the
geophysical scale using the data collected in Crescent City and Santa Cruz Harbors after the 2011
Tohoku Tsunami. The results from these validation tests and the model performance will be
discussed in Chapter 3. One final note on the sediment transport model; as the hydrodynamic
solver, the Method of Splitting Tsunamis (MOST) was utilized here, yet the developed sediment
model can be coupled with any two-dimensional hydrodynamic model with only a few
modifications.
The last research objective of this study is to develop a numerical tool to predict the motion of
large vessels under the influence of tsunami waves and currents inside ports. While the literature
is lavish on the maneuverability of ships in the presence of harmonic waves, there aren’t many
studies focusing on tsunami-induced ship motions. This endeavor requires an accurate
representation of flow field around the vessels, since the interaction between the flow and the deep
draft vessels gets amplified in restricted waters. For this purpose, the momentum equations of the
hydrodynamic model were modified with an additional horizontal pressure distribution term to
account for the effects of floating bodies which allows for the interaction of the flow and the
vessels. The large vessel transport model solves a suite of equations which considers inertial,

8
damping and viscous forces acting on a large vessel due to its’ movement as well as the forces
exerted by the flow. However, this research does not concentrate on how to solve the different sub-
problems related the theory of motion of ships such as the determination of hydrodynamic
coefficients associated with the inertial or damping forces, but aims to understand the extent of
interaction between the tsunami and large vessels and to clarify how the effect of tsunami currents
should be included in this analysis.
Understanding how the ships behave during a tsunami, where waves are relatively small, but
the currents are strong and erratic, can be used to mitigate some of the concerns regarding the civil
defense policies of port facilities. Ultimately, it is desired to significantly advance our ability to
provide quantitative guidance to ports and harbors in the development of resiliency and recovery
plans.
1.2 Numerical Model
The hydrodynamic modeling results presented here come from the application of the “method
of splitting tsunami” (MOST) numerical model (Titov & Synolakis, 1998). The MOST model has
been used extensively for tsunami hazard assessments in the United States and is currently used
for operational tsunami forecasting at the NOAA Pacific Marine Environmental Laboratory
(PMEL). Variants of the MOST model have been in constant use for tsunami hazard assessments
in California since the mid-1990s. The MOST solves the classical 2+1 nonlinear shallow water
(NSW) equations using a finite difference scheme.
ℎ
I
+(𝑢ℎ)
J
+(𝑣ℎ)
L
= 0 1-1
𝑢
I
+𝑢𝑢
J
+𝑣𝑢
L
+𝑔ℎ
J
=𝑔𝑑
J
−𝐷𝑢 1-2

9
𝑣
I
+𝑢𝑣
J
+𝑣𝑣
L
+𝑔ℎ
L
=𝑔𝑑
L
−𝐷𝑣 1-3
Where 𝜂(𝑥,𝑦,𝑡) = wave amplitude, d is the water depth, ℎ(𝑥,𝑦,𝑡) =𝜂(𝑥,𝑦,𝑡)+𝑑(𝑥,𝑦,𝑡) ,
𝑢(𝑥,𝑦,𝑡) and 𝑣(𝑥,𝑦,𝑡) are the depth-averaged velocities 𝐷(ℎ,𝑢,𝑣) is the drag coefficient
computed by the quadratic equation. More thorough information about the theoretical background
and the validation of MOST is provided in Titov and Synolakis (1998), and Burwell et al., (2007).

𝐷(ℎ,𝑢,𝑣) =𝑛
Q
𝑔ℎ
R
S
T
U
V𝑢
Q
+𝑣
Q

1-4
Where n is the Manning’s roughness coefficient.
As a final remark about the verification of the MOST results, (Lynett et al., 2014) conducted a
sensitivity analysis to understand the accuracy of MOST in resolving complex nearshore
hydrodynamics, where they have compared the MOST results with a high-order fully nonlinear
Boussinesq-type model (COULWAVE) and with available field data. They found out that MOST
matches the tsunami amplitude phase patterns for the first several hours and the amplitude
envelope for at least 24 hours after the first arrival of the tsunami. Likewise, MOST can also match
the current velocities extracted from eyewitness videos. Also, on average, current velocities
predicted by MOST are about 10-20 % higher than COULWAVE results, which leaves MOST on
the conservative side. As a result of these observations, they concluded that MOST was sufficiently
capable of modeling tsunami-induced currents using a 1/3 arcsec (~10 m) grid size. (Kalligeris et
al., 2016) made in situ measurements of current speeds in Ventura Harbor California during the
2015 Chilean tsunami. In their comparison of MOST results to measured data, they found out that
MOST in general matches the phase of the measurement well, and the velocity estimates by MOST
remain in the maximum envelope of the measurements

10
Figure 1-1: Comparison of time histories of (a) measured (colored dots) vs. modeled (curves)
velocities along the entrance channel of Ventura Harbor, b) Tide gauge recorded during the
tsunami vs. free surface elevation obtained by MOST (Kalligeris et al. 2016).

a)
b)

11
2. Effects of Tides and Tsunami Source Directionality on Tsunami-
Induced Currents
2.1 Overview
There are several factors that can influence the variability of the maximum tsunami-induced
currents in harbors, including the tsunami’s source, distance to the target area, and nearshore
features. In addition, other key considerations for obtaining a detailed and accurate description of
the maximum tsunami currents can be the tide levels and their interaction with the tsunami.
Only a few studies have investigated tide–tsunami interactions, with a focus on understanding
their influence on the maximum wave heights, runup, and inundation limits. This has mainly been
because, until recently, the context of tsunami hazards has been limited to overland flooding and
inundation. An example of these previous studies is the series of numerical experiments on
tsunami–tide interactions in Cook Inlet, Alaska, conducted by Kowalik & Proshutinsky, (2010),
which involved very strong tides and frequent tsunamis. They selected Anchor Point and
Anchorage as their pilot study areas and found that, at Anchor Point, tsunamis are likely to be
amplified by a factor of 4–6 relative to the tsunami magnitude, depending on the tide level when
the tsunami arrives. In addition, at Anchorage, a tsunami could be damped by 5% or amplified by
35%, depending on the phase of the tide when the tsunami arrived. Yet, they did note that these
results were very site and event specific and could not be generalized. Similarly, (Zhang et al.,
2011) investigated the dynamic effects of tides during the 1964 Alaska tsunami. Their results
revealed that the tsunami–tide interaction has considerable influences on the resulting run-up and
inundation. They also suggested that the inclusion of tides is important to accurately predict coastal
and, especially, estuarine inundation during tsunamis.

12
Mofjeld et al., (2007) performed a probabilistic analysis of the effects of tides on the maximum
tsunami wave heights at Seaside, Oregon. Their approach was based on an analytically analysis
of linearly superimposed exponentially decaying tsunamis and predicted the tides on the open coast
at Seaside. Although their study did not include nearshore modeling, they still determined that the
means, standard deviations, and probability distribution functions of the maximum tsunami heights
varied with the tides at the offshore region of Seaside.
Unlike the elevation-focused studies discussed above, a recent study by (Lee et al., (2015)
mainly focuses on effects of dynamic tides on tsunami propagation, while Shelby et al., (2016)
studies tide-tsunami interaction in Hudson River Estuary, where the effects of tide-tsunami
interaction on flow velocities are briefly investigated, but only for limited number of tide phases.
Nonetheless, their results clearly reveal that large differences can be observed in maximum current
speeds due to nonlinear interaction between tsunamis and tides.
In light of these previous studies, the aim here is to place another missing piece of this puzzle
and explore how the dynamic effects of tides influence tsunami-induced currents inside harbors
and bays. It is anticipated that the tidal influence on the maximum currents in a harbor or port
during a tsunami could be significant when the tidal elevation and currents are similar in magnitude
to those of the tsunami, and the interaction between the tides and tsunami-induced currents can be
non-linear. Therefore, an adequate demonstration of the physics of this problem is of paramount
importance.
While one focus of this research is on the relevance of tides on the current field, the influence
of the source location of the tsunami will also be examined. In their recent study, Borrero et al.
(2015) conducted a comprehensive sensitivity analysis of tsunami current hazards as a function of
the tsunami source location. They placed equivalent tsunami sources at equally spaced azimuthal

13
directions around the Pacific Rim, and simulated each source for selected New Zealand harbors,
and assessed which source was responsible for maximum currents at each grid node in modeled
harbors. Their findings suggest that the tsunami currents in harbors are more sensitive to the source
location than the tsunami heights. Additionally, they also inspected spectral properties of incoming
tsunamis in correspondence to harbor response for one of their modeled harbors.
This research approaches the same idea from a different perspective. Here, it is attempted to
determine whether the source location of the tsunami plays a major role in the predicted maximum
currents inside a harbor for a specific predicted local wave height, i.e. if two different sources
produce the same local wave height, do they also produce similar currents? The consideration of
this hypothesis arose from inspection of recent tsunami current measurements inside harbors, and
the observation that often the maximum currents occurred both early in the event and near the
same time as the maximum observed amplitude. These recent measurements include those in Santa
Cruz, California, where currents driven by the 2010 Chile tsunami were recorded (Lacy et al.,
2012). During the 2011 tsunami, far field current data was measured in Tauranga Harbor, New
Zealand (Borrero et al., 2015), Port Ayora in the Galapagos (Lynett et al., 2012), and in Hilo,
Kahului, and Honolulu harbors in Hawaii (Cheung et a., 2013). In all of these measurements, the
maximum de-tided tsunami currents are associated with either the first, second, or third waves,
and so are not likely to be a strong function of local resonance, as this process should require many
cycles to become dominant. Further motivation for this approach is the expectation that decisions
about advisories or warnings from a distant source are based primarily on the local height
predictions from National Oceanic and Atmospheric Administration (NOAA), which already
considers the open ocean propagation. This is important in issuing advisories/warnings for harbors
regarding maximum currents to be expected, to take necessary mitigation measures and to continue

14
harbor operations if possible, considering the direct relation between the current speeds and the
observed damage is considered (Lynett et al., 2014). Here, a foundation is being laid down for
understanding the dynamic effects of tides and wave directivity on current-based tsunami hazards
in a coastal zone by the application of numerical simulation tools for hazard mapping purposes.
In this study, MOST was used to model the tsunami waves from their source to a nearshore
region through nested grids as mentioned before. The model propagated the tsunami waves to the
shore, computing the wave amplitude, velocity, and overland inundation and captures the non-
linearity of the waves as they reach shallow water. Computations were stopped at the 5-m-depth
contour in all the parent grids, at which depth waves were reflected back (a solid wall boundary
condition was imposed). Runup and inundation computations were only performed in the
innermost grid, in which a bottom friction term was included in the momentum equations. For all
simulations presented here, the Manning’s “n” friction factor was 0.03.
The outermost grid at 4-arcmin (~7200 m) resolution covered the entire Pacific basin. Three
additional grids of increasingly finer resolution were derived from data obtained from NOAA’s
freely available National Geophysical Data Center [ngdc.noaa.gov/mgg/inundation/tsunami/],
specifically for tsunami forecasting and modeling efforts. The innermost nearshore grid with the
highest resolution was at the 1/3 arcsec grid, with boundary inputs for free surface elevation and
velocities from the previous MOST nested layer. Details of computational grids for all locations
are listed in Table 2-1. Also, Figure 2-1 maps the layout of all three-study areas in which
investigated harbors and tide stations are marked.

15

Figure 2-1: Layouts of the study harbors a) Crescent City Harbor, b) Half Moon Bay and c) San
Diego Bay.

Table 2-1: Grid setup for each of the study sites
Lon. Range (
o
E) Lat Range (
o
N) Nx Ny dx dt (sec)
Propagation 120 – 292 -74 – 62 2581 2879 4’ 10
Crescent City
Level 1 234.2592 – 236.3725 39.0350 – 41.9884 318 444 24” 2.5
Level 2 235.5617 – 236.1183 40.5242 – 41.9542 168 430 12” 1.5
Level 3 235.7142 – 235.9133 41.6284 – 41.8150 240 225 3” 1.5
Level 4 235.7655 – 235.8569 41.7165 – 41.7829 988 718 1/3” 0.25
Pillar Point
Harbor

Level 1 236.0308 – 238.5392 36.0200 – 38.9117 302 348 30” 2.5
Level 2 236.7758 – 238.4258 36.4333 – 38.5333 397 505 15” 1.5
Level 3 237.4396 – 237.5796 37.3896 – 37.5396 169 181 3” 1.0
Level 4 237.4600 – 237.5500 37.4600 – 37.5200 973 649 1/3” 0.25
San Diego
Bay

Level 1 241.0342 – 242.9692 32.3908 – 34.1509 383 353 18” 2.5
Level 2 242.1475 – 242.9450 32.4683 – 33.5408 320 430 9” 1.5
Level 3 242.6450 – 242.9142 32.5708 – 32.9317 324 434 3” 1.0
Level 4 242.7333 – 242.9037 32.5845 – 32.7500 1842 1788 1/3” 0.25

16
2.2 Tidal Time Series
Of particular interest in the present study was determining whether the tsunami-induced
currents in ports/harbors/bays are modulated by the existing background tidal currents, and if so,
the extent of this modulation. To answer these questions, sensitivity analysis results were compiled
from the following selected study areas: San Diego Bay, Crescent City Harbor, and Pillar Point
Harbor located inside Half Moon Bay. Several factors were considered when selecting these study
sites. First, we wanted to examine locations in Northern, Central, and Southern California, where
the tidal characteristics are modestly different. In addition, the physical characteristics of
individual harbors or bays can play major roles in amplifying the effects of tides. In this regard,
San Diego Bay and Pillar Point Harbor represent different types of basins.
San Diego Bay is a large, long, and narrow bay, which contains different oscillatory modes,
whereas Pillar Point Harbor is a small rectangular-shaped boat harbor with less complex internal
hydrodynamics. Finally, in Northern California, Crescent City Harbor was a natural candidate for
this study, because the most severe effects of the 2011 Tohoku tsunami were observed there.
Additionally, several past events have caused significant damage in Crescent City Harbor (Dengler
et al., 2008, Wilson et al., 2012, 2013). The tide data used for each study area were collected from
the NOAA Tides and Currents database, where 1-min (or 6-min, depending on the type of tide
station) water level data recorded during the 2011 Tohoku tsunami by local tide stations are
available for download. The details about the tide stations used are listed in Table 2-2.
Table 2-2: Tide stations used in this study.

Station Name Station ID Latitude (
o
N) Longitude (
o
W)
Crescent City 9419750 41.7450 235.8150
Pillar Point Harbor 9414131 37.5017 237.5183
San Diego Bay 9410170 32.7133 242.8267

17
Figure 2-2a shows the tidal time series used in this study to force the innermost grids, which
were recorded between March 11, 2011, and March 15, 2011, during the Tohoku-Oki Tsunami.
The very first step in the analysis part of this work was to reproduce the tide signal numerically
using MOST at the location of the tide gauge station. For this purpose, the innermost grid was only
forced by tides from the boundaries, and these input files were modified and re-run until the
numerically predicted tides matched the data measured at the tide gauge station. Once this was
achieved, the tsunami (Figure 2-2b) and tidal time series (Figure 2-2a) obtained numerically from
the MOST simulations were linearly superimposed. The linear superposition of tides and tsunamis
is valid in the open ocean since their amplitudes are very small compared to the water depth (Dean
& Dalrymple, 1991; Shelby et al., 2016). However, when the linear superposition of these two
waves is in the near field, as in this study, careful examination is important. The amplitude to depth
ratios for tide and tsunami should be less than or close to 0.1 for linear superposition to be valid,
where the shallow water nonlinearity parameter is used as a proxy for the validity of linear
superposition. In Table 2-3, the average depth where the tide and tsunami are combined is shown,
as are amplitude-to-depth ratios in each study site. As Table 2-3 reveals, tides in all three sites and
tsunamis in Half Moon and San Diego Bays are clearly linear at the superposition location. Only
in Crescent City is the incident tsunami amplitude weakly nonlinear. Here, the average depth is
relatively small due to shallow areas in the northernmost part of its eastern boundary, which puts
the tsunami in the weakly nonlinear regime for this part of the boundary. Nearer to the harbor
entrance, the depths at the superposition location are greater than 50 meters. Thus, it can be
assumed that the weakly nonlinear effects that might arise from a linear superposition of tide and
tsunami signals offshore of the harbor entrance are small.

18
Table 2-3: Amplitude to depth ratios for tide and tsunamis in all three study sites.
Crescent City Half Moon Bay San Diego Bay
Average Depth of Superposition 25 m 45 m 35 m
Tide only 0.064 0.027 0.034
Tsunami only 0.129 0.018 0.014

Figure 2-2: a) Input time series of tides with respect to mean low water (MLW) level for each
location extracted from tide gauge station measurements; b) Tsunami time series from Tohoku
2011 tsunami, predicted by MOST at each location. The study sites are shown in Figure 2-1.

Here, it must be noted that Fig. 2-2b shows a sample time series of the tsunami recorded at one
grid to illustrate the tsunami characteristics; however, MOST produces time series of free surface
elevation and flow velocities at each grid node along the seaward boundaries of the innermost
grids. A single tidal time series is superimposed with the tsunami since the extent of the innermost
boundaries were only a few kilometers, for which the spatial variation of tidal amplitudes will be
very small. Then, the tsunami signal is shifted in time to see how the tsunami currents would have

19
been affected if the tsunami had arrived at the maximum tide, minimum tide, or any other
intermediate tide level. For each study area, tsunamis were superimposed at 12 different tide
phases over a tide cycle. This illustrated the influence of the tidal currents on the tsunami currents
for various tsunami arrival times.
2.3 Wave Directionality Analysis
The approach for this analysis was to run hypothetical tsunami scenarios from three different
source regions along the Pacific and tune them in a way to ensure that they would create the same
maximum offshore wave amplitude near the selected study area. Then, the resulting current fields
could be compared. Furthermore, it is well known that when examining wave-driven free surface
flows, small changes (or errors) in the sea surface amplitude can lead to large changes (or errors)
in the fluid speed (Lynett et al., 2012, Borrero et al., 2015). Therefore, it was crucial in this study
to obtain the same offshore maximum free surface elevation within ±1% precision near the study
area for all three sources, so that their resulting current fields can be comparable.
Pillar Point Harbor, located inside Half Moon Bay, was picked as the pilot study area. Pillar
Point is a south-facing harbor located along the central coast of California and is exposed to the
effects of tsunamis from all over the Pacific Rim. The harbor, with its south-facing orientation,
should represent a difficult location to demonstrate that wave directionality is not an influential
parameter. The sources devised for this study were assumed to arise from the Alaskan–Aleutians
Subduction Zone, the Chile Subduction Zone, and the Mariana Subduction Zone. It must be noted
that these sources are completely hypothetical, or, in other words, they have no physical basis.
These sources are tuned to yield the maximum free surface elevation of 1 m offshore Pillar Point
Harbor. The point where the water surface elevations are compared and time series from three
different sources can be seen in Figure 2-3.

20

Figure 2-3: Time series from all three hypothetical sources predicted by MOST offshore Pillar
Point Harbor shown in Figure 2-1. As can be seen, the maximum wave amplitude is 1 m for all.

The model scenarios were devised by combining the NOAA-PMEL’s pre-defined unit sources
placed on fault plane segments positioned along the world’s subduction zones. Each unit source is
100 km long while 50 km wide and has 1-meter displacement that produces a Mw 7.5 earthquake.
Scaling and/or combining these unit sources, one can create more complex faulting scenarios.
Further information about the pre-defined unit sources in the Pacific can be found in (Gica et al.,
2010)
The Alaska scenario is the combination of 12 unit sources and represents a Mw ~ 8.6
earthquake with a uniform slip of 3.75 meters, and a rupture area of 1200 km by 100 km in length
and width respectively. The Mariana scenario, which is an Mw ~ 8.9 earthquake, is 1600 km long
and 100 km wide and has 8 meters of uniform slip. Lastly, the Chile scenario corresponds to an
Mw ~ 9.2 earthquake, and 1400 km in length while 100 km in width, with the uniform slip of 23
meters. Tables A1 to Table A3 provide all the details of the source parameters used for each
scenario, each creates a 1-meter amplitude near the study area.

21
2.4 Model Results and Discussion
2.4.1 Effect of Tides on Tsunami Currents
Figures 2-4 and 2-5 show model results of the simulation-based information, which should be
useful for understanding the effects of tides on tsunami currents.

Figure 2-4: Maximum current speeds predicted by MOST for a) Tsunami only, b) when tsunami
arrives at high tide, c) when tsunami arrives at low tide, and d) when tsunami arrives at mid-high
tide. The vectors show the direction of the maximum detided currents.

When interpreting these results, the following should be kept in mind: the results presented
here were de-tided, which means the velocity time series of the tide-only simulations were
subtracted from the time series obtained from the tsunami-plus-tide simulations. This yields a time

22
series that includes the tsunami signal and any alteration of this signal from the tide-tsunami
interaction. Figure 2-4, plots the maximum current speeds at Crescent City, predicted from the
2011 Tohoku tsunami source, superimposed with the tide signals for high, low, and mid-high tides,
along with the results for the tsunami-only case. Although it is difficult to quantify the effects of
the tides on the maximum currents from Figure 2-4, the differences in the overall current fields
compared to the tsunami-only and various tsunami-plus-tide cases is evident. This is more
profound near the entrance of the outer harbor as can be seen in Figure 2-5, which is almost the
same plot as Figure 2-4 except the vectors showing the direction of the tidal currents.

Figure 2-5: Maximum current speeds predicted by MOST for a) Tide only, b) when tsunami arrives
at high tide, c) when tsunami arrives at low tide, and d) when tsunami arrives at mid-high tide. The
vectors show the direction of the tidal currents at the time of the maximum currents occur.

23
In terms of quantification of the tide effects on currents, one can refer to Figure 2-6, in which
the normalized means of the predicted maximum currents at each location are plotted against the
corresponding tide level. The mean of the maximum currents is the spatial average of the maximum
detided currents calculated at each grid node. Subsequently, they were normalized by the mean
maximum currents obtained from tsunami + max tide case. For example, the value corresponding
to “Max” in Figure 2-6 shows the mean of the maximum currents estimated in that particular bay
or harbor, if the tsunami would arrive during the high tide. The same logic applies for “Min,”
“Mid-High,” “Mid-Low,” or for any other intermediate tide level. This metric makes it possible to
estimate the scale of the tidal effects on tsunami currents.

Figure 2-6: Normalized maximum mean currents plotted against tide level at a) San Diego Bay, b)
Half Moon Bay, and c) Crescent City Harbor. The tide elevations are also shown.

24
Figure 2-6a shows a plot of the variation in the mean of the maximum speeds in San Diego
Bay. The focus was on two harbors: the America’s Cup Harbor located in the north, and the Yacht
Basin in the south of Shelter Island (Figure 2-1c), in which the strongest currents were witnessed
during the 2011 Tohoku tsunami (Wilson et al., 2013). The highest modulation in the mean-
maximum currents due to tides was estimated in the America’s Cup Harbor, where the difference
between the lowest and highest mean-maximum currents was over 25%, with the same metric
showing a 15% change in the Yacht Basin. However, these maximum tide-tsunami interaction
effects in the two harbors occur at very different tidal phases, indicating that the effect is highly
localized, even with two harbors in close proximity within the same bay system.
In Crescent City Harbor, the difference between the smallest and largest mean-maximum
currents is on the order of 20% (Figure 2-6c). It could be expected that greater effects from tides
would be observed here, as the tidal range is larger; however, the tsunami amplitude is also larger,
and these two scales compete against each other.
Finally, in Pillar Point Harbor, this margin from largest to smallest is around 10%, as shown
in Figure 2-6b, which agrees well with a priori considerations about the scale of the tidal effects
on the currents in Pillar Point. Looking across the various scaled mean-maximum speed profiles,
there is no clear correlation between a specific tidal phase and the creation of the greatest currents.
Thus, this comparison sheds some insight into the oft-asked question of when during the tidal cycle
should one expect the tsunami currents to be greatest: it appears there is no good answer as local,
site-specific effects control the response. However, what this comparison does show rather clearly
is that for a tsunami with height on the order of the tidal range, the variability of tsunami arrival
time with respect to the tidal phase is likely to impact the maximum currents with uncertainty in
the range of +/- 25% or less.

25
2.4.2 Wave Directionality Simulations: Maximum Water Surface Elevations and Current
Speeds
The results obtained in Pillar Point Harbor from the hypothetical tsunami simulations are
presented in Figures 2-7 to 2-11. Maximum computed tsunami wave heights for the propagation
grid across the Pacific, and waves propagating towards the Central California are shown in Figure
2-7, which gives an idea about the initial directions of all scenario tsunamis relative to the Half
Moon Bay. It can be seen that the three sources are approaching Pillar Point harbor from a spread
of roughly 180 degrees, and thus these sources do test the limits of offshore directionality. Note,
however, that while in the offshore region the wave fronts show large differences in angles of
incidence, by the time these wave fronts have refracted into the shallow water offshore of the
harbor, these differences are much reduced. It is indeed this observation that motivated this part of
the study.
The probability distribution of the maximum speeds given in Figure 2-8a suggests a good
agreement between all three cases. The peaks of the curves lie at around 0.7 m/s, and all three
distribution shapes are similar and are characterized by slowly decaying tails at the high-velocity
end. On the other hand, the probability distributions of the free surface elevations display arguably
more source-dependent variability, with peaks in the range of 0.6–0.8 m. However, it must be
reiterated that the velocity distributions contain statistically significant probability values
extending beyond seven times the peak value of 0.7 m/s, while free surface elevation distributions
only extend 0.5-1.0 beyond the peak. The velocity distribution for a tsunami event is expected to
be much broader, in a scaled sense, than its free surface elevation distribution. This has clear
implications for uncertainty quantification in tsunami speed predictions.

26

Figure 2-7: Maximum computed tsunami heights across the Pacific (left panel), and the waves
approaching the California (right panel) from a) Alaska b) Chile and c) Mariana scenarios. Black
dots mark the location of the Half Moon Bay.

27

Figure 2-8: Probability distributions from all three sources of a) maximum current speeds and b)
free surface elevations.

Taking the model data from all three scenarios, a scatter plot of maximum currents versus
water depth as shown in Figure 2-9. From the figure, it is seen that the largest spread in simulated
maximum currents exist at depths of 5-10 m, which is due to large eddies and jets coming in and
out from the outer harbor. Despite this large physical variability, the envelope of maximum
currents from all scenarios shows good agreement with each other. At depths greater than 15 m
the current speeds show less variability in general, and similar for all cases. The single point peaks
in current speeds (speeds higher than 4 m/s) observed at depths less than 5 m are due to model
instabilities which happen at locations like breakwater tips, sharp corners etc. If these points are

28
neglected (they would typically be filtered out when producing hazard maps), the maximum
currents lay in the same range for this depth region as well. The differences in the predicted
maximum current speeds for each scenario are examined further in Figure 2-10. In this figure, the
simulated maximum currents from the Alaska and Chile scenarios are plotted against Mariana
scenario for every grid point in the computation domain. The figure reveals that the correlation
between the maximum current speeds gets higher as the observed speed increases. For both Alaska
and Chile scenarios, the more than 98% of the grid points fall within +/- 1 m/s/ of the current
speeds observed in Mariana scenario, for all current speeds.

Figure 2-9: Scatter plot of maximum simulated current speeds as a function of water depth of all
three scenarios.

Figure 2-11a shows the maximum current speeds from all three hypothetical cases predicted
by MOST. The model results suggest that high current speeds driven by eddies formed both within

29
and outside of the Pillar Point harbor near the entrance, as a result of the strong advection of the
flow. At the entrance to the inner harbor, however, the maxima plot shows relatively slower current
speeds. Although the maximum current plots give an idea about the spatial variability of the
current patterns, it is not easy to perceive localized scenario-based differences. To refine the local
variations in the maximum speeds from different sources, Figure 2-11b and 2-11c show the
absolute and relative differences in the Alaska and Chile scenarios with respect to the Mariana
case.

Figure 2-10: Correlation of the maximum currents speeds of Mariana scenario versus maximum
current speeds from a) Alaska scenario and b) Chile scenario predicted at every grid point in the
innermost grid. The solid line indicates the 1:1 relation and the two dashed lines indicate 1 m/s
error bounds.

The largest deviations are observed outside of the harbor, which can be attributed to energetic
eddies following different paths. Nevertheless, the error patterns inside both harbors show a great
similarity. In the vicinity of the entrance of the outer harbor, the relative error is around 20% for
both the Alaska and Chile scenarios where the maximum speeds occur, and 10%–15% in the inner

30
harbor. This analysis clearly shows that for the harbor and sources examined, the effect of offshore
directionality and tsunami frequency content has a very weak effect on the maximum currents
experienced in the harbor.

Figure 2-11: a) Maximum current speeds estimated by the model in Pillar Point Harbor, b) absolute
difference in maximum current speeds from Alaska and Chile cases with respect to Mariana case,
and c) relative (%) difference in maximum current speeds from Alaska and Chile cases with respect
to Mariana case.

31
The much more important dependency on maximum currents is on the near-harbor amplitude
of the wave, indicating that currents in a harbor from a tsunami generated by a large far-field
earthquake may be reasonably predicted with only information about the predicted tsunami
amplitude. Secondary effects, such as the duration of the maximum current and the persistence of
currents in the harbor, are likely to have stronger source dependence.
2.5 Conclusions
Observations after recent transoceanic events in California have revealed that even though
tsunamis from distant sources do not cause significant inundation and overland flow, they can still
induce strong currents with erratic flow patterns, which can be damaging to infrastructure and
vessels in ports and harbors. In this study, a sensitivity analysis was conducted to understand the
effects of tides and wave directionality on localized tsunami-induced currents. To this end, it has
addressed the fundamental question of whether the inclusion of tides in tsunami simulations has
any significant impact on the tsunami-generated currents. It has been found that the tide-induced
variability in the tsunami currents ranged from 10% to 25%, which is a variability that is, currently,
lower than a typical forecast accuracy.
However, it should be kept in mind that these results do not provide a general rule in
quantifying the specific tide effects on tsunami currents for a generic location. For the same
tsunami event and tide signal, two neighboring harbors in San Diego Bay reacted dissimilarly.
Likewise, the pattern of tide effects for Pillar Point Harbor and Crescent City Harbor also did not
agree, which showed the influence of the physical characteristics of individual harbors, such as the
bathymetric features or existing infrastructure, as well as the strength of the tsunami and tide
signals, on these types of analyses. Therefore, to better accommodate the significance of tides in a
particular area of interest in relation to tsunami currents, it is strongly recommended that similar

32
analyses be conducted and assessed on a site-specific basis, which will give the most accurate
results. Finally, the findings regarding the role of the wave source region on the localized tsunami
currents were summarized, with the hope of determining whether decisions about maritime
advisories or warnings might be made based only on the local wave height predictions. It has been
found that statistically, the source location does not play a significant role in the prediction of
maximum speeds and that spatial variability in maximum currents is dominated by the presence of
eddies, for which deterministic simulation is of limited use.

33
3. Sediment Transport
The overall process of sediment transport can be summarized in three phases: erosion,
transport, and deposition. Each phase can be approached numerically by the assembly of an
integrated model compounded of three separate models; one model each for the hydrodynamics,
the sediment transport, and the morphodynamic evolution (Rakha et al.,1997; Cao et al., 2004;
Díaz et al., 2008; Kim & Lee, 2011). Representing the underlying physics and coupled within an
integrated model, the developed tool is capable of simulating sediment processes in dynamic
environments where a shallow flow field interacts with the sediments in the sea bottom, resulting
in water depth changes.
The exploration of sediment transport predictions using computer models increased rapidly
over the last 30 years with improved computer power. This is also echoed in research communities
where various sediment transport and morphodynamic models have been developed. Earlier
studies provide a thorough overview of these existing models, in terms of their accuracy,
applicability, and limitations (e.g. Merritt et al., 2003; Papanicolaou et al, 2008 & Sugawara et al.,
2014). One common problem that is emphasized in all these studies regarding the existing
sediment transport models is the high level of empiricism they involve since their formulations are
highly tied to experimental and field calibration. Therefore, it should be admitted that sediment
transport models incorporate a certain degree of simplification to be computationally feasible, and
at the present stage, no reliable and comprehensive theoretical formulas can describe the two-phase
phenomenon of sediment and flow (Papanicolaou et al. 2008).
Additionally, Merritt et al. (2003) state that simpler models can perform better in terms of
robustness, and hence provide more stable results than more complicated models. As the model

34
gets more sophisticated, or in other words, as the number of processes considered and related
parameters increases, models start to run the risk of imposing higher uncertainties regarding the
model inputs (Wheater et al., 1993). The advantage of using a more complicated model, which
provides a more realistic realization of the process, may well be offset due to these uncertainties.
Moreover, there are some other studies in the literature that argue that empirical models can be
more accurate than models with more complicated structures, when used within the developed
framework, (Ferro & Minacapilli, 1995; Letcher et al., 1999; Perrin et al., 2001). According to
Mettier et al, 2003, and Steefel, C.I. & Van Cappellan, 1998, the main consideration in determining
a model’s value is its simplicity relative to its explanatory power, as long as it can explain the
observations sufficiently.
In this proposed research, following the above discussion, erosion and deposition phases are
determined using empirical relations, whereas the transportation of entrained particles is estimated
using an advection-diffusion scalar transport equation, which is evaluated numerically.
A one-dimensional (1D) sediment transport model, similar to the one used in this present study,
was also utilized by Cao et al. (2004). They coupled a depth-averaged transport equation with a
shallow water model (SWE) and solved the system using a weighted averaged flux (WAF) method
together with Harten-Lax-Van Leer (HLLC) approximate Riemann solver to estimate the change
in a sediment bottom due to dam-break flows. Then, Wu & Wang, (2007) tackled the same 1D
problem with a first-order accurate FVM. Similarly, Simpson & Castelltort, (2006) used a first-
order accurate approximate Riemann solver but solved the problem in two-dimensions (2D). Also,
Xia et al., 2(010) considered the transport of bed load as well as the graded sediment effects in 2D
and utilized a second-order accurate monotone upstream-centered scheme. Later, Kim & Lee, 2011
included the effects of variable density and solved a two-dimensional version of a similar transport

35
equation using a fourth-order accurate FVM. All these studies obtained reasonable results in
comparison to laboratory experiments investigating dam break flows. However, neither tested their
model with large-scale problems. Additionally, these studies compute the source and sink terms
due to erosion and deposition of the sediment particles in the scalar transport equation by using
empirical functions.
3.1 Scalar Sediment Transport Model
The transport model used in this research is given in Eqn. 3-1, which is a depth-averaged
advection-diffusion equation. The exchange of sediment between the flow and the bottom is
determined by two individual mechanisms, i.e., entrainment of the bed sediment due to turbulence
and sediment deposition under gravitational action.

𝜕𝐻𝐶
𝜕𝑡
+
𝜕𝐻𝐶𝑈
𝜕𝑥
+
𝜕𝐻𝐶𝑉
𝜕𝑦
=
𝜕
𝜕𝑥
\𝐾
^
𝐻
𝜕𝐶
𝜕𝑥
_+
𝜕
𝜕𝑦
\𝐾
^
𝐻
𝜕𝐶
𝜕𝑦
_+𝐸−𝐷 3-1
Here, 𝐻 is the flow depth while 𝑈 and 𝑉 are the depth-averaged flow velocities. C symbolizes
the depth-averaged sediment concentration and 𝐾
^
is the sediment diffusion coefficient in the
horizontal plane. In this study, 𝐾
^
is assumed to be the same as the flow eddy viscosity given in
Smagorinsky’s turbulent eddy viscosity model (Rakha et al. 1997). Erosion (E) and deposition (D)
are additional source and sink terms accounting for the addition and removal of bed material to
and from the water column.

𝐾
^
= (𝐶
`
∆𝑥∆𝑦)
Q
b
2𝜎
de
𝜎
de

3-2
In Eqn. 3-2 𝐶
`
is a constant and set to 0.2 (Chen et al., 1999). ∆𝑥 and ∆𝑦 are the grid sizes in
the x and y directions respectively, and 𝜎
de
is the strain tensor in the horizontal plane. This approach

36
involves a non-capacity formulation, in which deposition and erosion processes are assumed as
independent empirical functions and modulate the sediment transport rate by transferring mass
across the bottom boundary of the flow. Finally, the change in the bottom morphology is evaluated
using the following first order partial differential equation;

𝜕𝑧
g
𝜕𝑡
=
∇.𝒒
1−𝑝

3-3
Where 𝑧
g
is the bottom elevation, 𝑝 is the porosity of the bottom, and 𝒒 is the net sediment
flux across the bottom boundary and is simply the difference between the erosion, 𝐸, and the
deposition, 𝐷, rates.
3.2 Erosion and Deposition Rates
Sediment exchanges between a sediment boundary layer and flow fields are estimated through
erosion and deposition fluxes. Various types of formulas have been developed for erosion and
deposition estimates under either steady or unsteady flows and are mostly based on experimental
data (see for instance Castro D´ıaz et al., 2008; Abderrezzak & Paquier, 2010). Among those, an
empirical formula proposed by Cao et al. (2004) is implemented in the present model as it has
shown accurate predictions over shallow flows in recent studies (Simpson and Castelltort; 2006,
Xia et al. 2010; Kim and Lee, 2011).
The calculation of the entrainment and deposition rates continues to be one of the pivotal
components of computational models for sediment transport and morphological evolution. The
following scaling relations for the determination of erosion and deposition fluxes are introduced
in Cao (1999).

37

𝐸 =
𝜑(𝜃−𝜃
k
)√𝑢
Q
+𝑣
Q
ℎ(𝑑
mn
)
Rn.Q
𝑖𝑓 𝜃 ≥𝜃
k

3-4

𝐷 = 𝛼𝐶𝑤
s
(1−𝛼𝐶)
t

3-5
In Equations 3-4 and 3-5, 𝐸 and 𝐷 are the erosion and deposition rates respectively. 𝜃 is the
Shields parameter and calculated using Eqn 3-6. 𝜃
k
is the critical Shields parameter set for the
initiation of the motion of the particles, and is assumed as 0.045. The porosity of the sediments in
the sea bottom are represented by 𝑝. A coefficient 𝛼 is used to determine the depth-averaged
concentration and is usually larger than unity, since the near-bed concentration, 𝐶, cannot be larger
than 1−𝑝 physically (Cao, 2004), and it is calculated by Eqn. 3-7 empirically. 𝑑
mn
is the median
grain size of the sediment particles, and assumed to be constant for the entire solution domain.

𝜃 =
𝑢
g
Q
(𝑆𝐺−1)𝑔𝑑
mn

3-6

𝛼 =min (2,
1−𝑝
𝐶
)
3-7
The fall velocity of the sediment particles, 𝑤
s
, in still water in Eqn. 3-5 can be computed by
solving Equations 3-8 to 3-10 iteratively (Ponce, 1989).

𝑤
n
=z
4
3
𝑔𝑑
mn
𝐶
|
𝜌
`
−𝜌
~
𝜌
`

/Q

3-8

38

𝐶
|
=
24
𝑅𝑒
+
2.6
𝑅𝑒
5.0

1+
𝑅𝑒
5.0

.mQ
+
0.411
𝑅𝑒
263,000

R.S
1+
𝑅𝑒
263,000

R.n
+
𝑅𝑒
n.
461,000

3-9

𝑅𝑒 =
𝑤
n
𝑑
mn
𝜈
3-10
3.3 Numerical Methods
In this section, numerical methods that are utilized to solve the equations to estimate the
transport of the sediment particles (Eqn. 3-1), and the evolution of the bottom morphology (Eqn.
3-3) will be discussed in detail.
For the time integration of the conserved variable in Eqn. 3-1, 𝐻𝐶, the third-order accurate
Adams-Bashforth predictor scheme and the fourth-order accurate Adams-Moulton corrector
scheme are utilized, whereas the spatial representation of the numerical grid is cell-averaged using
a Finite Volume scheme. At the predictor step, values at a new time step (𝑛+1) are evaluated
explicitly using the known values at time steps, (𝑛−2), (𝑛−1) and 𝑛 by;

(𝐻𝐶)

= (𝐻𝐶)

+
Δ𝑡
12
(23𝐹

−16𝐹
R
+5𝐹
RQ
) 3-11
Then the values calculated at all four time steps including (𝑛+1) are input to the implicit
corrector step given in Eqn. 3-12, to improve the prediction of the variables in the step (𝑛+1).

(𝐻𝐶)

= (𝐻𝐶)

+
Δ𝑡
24
(9𝐹

+19𝐹

−5𝐹
R
+𝐹
RQ
) 3-12
Where; 𝐹 is given by;

39

𝐹 = −\
𝜕𝐻𝐶𝑈
𝜕𝑥
+
𝜕𝐻𝐶𝑉
𝜕𝑦
_+
𝜕
𝜕𝑥
\𝐾
^
𝐻
𝜕𝐶
𝜕𝑥
_+
𝜕
𝜕𝑦
\𝐾
^
𝐻
𝜕𝐶
𝜕𝑦
_+𝐸−𝐷 3-13
This corrector iteration is repeated until the calculated values of the variables converge within a
given tolerance.
A finite volume method using the fourth-order compact MUSCL-TVD (Monotone Upstream
centered Scheme for Conservation Law – Total Variance Diminishing) scheme is implemented to
solve the leading-order conservative terms of the transport equation (Yamamoto & Daiguji, 1993;
Kim et al, 2008;). This method shows a great capability of capturing shock fronts without
numerical dispersion and is proven to be highly stable and accurate (Toro, 2001). The first step of
this method is to construct a Riemann problem by determining the interface values of flux terms
around each cell. Once conserved variables are reconstructed at each interface of piecewise
volume; final values are determined using the HLL Riemann solver introduced by Toro (2002).
During this process, reconstructed variables are usually limited by slope limiter functions in order
to avoid under- or over-shooting. The cell interface values at 𝑖±1/2 are calculated as follows:

(𝐻𝐶)
d/Q

= (𝐻𝐶)
d
+
1
6
Δ
∗
(𝐻𝐶)
dR/Q
+2Δ
∗
(𝐻𝐶

)
d/Q

3-14

(𝐻𝐶)
d/Q

= (𝐻𝐶)
d
−
1
6
2Δ
∗
(𝐻𝐶)
d/Q
+Δ
∗
(𝐻𝐶

)
dT/Q
3-15
The subscript i refer to a computational cell and superscripts L and R indicate the left and the
right states of the variable HC at the i
th
cell respectively. Also, Δ
∗
(𝐻𝐶)
dR/Q
, 2Δ
∗
(𝐻𝐶

)
d/Q
,
2Δ
∗
(𝐻𝐶)
d/Q
, and Δ
∗
(𝐻𝐶

)
dT/Q
is determined by the minmod flux limiter as;

40

Δ
∗
𝐻𝐶
dR

Q
=𝑚𝑖𝑛𝑚𝑜𝑑Δ
∗
(𝐻𝐶)
dR

Q
,2Δ
∗
(𝐻𝐶)
d

Q

3-16

Δ
∗
𝐻𝐶

d

Q
= 𝑚𝑖𝑛𝑚𝑜𝑑Δ
∗
(𝐻𝐶)
d

Q
,Δ
∗
(𝐻𝐶)
dR

Q

3-17

Δ
∗
𝐻𝐶
d

Q
= 𝑚𝑖𝑛𝑚𝑜𝑑Δ
∗
(𝐻𝐶)
d

Q
,Δ
∗
(𝐻𝐶)
d
T
Q
3-18
Δ
∗
𝐻𝐶

d
T
Q
= 𝑚𝑖𝑛𝑚𝑜𝑑Δ
∗
(𝐻𝐶)
d
T
Q
,Δ
∗
(𝐻𝐶)
d

Q
3-19

Δ
∗
𝐻𝐶
d

Q
= Δ(𝐻𝐶)
d

Q
−
1
6
Δ
T
𝐻𝐶
d

Q
3-20

Δ
T
𝐻𝐶
d

Q
= Δ𝐻𝐶
dR

Q
−2Δ𝐻𝐶
d

Q
+Δ𝐻𝐶
d
T
Q

3-21
Where;
𝑚𝑖𝑛𝑚𝑜𝑑(𝑖,𝑗) = 𝑠𝑖𝑔𝑛(𝑖)max [0,min {|𝑖|,𝑠𝑖𝑔𝑛(𝑖)𝑗} 3-22

Δ𝐻𝐶
dR

Q
=𝑚𝑖𝑛𝑚𝑜𝑑Δ(𝐻𝐶)
dR

Q
,2Δ(𝐻𝐶)
d

Q
,2Δ(𝐻𝐶)
d
T
Q

3-23

Δ𝐻𝐶
d

Q
=𝑚𝑖𝑛𝑚𝑜𝑑Δ(𝐻𝐶)
d

Q
,2Δ(𝐻𝐶)
d
T
Q
,2Δ(𝐻𝐶)
dR

Q
3-24

Δ𝐻𝐶
d
T
Q
=𝑚𝑖𝑛𝑚𝑜𝑑Δ(𝐻𝐶)
d
T
Q
,2Δ(𝐻𝐶)
dR

Q
,2Δ(𝐻𝐶)
d

Q

3-25
Where;
𝑚𝑖𝑛𝑚𝑜𝑑(𝑖,𝑗,𝐾)= 𝑠𝑖𝑔𝑛(𝑖)max [0,min {|𝑖|,𝑠𝑖𝑔𝑛(𝑖)𝑗,𝑠𝑖𝑔𝑛(𝑖)𝑘} 3-26

41
Subsequently, the HC values are constructed at all interfaces, and the resulting formulation
becomes a Riemann problem, which can be solved by the HLL approximate Riemann solver
introduced which in Toro, (2002) and discussed fully in Kim et al., (2009).

𝐹
d/Q
=
𝐹

, 0 ≤𝑆

𝐹
∗
, 𝑆

≤0 ≤𝑆

𝐹

, 0 ≥𝑆

3-27

𝐹
∗
=
𝑆

𝐹

−𝑆

𝐹

+𝑆

𝑆

(𝑈

−𝑈

)
𝑆

−𝑆

3-28
The wave speed estimates, S
£
and S
¤
are obtained by;

𝑆

= 𝑈

+𝑎

𝑞

𝑆

= 𝑈

−𝑎

𝑞

3-29
In which, 𝑎

and 𝑎

are wave celerities, 𝑈

and 𝑈

are the depth-averaged velocities at the left
and right interfaces of the computation cell. The flux terms, q
,
, in Eqn. 3-29 are given by;

𝑞
,
⎩
⎨
⎧
ª
1
2
(𝐻
∗
+𝐻
,
)𝐻
∗
(𝐻
,
)
Q
, 𝐻
∗
> 𝐻
,
1, 𝐻
∗
≤ 𝐻
,
3-30

𝐻
∗
=
1
g
z
1
2
(𝑎

+𝑎

)+
1
4
(𝑈

+𝑈

) 3-31
Finally, in Eqn. 3-30, the case when 𝐻
∗
>𝐻
,
indicates a shock, while the opposite means a
rarefaction.

42
3.4 Validation of the Model
3.4.1 Laboratory Experiment on 2D Dam Break Flow
The developed model was first tested with a laboratory experiment, which analyzed a partial
dam-breach over a mobile bed, conducted in the Hydraulics Laboratory of Tsinghua University,
China (Xia et al. 2010). The experiment was performed in a channel that is 18.5 m long and 1.6
m wide, and the 4.5 m of the downstream bed was erodible and made up of coal ash (See Figure
3-1). The median particle size, 𝑑
mn
, of the coal ash was 0.135 mm with a specific gravity (SG) of
2.248. Initial water depths were set to 0.4 m in the upstream and 0.12 m in the downstream
respectively. A 20 cm wide breach formed instantaneously at t=0 at the center of the dam located
at X=2m, and the generated shock wave propagated towards to downstream, significantly eroding
the bottom. Bed levels at two cross sections at X=2.5 m, and X=3.5 m were measured 20 s after
the start of the experiment, using an ultra-acoustic topographic surveying meter. However, there
weren’t any measurements of water depth and velocity variations made due to the use of coal ash.

Figure 3-1: Experimental setup of 2D dam-break flow through a partial breach.

43
When re-creating the physical model in the computations, the grid size was selected as 0.01 m
and the time step set to 0.002 sec. Figure 3-2 shows the differences in bottom profile between the
simulated and measured data. As in the experiment, the water level in the reservoir and the
downstream are set to 0.4 m and 0.12 respectively. While the sediment transport model
underestimated the erosion in CS-1, it captures the changes in the flanks quite well. Whereas in
CS-2 it significantly overshot both the eroded section in the center, and predicted deposition in the
flanks of the channel. The main reason for the poor model performance in this test are that
numerical instabilities arose in the hydrodynamical module due to insufficient numerical
dissipation and/or lack of higher-order physics. For example, the role of turbulence becomes more
profound in these types of problems where horizontal and vertical circulations formed by the jet-
like flows occur, as discussed in Xia et.al., (2010).

Figure 3-2: Comparisons of bottom profile calculated by model to experimental data at a) Cross-
section 1, b) Cross-section 2.

44
3.4.2 Large Scale Tests in Crescent City and Santa Cruz Harbors
As discussed in detail in the introduction, the 2011 Tohoku Tsunami caused severe
sedimentation problems in some California harbors. Among them, Santa Cruz and Crescent City
harbors were investigated by CGS, and the outcomes of the survey were reported in Wilson et al.
(2012). These survey results provided useful data to test the present sediment transport model in
the large scale.
In this study, MOST was used to propagate tsunami waves from their source to a nearshore
region through nested grids. The model propagated the tsunami waves to the shore, computing the
wave amplitude, velocity, and overland inundation and captured the non-linearity of the waves as
they reached shallow water. Computations were stopped at the 5-m-depth contour in all the parent
grids, at which depth the waves were reflected back (a solid wall boundary condition was
imposed). Runup, inundation, and sediment transport computations were only performed in the
innermost grid, in which a bottom friction term was included in the momentum equations. For all
the simulations presented here, the Manning’s “n” friction factor was 0.03.
Similar to the analysis presented in the previous chapter, the outermost grid at the 4-arcmin
resolution covered the entire Pacific basin, and three additional grids of increasingly finer
resolution were derived from data obtained from NOAA’s database. The innermost nearshore grid
with the highest resolution was at the 1/3 arcsec grid, with boundary inputs for free surface
elevation and velocities from the previous MOST nested layer. Details of computational grids for
all locations are listed in Table 3-1.

45
Table 3-1: Nested grid system for each study site.
Lon. Range (
o
E) Lat Range (
o
N) Nx Ny dx dt (sec)
Propagation 120 – 292 -74 – 62 2581 2879 4’ (~7200 m) 10
Crescent City
Level 1 234.2592 – 236.3725 39.0350 – 41.9884 318 444 24” (~720 m) 2.5
Level 2 235.5617 – 236.1183 40.5242 – 41.9542 168 430 12” (~360 m) 1.5
Level 3 235.7142 – 235.9133 41.6284 – 41.8150 240 225 3” (~90 m) 1.5
Level 4 235.7821 – 235.8375 41.7166 – 41.7536 600 400 1/3” (~10 m) 0.25
Santa Cruz
Level 1 236.0308 – 238.5392 36.0200 – 38.9117 302 348 30” (~900 m) 2.5
Level 2 236.7758 – 238.4258 36.4333 – 38.5333 397 505 15” (~450 m) 1.5
Level 3 237.9367 – 238.2083 36.8233 – 37.0442 327 266 3” (~ 90 m) 0.75
Level 4 237.9893 – 238.0077 36.9407 – 36.9776 200 400 1/3” (~10m) 0.25

The empirical parameter (𝜑) in Eqn. 3-4 was set to 5x10
-5
as suggested in previous studies
(Cao et al. 2004; Simpson and Castelltort, 2006; Wu and Wang, 2007; Xia et al. 2010; Kim and
Lee, 2012). As stated earlier, one major shortcoming of sediment transport models is the use of a
single grain size (d50) for the entire domain. To test the model’s sensitivity against d50, four
simulations have been run with different grains sizes varying from 0.15 mm to 0.35 mm. These
values were selected based on the description of tsunami deposits given in Wilson et al. (2012), in
which they are classified as “silty sand”. Last but not least, the fall velocity (wo) of the sediment
particles was calculated using Equations 3-8 to 3-10.
The results of sediment transport modeling are given in Figure 3-5 , in which the bathymetric
changes with respect to different d50 are plotted for Crescent City and Santa Cruz harbors
respectively. Figure 3-3 immediately suggests that the results in Crescent City harbor should be
examined in two sections. In the first section covering the small boat basin (SBB) and its’ the
entrance, the numerical model captures the overall process well. The zonation of scour region and
sedimentation areas match the field survey data for all d50, however, the amount of the transported
material decreases for larger d50 values.

46

Figure 3-3 Comparison of bathymetric changes in Crescent City Harbor obtained by a) model with
d50 = 0.15 mm, b) model with d 50 = 0.20 mm, c) model with d 50 = 0.25 mm, d) model with d 50 =
0.35 mm to e) field data provided by CGS.
0 250 500 750 1,000 125
Meters
0 250 500 750 1,000 125
Meters
0 250 500 750 1,000 125
Meters
0 250 500 750 1,000 125
Meters
0 250 500 750 1,000 125
Meters
Legend
High : 1.5
Low : -4
1.5
Elevation
(meters)
1.0
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
-3.0
-3.5
Legend
High : 1.5
Low : -4
1.5
Elevation
(meters)
1.0
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
-3.0
-3.5
Legend
High : 1.5
Low : -4
1.5
Elevation
(meters)
1.0
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
-3.0
-3.5
a) b)
c) d)
e)

47
On the other hand, in the entrance of the main harbor, the numerical model over-predicts the scour
area compared to the field survey results, when d50 is 0.15 mm. Although, the numerical model
results start getting closer to field measurements in terms of the patterns of scour and deposition
with increasing values of d50, the amount of erosion near the western jetty of the main entrance is
underpredicted when d50 is greater or equal to 0.25 mm. One conclusion that can be drawn from
this analysis is the importance of using variable d50 in large domains such as Crescent City harbor
where the sediment characteristics are likely to be more variable.
Finally, Figure Figure 3-4 shows the scatter plot of the field measurement data vs model
predictions. The mean regression line between the data and the standard deviation bounds are
shown with yellow line in Figure 3-4a in addition to the point cloud.

Figure 3-4: a) Scatter of measured data vs. model results in Crescent City. Red lines are the ± 1
standard deviation from the 1:1 regression line shown with black. Yellow line indicates the mean
of the data calculated within 0.5-meter bins. The colorbar indicates the density of the point cloud;
b) Change in mean of the data calculated within 0.5 meter bins with the empirical parameter 𝜑.

48
Since the mean regression line lies above the 1:1 regression line for the positive values, it indicates
that the model yielded lower deposition predictions in an averaged sense. Contrarily, the mean
regression line remains below the 1:1 regression line for the negative values, implying an overshot
in erosion values in the model. A sensitivity analysis has been performed to check if the model
predictions can be improved by changing the values of the empirical parameter, 𝜑 between 4𝑒−5
and 8𝑒−5, and the mean regression lines calculated for each case are plotted in Figure 3-4b. As
can be seen from the figure, model predictions of erosion and deposition up to 1 meter improve
significantly and gets closer to 1:1 regression line in the mean sense with larger values of phi. On
the other hand, despite changing values of 𝜑, the deposition amounts higher than 1 meter are
drastically and consistently underpredicted by the model. Nonetheless, the sensitivity analysis
showed that the model can provide good estimates of the bathymetric changes in the geophysical
scale if calibrated properly to assign the optimum values to parameters like, 𝑑50 or empirical
parameter, 𝜑.
In Santa Cruz harbor, the survey area only covers the vicinity of the harbor entrance: the region
that is in between two jetties. Therefore, simulated results can be compared in a very limited area.
The model results in Santa Cruz are consistent in terms of the sediment transport configuration for
all grain sizes, and only the amount of the sediment transported changes. In comparison to the field
survey data, the model misses the sedimentation along the eastern jetty. The survey results show
up to 1.5 – 2 meters of deposition inside the eastern jetty, wherein the model estimates an erosion
of 1.5 meters, which also decreases as the size of the sediments increase. Some of the deposition
areas along the western jetty are predicted to be scoured by the model; here it must be noted that
Wilson et al. (2012) mentioned a mismatch between the NOAA and harbormaster surveys at these
locations. They mention in their paper that “Comparison of the NOAA post-tsunami versus

49

Figure 3-5: Comparison of bathymetric changes in the entrance of Santa Cruz Harbor obtained by
a) model with d50 = 0.15 mm, b) model with d50 = 0.20 mm, c) model with d50 = 0.25 mm, d)
model with d50 = 0.35 mm to e) field measurements made by CGS.
Legend
High : 1.5
Low : -4
2.5
Elevation
(meters)
2.0
1.5
1
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
Legend
High : 1.5
Low : -4
2.5
Elevation
(meters)
2.0
1.5
1
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
Legend
High : 1.5
Low : -4
2.5
Elevation
(meters)
2.0
1.5
1
0.5
0
-0.5
-1.0
-1.5
-2.0
-2.5
a) b)
c) d)
e)

50
Harbormaster pre-tsunami data shows similar scour patterns and is suggestive of some shoaling
along the eastern jetty. However, scour appears more pronounced in this comparison than it does
when comparing the pre- and post-tsunami Harbormaster surveys and no indication of shoaling
along the western jetty is observed.” This description of bathymetric changes measured by NOAA
demonstrates a better agreement with the model results. The discrepancy between two
measurements collected on the same date may imply the model performance may be higher than
anticipated in Santa Cruz. Also, the difference between pre-tsunami bathymetries used by the
model and made available to CGS were different, which is also a factor contributing to the
prediction errors of the model.
3.5 Conclusions
The objective of this chapter was to develop a numerical tool to predict changes in the
nearshore bathymetry due to sediment transport process taking place during a tsunami.
Hydrodynamic inputs for the sediment transport model were gathered from a classical NSW
solver, MOST, yet any 2HD hydrodynamic solver can be coupled with the sediment transport
model easily if desired.
This tool is first tested using data given in Xia et.al., (2010), where the results of an
experimental 2D dam-break problem are published. Overall, the model a reasonable but not
excellent realization of the case, mainly due to the physical assumptions of the hydrodynamic
model. Given the physics involved in an instantaneous dam breach problem, a leading order NSW
model lacking higher-order dissipation and diffusion terms, like MOST, may have errors.
Next, the model is tested at the geophysical scale in Crescent City and Santa Cruz harbors, and
the model predictions in depth changes within these two locations were compared to the field

51
measurement data collected after the 2011 Tohoku Tsunami and provided by Wilson et.al., (2012).
Unlike the realization of the 2D dam break problem, model results in these two locations was
found to be in better agreement with the field measurement data). Particularly in Crescent City,
the model results captured the overall process fairly well for a certain combination of the input
parameters.
Finally, two sensitivity analyses have been performed to evaluate the effect of the median grain
size of the sediments, 𝑑50, and empirical parameter, 𝜑, on the predicted bathymetric changes in
Crescent City Harbor after the 2011 Tohoku Tsunami. It has been found that while the model
provided better estimates with smaller 𝑑50 in the small boat basin, bigger 𝑑50 values matched
better in the outer harbor, which implies that in a relatively large harbor like Crescent City,
sediment size may vary considerably and using a single value throughout the computational
domain may over simplify the problem. Likewise, model performance was improved remarkably
as the value of the empirical parameter, 𝜑, increased. Consequently, the outcomes of these
sensitivity analyses underlined once again the significance of the initial estimation of parameters
on the produced results, which is a major issue that is inherent to all sediment transport modeling
studies.
Nonetheless, it can be concluded that, with proper calibration of the input parameters, the
model can provide reasonable estimates of expected morphological changes induced by tsunamis
at the geophysical scale, and results can be used to provide quantitative guidance to harbor
authorities, which was the ultimate goal in developing this tool.

52
4. Modeling the Transport of the Large Vessels with the Effect of
Tsunami Currents
This chapter focuses on the development of a transport model to estimate the motion of large
vessels in ports under the effect of the tsunami-induced currents. The assessment of seakeeping
and maneuverability of the large vessels has always been of high interest of researchers whether
they are in unrestricted (water depth >> draft of the ship) waters, or in ports and canals where they
are restricted in both depth and breadth. However, the recent incidents in Port Salalah, Oman and
some ports in Eastern India during the 2004 Indian Ocean Tsunami, and the reports of similar
issues from Japan and California ports and harbors during the 2011 Tohoku Tsunami showed that
the maneuvering capabilities of vessels can become very limited when they are captured by strong
currents, and within eddies in particular. The vessels that break free and drift-off in and around the
port during tsunamis generally exhibits random behaviors, and thus poses a great risk of hazard to
the port and its vicinity.
The motivation for this study comes from the desire to address the issues based on the
observations discussed above and to be able to provide quantitative guidance to authorities to help
them develop mitigation measures against these types of hazards. For this purpose, a large vessel
transport model is two-way coupled with a hydrodynamic numerical model allowing the
interaction between the flow and the vessel to be accounted for. This interaction was ignored in
most of the earlier published efforts (e.g. Tomita and Honda, 2010; Sakakibara et al., 2010;
Kobayashi et al., 2011) that study large vessel motion during tsunamis. It will be discussed in this
chapter that the presence of deep draft ships in shallow and/or restricted waters has significant
influence on the flow field that shouldn’t be neglected. The interaction between flow and the vessel

53
escalate since the draft of the vessel in shallow waters reduces to the effective flow depth
considerably. Such a decrease in depth can cause the current speed and consequently the current
force exerted on the vessel to increase. As in previous chapters, hydrodynamic inputs are
determined using MOST, and the gathered information is used to estimate the inertial damping
and viscous forces due to the movement of the ship based on a linear ship motion model.
4.1 Inclusion of Ships as Pressure Disturbance
There are several effective techniques exist to model the fluid structure interaction. Immersed
Boundary (IB) method, probably the most popular one, discussed in Peskin (2002) and since then
has been used extensively in studies dealing with fluid structure interaction. In the IB method, the
fluid is represented in Eulerian coordinates while structures are in a Lagrangian coordinate system
with moving grid points over the computation mesh. More recently, Kamrin et al., (2012)
established a new method called Reference Map Technique (RMT), in which both the fluid and
the solids are defined in the Eulerian reference frame. Reference map is an inverse motion function,
and it is employed to map the deformed surface of the structure into a set of Eulerian points.
However, both methods postulate the fluid structure interaction by enforcing a no-slip condition
over the boundary of the solids. Yet, this aspect of IB and RMT methods becomes a limitation in
2HD models when attempting to model the interaction between the fluid and the floating objects.
An object with a no-slip boundary condition in a 2HD model will block the flow that is supposed
to pass underneath the floating object, hence the essence of the problem will change. Due to this
fact, the fluid structure interaction in this study is attained by including the vessels as free surface
pressure disturbances.
Waves generated by moving pressure disturbances have been of interest to various studies to
analyze and model the different physical processes. Since the literature in this field is quite prolific,

54
some select examples of experimental, analytical and numerical efforts concentrating on the waves
generated by moving pressure disturbance are summarized below.
In one of the earlier studies, Greenspan (1956) developed analytical solutions for the edge
waves created by moving atmospheric pressure distributions, while, An et al. (2012) looked in to
the same problem numerically. They used linear shallow water equations forced with horizontal
pressure distribution terms to predict the edge waves generated by hurricanes moving parallel to
the shoreline. Also recently, Vennell (2010) dealt with the resonance of trapped ocean waves
created by moving atmospheric disturbances theoretically for different bottom topography
conditions.
Ertekin et al. (1986) presented a derivation of Green-Naghdi (G-N) equations with horizontal
pressure terms and solved these equations numerically to analyze the waves created by a free
surface pressure moving at a constant speed in a shallow channel with finite width and flat bottom.
Their aim was to investigate the characteristics of the waves generated by this moving pressure
disturbance for various Froude numbers and maximum pressure (𝑃
t°J
) values. Lee et al. (1989)
on the other hand, carried out a more comprehensive study and investigated the generation and
propagation of the solitary waves generated by pressure disturbances theoretically and numerically
by employing the Boussinesq-type model derived by Wu (1987) and KdV model presented in Lee
(1989). They also compared their theoretical and numerical findings with the data that they
gathered through series of controlled experiments. Liu and Wu (2004) modified a fully nonlinear
and weakly dispersive depth-integrated Boussinesq-type model to include the effects of moving
free surface pressure. Their focus was the waves generated by a pressure distribution moving in
trapezoidal and rectangular channels, where they also compared and discussed their findings with
existing results attained by the forced KdV, Green-Naghdi equations and the Boundary Integral

55
Equation (BIEM) formulation. Likewise, Bayraktar and Beji (2013) and David et al. (2017) uses
Boussinesq-type with horizontal pressure terms to exclusively model the ship-borne waves; while
the former is for ships moving over constant depth, the latter has capability to compute the waves
over irregular bathymetry.
Following the literature, the momentum equations in x and y directions in MOST can be
modified by adding free surface pressure terms defined in the horizontal plane to force these
equations.
𝑢
I
+𝑢𝑢
J
+𝑣𝑢
L
+𝑔ℎ
J
= 𝑔𝑑
J
−𝐷𝑢−
1
𝜌
𝜕𝑃
n
𝜕𝑥
4-1
𝑣
I
+𝑢𝑣
J
+𝑣𝑣
L
+𝑔ℎ
L
=𝑔𝑑
L
−𝐷𝑣−
1
𝜌
𝜕𝑃
n
𝜕𝑦
4-2
In Eqns. 4-1 and 4-2, 𝑃
n
(𝑥,𝑦) is the free surface pressure disturbance; then the initial
conditions will look like;
𝑢 =0, 𝑣 =0, 𝜂 = −
𝑃
n
𝜌𝑔

4.1.1 Waves Generated by Moving Pressure Disturbances
The wave train created by a travelling ship consists of a long leading bow and a trailing long
hull wave, and the series of diverging waves inbetween. Physically, these diverging waves are
short and dispersive, therefore the nonlinear shallow water equations that lack the relevant
dispersion terms cannot produce the free surface generated by the given ship completely. Hence,
only the leading bow wave in front and the stern wave behind the vessel can be produced
adequately by an NSW model like the one employed in this study. However, the main purpose of

56
this study, as mentioned earlier, is to be able to capture the variations induced in flow like the
acceleration/deceleration of flow or eddies generated due to flow separation from the ship hull etc.
rather than focusing on the waves generated by vessels. Nonetheless, to assess the model’s
performance in producing the leading bow wave, our model results are compared to numerical
results published in Ertekin et al., (1986).
As discussed briefly above, Ertekin derived a set of fully nonlinear, depth-integrated Green-
Naghdi equations to study the soliton generation by moving pressure disturbances. The
dimensionless pressure disturbance 𝑃
n
(𝑥,𝑦) that he took is rectangular with smooth corners as
given by Eqn. 4-3. It is first defined inside the quarter of a rectangle and then the pressure
distribution is completed into a full rectangle using the symmetry.
𝑃
n
(𝑥,𝑦) = 𝑝
t°J
𝑓(𝑥)∗𝑔(𝑦) 4-3
𝑓(𝑥)=
⎩
⎨
⎧1, 0 <𝑥 <
𝛼𝐿
2
cos
Q
z
𝜋(𝑥−𝛼𝐿/2)
(1−𝛼)𝐿
,
𝛼𝐿
2
< 𝑥 <
𝐿
2

𝑔(𝑦)=
⎩
⎨
⎧1, 0< 𝑦 <
𝛽𝐵
2
cos
Q
z
𝜋(𝑦−𝛽𝐵/2)
(1−𝛽)𝐵
,
𝛽𝐵
2
<𝑦 <
𝐵
2

The numerical tests were carried out for Froude Number, 𝐹𝑟, of 1.0 and 1.2 with 𝑝
t°J
=0.3.
Moreover, for 𝐹𝑟 =1.2, additional simulations were done using 𝑝
t°J
values of 0.1 and 0.15. In
all these tests water depth, ℎ
n
, was constant and set to unity, while the dimensions of the pressure
distribution and the channel were laid out using the dimensionless relations presented below as
they appear in Ertekin’s paper:

57
2𝑏
ℎ
n
=8 ,
𝐵
2𝑏
= 0.5 ,
𝐿
𝐵
=2 , 𝛼 =0.7 , 𝛽 =0.4
In Figure 4-1 and Figure 4-2, the snapshots of the free surfaces taken at different times during
the simulations from each model are shown, which provide qualitative comparisons of the results.

Figure 4-1: The snapshots of the free surface at various instants for 𝐹𝑟 =1.0 and 𝑝
t°J
=0.3. Left
panel is taken from Ertekin et al., (1986); right panel shows the results obtained by MOST.

As Figure 4-1 demonstrates, the initial leading and trailing wave profiles are analogous for
both models at the early stages of the simulations with 𝐹𝑟 =1.0 and 𝑝
t°J
=0.3. However, after
UT/h
o
= 10
UT/h
o
= 30
UT/h
o
= 50
UT/h
o
= 60

58
𝑈𝑇/ℎ
n
=50, the soliton formation starts in Ertekin’s model, and appears to detach from the
pressure disturbance at 𝑈𝑇/ℎ
n
=80. On the other hand, in MOST the leading wave turns into a
bore with a steep front, displaying a dispersive character owing to numerical dispersion, but it still
is attached to the pressure disturbance. The depression waves generated in both models on the
other hand, have similar patterns; they both follow the pressure disturbance and get longer in time.

Figure 4-2: The snapshots of the free surface at various instants for 𝐹𝑟 =1.2 and 𝑝
t°J
=0.1. Left
panel is taken from Ertekin et al., (1986); right panel shows the results obtained by MOST.

UT/h
o
= 10
UT/h
o
= 30
UT/h
o
= 50
UT/h
o
= 80

59
On the other hand, results illustrated in Figure 4-2 come from the case where 𝑝
t°J
=0.1,
which is relatively a weaker disturbance. Therefore, even though the pressure disturbance was
faster (𝐹𝑟 = 1.2) in this case no solitons were generated, thus qualitatively MOST does a better
job in creating a similar free surface as Ertekin et al.’s throughout the simulation. In both models,
the leading waves are similar and preserve their shapes during the simulation, but in MOST the
flanks are slightly more amplified. As in the simulation with 𝑝
t°J
=0.3, the depression waves
elongate in both models. Although these comparisons are predominantly visual, it can be said that
in general the agreement of the created free surface elevations between MOST and Ertekin’s model
is sufficiently good.
A quantitative evaluation of the present model performance with respect to Ertekin’s model is
given in Figure 4-3 andFigure 4-4, in which free surface elevations generated by MOST are
compared to Ertekin’s data for different initial conditions. From Figure 4-3a, it can be seen that
the leading crest produced by MOST provides good agreement with Ertekin’s results. Both the
wave amplitude and the phase are captured well. However, MOST results, with the first trough,
start to phase out and diverge from Ertekin’s data. This difference is particularly clear after the
second wave and can be attributed to the lack of dispersion and numerical breaking in MOST. Free
surface elevations are further compared at a later time for the same case in Figure 4-3b. Although
it can be argued that the crest elevations of the leading waves look similar, the wave profiles are
significantly different. In Ertekin’s data, it is apparent that the solitons are formed and separated
from the pressure disturbance, whereas in MOST, the leading crest gets longer and with a steep
breaking front. There is waviness in the depression region in both models, but the profiles do not
match. This difference is somewhat expected due to the high nonlinearity (𝑎/𝑙 ≈0.5) of the
leading wave in Ertekin’s model, which can be handled numerically by a fully nonlinear dispersive

60
model, whereas the leading solitary wave is transformed in to a breaking bore in the MOST
simulation. On the other hand, Lee et al., (1989) states that during their laboratory experiments
they also have observed breaking of the leading waves when the peak pressure (𝑝
½¾¿
) of the
disturbance was greater than 0.15, although the numerical results of Ertekin do not indicate that.

Figure 4-3: Comparison of free surface elevations calculated in MOST with Ertekin’s model at a)
65ℎ
s
and b) 130ℎ
s
ahead of the initial position of the pressure disturbance for 𝑝 =0.3 and 𝐹𝑟 =
1.0.

It must be also noted here that MOST doesn’t have an implicit wave breaking model; the numerical
dissipation inherently mimics the wave breaking evolution effectively at the breaking point (Titov
and Synolakis, 1998).
a)
b)

61
Last but not least, water levels from Ertekin for 𝐹𝑟 =1.2 and 𝑝
t°J
=0.15 are compared in
Figure 4-4 with MOST results for 𝐹𝑟 =1.05 and 𝐹𝑟 =1.2. Initially, the leading wave amplitudes
from MOST simulations for different Froude Numbers agree well with Ertekin data, but as the
simulation evolves, the wave amplitude and consequently the nonlinearity of the leading wave in
Ertekin’s model increases. While the leading wave generated by the disturbance moving with 𝐹𝑟 =
1.05 matches reasonable well for all times with Ertekin data, the wave in 𝐹𝑟 =1.2 breaks almost
immediately and therefore the amplitude of the leading wave in this MOST simulation remains
significantly lower than the Ertekin’s. Although there aren’t any data provided by Ertekin for the
case when 𝐹𝑟 =1.05, nonetheless an experiment for this case is conducted in MOST to see if it
can handle the waves generated in slightly supercritical conditions, since the mechanism of the
soliton generation in Ertekin’s model is similar for the pressure disturbances moving with 𝐹𝑟 <
1.2.
For typical port conditions, where the minimum water depth is around 15 meters, 𝐹𝑟 = 1.1
corresponds to a vessel velocity of ~13 m/s, which is very unlikely to be gained by a conventional
vessel due to tsunami currents. Therefore, from these comparisons, it can be affirmed with some
foundation that the MOST model is able to produce a realistic estimate of the leading wave height
and phase generated by the moving pressure disturbances sufficiently and is useful in this analysis
despite this limitation in Froude Number.

62

Figure 4-4: Time series comparison of free surface elevations at different instances during the
simulation for 𝐹𝑟 =1.05 & 1.2 in MOST, and 𝐹𝑟 =1.2 in Ertekin, for 𝑝 =0.1.

UT/h
o
= 10
UT/h
o
= 30
UT/h
o
= 50
UT/h
o
= 70

63
4.1.2 Benchmark Test in OpenFOAM
Of particular interest to this study, as stated earlier, is to develop a model that can sufficiently
resolve the flow around vessels. Interaction of the crossflow with the vessels results in acceleration
in the flow passing underneath, or the wakes generated behind the vessel. They are primarily driven
by the total drag acting on the vessel, which has 2 major components: i) friction drag due to viscous
resistance of the water, ii) form drag due to shape of the ship hull (Sorensen, 1973). The effect of
the form drag is attained by including pressure gradient terms (𝜕𝑝/𝜕𝑥,𝜕𝑝/𝜕𝑦) in the momentum
equations. Here it is further argued that the friction drag can also be encompassed by increasing
the shear stress locally at the grid points at which pressure gradients are non-zero.
To prove the above hypothesis, a benchmark test was set-up using an open source CFD
software, OpenFOAM, to assess the proposed model’s ability to generate the wake field behind a
floating object for an idealized and simplified case. As OpenFOAM retains numerous solvers that
can address a wide variety of CFD problems, IHFOAM was selected to be used for this benchmark
problem. IHFOAM is a two-phase, incompressible fluid solver based on Reynolds Averaged
Navier-Stokes (RANS) equations using the volume of fluid (VoF) method and offers realistic wave
generation and absorption boundary conditions (Higuera et al., 2014). The governing equations
that the IHFOAM solves are given by;
𝜕𝑢 À
d
𝜕𝑥
d
=0 4-4
𝜕𝜌𝑢 À
d
𝜕𝑥
d
+𝑢 À
e
𝜕𝜌𝑢 À
d
𝜕𝑥
d
=−𝑔
d
𝑥
e
𝜕𝜌
𝜕𝑥
e
+𝜇
ÂÃÃ
𝜕
𝜕𝑥
e

𝜕𝑢 À
d
𝜕𝑥
e
+
𝜕𝑝
𝜕𝑥
d
4-5

64
In Eqn. 4-5, 𝜇
ÂÃÃ
is the effective coefficient of viscosity, which combines the effects of
dynamic molecular viscosity and viscous diffusion and/or dissipation due to turbulence in the flow
(𝜇
ÂÃÃ
=𝜇+𝜌𝜈
IÄÅg
). The dynamic viscosity of water is, 𝜇 =1𝑒−3 𝑘𝑔/𝑚⋅𝑠, while the
additional viscosity due to turbulence, 𝜈
IÄÅg
, is estimated based on the selected turbulence model
in OpenFOAM. This second order diffusion/dissipation term is one of the key differences between
OpenFOAM and MOST when it comes to physics included in the underlying theory of the models.
In simulations where turbulence-controlled features such as large eddies and wakes are generated
in the flow, the significance of these terms gets large (Kim et al., 2009). Therefore, OpenFOAM
is expected to provide results closer to reality and it can set a benchmark to assess the performance
of MOST with the modifications made.

Figure 4-5: Schematic view of the benchmark set-up. a) Top-view, showing the block and the
numerical gauges where the time series are extracted; b) Side view of the set-up (vertically
exaggerated).

G1 G2 G3 G4 G5 G6 G7 G8
(G9-G13)
(G14-G18)
(")
($)
0.1 1 2 3 7 8 9 9.9 4 4.5 5 5.5 6
1
0.75
0
-0.75
-1
-0.5
0.5
1 2 3 7 8 9 4 5 6 0
0
10
5
1
2
3
4
y (m)
z (cm)
x (m)

65
The test was conducted in a 10-meter-long and 2-meter-wide shallow tank, where a floating
block is placed in the center of the channel (Figure 4-5). A steady upstream flow discharged into
the tank, where the water is initially at rest, with a depth averaged velocity (𝑈) of 0.1 𝑚/𝑠. The
water depth (ℎ
n
) in the tank was 0.05 m, which translates in to 𝑅𝑒 ≈ 5000. The floating block
was 2 meters in length and 1 meter wide, with a draft of 0.04 m. The draft to water depth ratio
(𝐷/ℎ
s
) of the block is 0.8. Since the block remains immersed during the entire simulation, the
vortices in MOST were physically generated through the gradients in the pressure and the bottom
stress.
The time evolution of the vortices generated by the present model and OpenFOAM in the lee
of the floating block are compared using different measures. The first plots shown in Figure 4-6,
are the snapshots of the velocity vector maps at different times during the simulations. According
to the figure, the flow starts to separate from the right and left tips of the block at around 𝑡 =10 𝑠.
The separation regions grow as the jets supply momentum through the gaps at the sides of the
block and two symmetric eddies with positive and negative vorticities are generated on the right
and the left corners of the block. MOST captures the generation of the eddies well at this initial
stage when compared to OpenFOAM. As the flow progresses, the eddies are carried away from
the block until they detach by 𝑡 =20 𝑠. They are further advected downstream until 𝑡~30 𝑠. Later,
the vortices in both MOST and OpenFOAM become stationary and their centers stay at around
𝑥 =6.5 𝑚 until the end of the simulations. Since the vortices are limited by the channel walls
along 𝑦 direction, the radius growth continues only along 𝑥 in time.
Next, the vorticity (𝜔) field maps and swirl strength (𝜆
kd
) of the vortices are compared in
Figure 4-7 and Figure 4-8 respectively. Swirl strength conveys the local frequency of the rotation

66
(𝑇 = 2𝜋/𝜆
kd
) and can be calculated for the vortices that experience swirling motion by using the
local velocity gradient tensor given in 4-6.
𝐃
Q|
=
⎣
⎢
⎢
⎢
⎡

𝜕𝑢
𝜕𝑥
𝜕𝑢
𝜕𝑦
𝜕𝑣
𝜕𝑥
𝜕𝑣
𝜕𝑦

⎦
⎥
⎥
⎥
⎤
4-6
The tensor given in Eqn. 4-6 has either two real eigenvalues or pair of complex conjugate
eigenvalues (Adrian et al., 2000). The swirl strength, 𝜆
kd
, is zero unless the vortex undergoes
swirling motion, otherwise 𝜆
kd
is equal to the imaginary part of the complex conjugate pair. The
difference between vorticity and the swirl strength is that the shear layers are excluded in the
former, so the boundaries defined by the 𝜆
kd
= 0 contains only the core of the vortex (Adrian et
al., 2000).
From Figure 4-7 and Figure 4-8, it can be seen that at early times, models show similar
behaviors. As vortices are advected downstream, after 𝑡 =20 𝑠, the energy becomes more
concentrated near the center in MOST, whereas the it gets more diffused in OpenFOAM. This
becomes more apparent as time evolves, and the vortices in OpenFOAM grow faster than their
counterparts in MOST. This fact is seen clearer in swirl strength maps; in MOST, the vortices and
the jets coming through the sides of the block gets more unstable in time than the ones in
OpenFOAM. Also, smaller eddies appear starting from 𝑡 ≈ 30 𝑠 in MOST, which don’t merge
with the main vortices and spread towards the downstream in an erratic fashion.

67

Figure 4-6: Velocity vector maps at various times during vortex generation and evolution. Left
panel MOST; right panel OpenFOAM results.

68

Figure 4-7: Vorticity (𝜔) maps at various times during vortex generation and evolution. The red
lines represent the 𝜔 =0 contour. Left panel MOST; right panel OpenFOAM results.

69

Figure 4-8: Swirl strength (𝜆
kd
) maps at various times during vortex generation and evolution. The
red lines represent the 𝜆
kd
=0 contour. Left panel MOST; right panel OpenFOAM results.

70
Conversely, in OpenFOAM the structure of the eddies is more coherent, while their energy
consistently decays due to viscous dissipation until the end of the simulation. This difference
observed between models can be explained by the lack of a viscous turbulence model in MOST
that limits the horizontal and vertical mixing, which also prevents cores to grow further in time.
However, despite these differences, given the limitations of MOST due the its fundamental theory,
these comparisons show its’ ability to provide reasonable realization of the case when flow passes
a floating object compared to a more sophisticated model in the larger extent.

Figure 4-9: Comparison of the time series of streamwise velocities, 𝑢 (𝑚/𝑠), estimated in MOST
(blue curves) and OpenFOAM (red curves), for gauges 1-8 (see Figure 4-5 for the exact locations
of the numerical gauges).
Gauge #1 Gauge #2
Gauge #3 Gauge #4
Gauge #5 Gauge #6
Gauge #7 Gauge #8

71

Figure 4-10: Comparison of the time series of streamwise velocities, 𝑢 (𝑚/𝑠), estimated in MOST
(blue curves) and OpenFOAM (red curves), for gauges 9-18 (see Figure 4-5 for the exact locations
of the numerical gauges).

The time series of the streamwise velocities, 𝑢 (𝑚/𝑠), estimated by models are also used to
validate the results in MOST. The numbering and the locations of the gauges are shown in Figure
4-5. In Figure 4-9, the velocities extracted at various locations in the front and the lee of the floating
block are compared. The comparison indicates that, in front of the block, the velocities from both
models agree with each other starting 𝑡~15 𝑠 until 𝑡~25 𝑠. After this time, the velocities stay
lower in OpenFOAM since more energy is reflected back off the face of the block than it does in
MOST. However, behind the breakwater, the velocities measured in OpenFOAM are ~30% lower
than the velocities extracted from MOST. The reason for this deficiency in velocities in
Gauge #9 Gauge #10
Gauge #11 Gauge #12
Gauge #13 Gauge #14
Gauge #15 Gauge #16
Gauge #17 Gauge #18

72
OpenFOAM can be due to the more aggressive energy absorbing in the outlet boundary in
OpenFOAM simulations. On the other hand, velocity predictions on both sides of the block, in
general, agrees well as shown in Figure 4-10. The velocities between 𝑡 = 10𝑠 and 𝑡 =20 𝑠 are
almost identical, which is the critical window in generation and separation of the vortices. Like lee
side velocity comparison, after 𝑡~20 𝑠, also here the velocities in OpenFOAM stays ~20% below
the MOST predictions.
4.2 Equations of Motion
The traditional approach is the use of two reference frames to define the equations of motion
(Figure 4-11). One of which is the fixed to earth frame 𝑂(𝑥,𝑦,𝑧) that may be coincided with the
origin of the computational grid and the other one is a body-fixed frame 𝑂
n
(𝑥
n
,𝑦
n
,𝑧
n
), which
moves and rotates with the ship. For surface vessels, the most commonly adopted position for the
body-fixed frame is on the center of gravity (CG) of the ship, which gives hull symmetry about
the x0-z0 plane and approximate symmetry about the y0-z0 plane, while the origin of the z0 axis is
relative to the still water surface (Price & Bishop, 1974). In the standard notation given in SNAME
(Society of Naval and Maritime Engineers) (1950) and ITTC (International towing tank
Conference) (1975), the variables describing the position and orientation of a ship are defined
relative to the inertial (earth-fixed) frame, and the coordinates are noted as [𝑥 𝑦 𝑧]
Ñ
and
[𝜙 𝜃 𝜓]
Ñ
. On the other hand, the linear [𝑢 𝑣 𝑤]
Ñ
and angular [𝑝 𝑞 𝑟]
Ñ
velocities, forces
[𝑋 𝑌 𝑍]
Ñ
and moments [𝐾 𝑀 𝑁]
Ñ
are expressed in body-fixed coordinate system.

73

Figure 4-11: Name of the translatory and the angular displacements of a ship, as presented in
SNAME (1950).

Here, we can define two vectors, one of which represents the position/orientation (Eqn. 4-7)
with respect to earth-fixed coordinate system, whilst the other expresses the linear/angular
velocities Eqn. 4-8 according to the body-fixed frame, and consequently, the rate of change of the
position/orientation vector can be obtained using the relation given in Eqn. 4-9.
𝑛 = [𝑥 𝑦 𝑧 𝜙 𝜃 𝜓]
Ñ
4-7
𝜈 = [𝑢 𝑣 𝑤 𝑝 𝑞 𝑟]
Ñ
4-8
𝑛̇ =𝑅(𝑛)𝜈 4-9
Where, 𝑅(𝑛) is a transformation matrix based on the Euler angles (𝜙,𝜃,𝜓) given as in
(Fossen, 1994);

74
𝑅(𝑛)= z
𝑅

(𝜙,𝜃,𝜓) 0
TJT
0
TJT
𝑅
Q
(𝜙,𝜃,𝜓)
4-10

𝑅

(𝜙,𝜃,𝜓)
= Ú
cos(𝜓)cos(𝜃) −sin(𝜓)cos(𝜙)+cos(𝜓)sin(𝜃)sin(𝜙) sin(𝜓)sin(𝜙)+cos(𝜓)cos(𝜙)sin(𝜃)
sin(𝜓)cos(𝜃) cos(𝜓)cos(𝜙)+sin(𝜙)sin(𝜃)sin(𝜓) −cos(𝜓)sin (𝜙)+sin(𝜃)sin(𝜓)cos(𝜙)
−sin (𝜃) cos(𝜃)sin (𝜙) cos(𝜃)cos(𝜙)
Û
4-11

𝑅
Q
(𝜙,𝜃,𝜓)= Ú
1 sin(𝜙)tan (𝜃) cos(𝜙)tan (𝜃)
0 cos(𝜙) −sin(𝜃)
0 sin(𝜙)/cos(𝜃) cos(𝜙)/cos(𝜃)
Û 4-12
Then, following Newton’s second law, the equations of motion of the vessel with respect to
the body-fixed coordinate system can be written in the vector form as in Eqn.4-13. In Eqn. 4-13,
𝑀 represents the mass and inertia matrix due to the rigid body dynamics, and 𝐹 is all the forces
and moments arose due to flow or the motion of the hull in water.
𝑀𝜈̇ =𝐹(𝜈̇,𝜈,𝑛) 4-13
Generally, in comparison to the other motions, pitch and heave are neglected for conventional
ships, which leaves us with 4-degrees of freedom (dof); surge, sway, yaw and roll (e.g. Perez &
Blanke, 2002; Skejic and Faltinsen, 2008). Additionally, as an initial assumption the motion of
ships is restricted in the horizontal plane, thus, the roll motion is also ignored in this study.
Subsequently, the equations of motions in 3-dof (surge, sway, and yaw) are expressed in Eqn. 4-
14 expanding the Newtonian approach given in Eqn. 4-13.
𝑚(𝑢̇ −𝑣𝑟) =𝑋
^Ä
+𝑋
~°ÝÂ

𝑚(𝑣̇ +𝑢𝑟)=𝑌
^Ä
+𝑌
~°ÝÂ

𝐼
ßß
𝑟̇ =𝑁
^Ä
+𝑁
~°ÝÂ

4-14

75
Where, 𝑚 is the mass, and 𝐼
ßß
is the moment of inertia of the vessel. The external forces acting on
the hull due to movement of the vessel and the flow passing through will be described in more
detail below.
The hull forces are the sum of all the hydrodynamic forces exerted on the hull due to its’ own
motion. They are also called the hydrodynamic derivatives and appear in the series expansion of
the hull forces;
𝑋
^Ä
≅ 𝑋
Ä̇𝑢̇ +𝑋
ÄÄ
𝑢
Q
+𝑋
Å|Ä|
𝑟|𝑢|+⋯ 4-15
Where,
𝑋
Ä̇ =
𝜕𝒇
ã
𝜕𝑢̇ , 𝑋
ÄÄ
=
𝜕
Q
𝒇
ã
𝜕𝑢
Q
, 𝑋
Å|Ä|
=
𝜕
Q
𝒇
ã
𝜕𝑟𝜕|𝑢|

The above coefficients represent a different component of the inertia forces. For instance, 𝑋
Ä̇
means the force in 𝑥
s
-direction due to the acceleration 𝑢̇ in 𝑥
s
-direction (Figure 4-11) and so forth.
But due to blind expansion of the series, higher order non-linear derivatives do also appear, and
the basic difficulty in this type of analysis is to determine the values of these nonlinear coefficients
for each hull shape and vessel type (Morgan, 1978). Following Newton’s second law which states
that the hull forces depend linearly on the acceleration only, we retain the first order acceleration
terms, and neglect the cross derivatives including both velocities and accelerations (Fossen, 1994).
Consequently, in this study, only the linear inertia terms are kept. The first of which are the added
mass and inertia terms arising due to the motion of the fluid associated with the accelerations
produced by the ship, which also reflects the build-up of the kinetic energy of the fluid as the hull
moves through it. Furthermore, the first order viscous damping terms are also important, especially

76
in the case of free motions, where the vessel drifts off (Bandyk & Beck, 2009). The resulting
equations of motion are:
𝑋
^Ä
= 𝑋
Ä̇𝑢̇ −𝑌
Ý̇𝑣𝑟+𝑋
ÝÅ
𝑣𝑟
𝑌
^Ä
= 𝑌
Ý̇𝑣+𝑋
Ä̇𝑢𝑟+𝑌
Ý
𝑣+𝑌
Å
𝑟
𝑁
^Ä
= 𝑁
Å̇𝑟̇ +𝑌
Å̇𝑣̇
4-16
In the above equation, the terms with coefficients 𝑋
Ä̇ , 𝑌
Ý̇ and 𝑁
Å̇ are the forces due to added mass
and inertia. The coefficients 𝑌
Ý
,𝑌
Å
and 𝑁
Ý
,𝑁
Å
that appear in sway force and yaw moment account
for the viscous damping. Some common methods in estimating these hydrodynamic parameters
are include CFD modeling or running controlled laboratory experiments and deriving the
coefficients by the analysis of the acquired data (Blanke & Jensen, 1997). If the Equations 4-14
and 4-16 are merged and arranged, they can be expressed in the matrix form as 𝑨∙𝑼
̇ =𝒃, where
the matrix and the vectors are;
è
𝑚+𝑋
Ä̇ 0 0
0 𝑚+𝑌
Ý̇ 𝑁
Ý̇ 0 𝑌
Å̇ 𝐼
ßß
+𝑁
Å̇ éè
𝑢̇ 𝑣̇ 𝑟̇ é =Ú
(𝑚+𝑌
Ý̇ +𝑋
ÝÅ
)𝑣𝑟+𝑋
~°ÝÂ
−(𝑚+𝑋
Ä̇ )𝑢𝑟+𝑌
~°ÝÂ
𝑁
~°ÝÂ
Û 4-17
External forces due to flow conditions are calculated using empirical relations. The governing
equations to determine the transverse and longitudinal forces arising from the tsunami currents on
vessels are acquired from aU.S. Army Corps of Engineers report (2005) on mooring line design
guidelines. The approach is intended to be first order, and thus differential loads are not treated in
this phase of the analysis. These equations have also been employed by Keen et al., (2017) and
found to be sufficient for this type of analysis. Hence, the static transverse drag force, 𝐹
Lk
acting
on the vessel is calculated using;

77
𝐹
Lk
= 0.5 𝜌
~
𝐿
~
𝑇 𝐶
Lk
𝑉
k
|𝑉
k
| sin (𝜃
k
) 4-18
Where, 𝐹
Lk
= Transverse current force (newtons); 𝜌
~
= mass density of water (𝑘𝑔/𝑚
T
), 𝑉
k
=
current velocity (m/s), 𝐿
~
= vessel waterline length (m); T = average vessel draft (m); 𝐶
Lk
=
transverse current force drag coefficient and 𝜃
k
= incident current angle (degrees). The drag
coefficient, 𝐶
Lk
which includes the effects of water depth and vessel’s underkeel is given by;
𝐶
Lk
=𝐶
n
+(𝐶

−𝐶
n
)\
𝑇
𝑑
_
Q

4-19
In Eqn. 4-19, 𝑑 is the water depth, 𝐶

is shallow water current force drag coefficient with the
recommended value of 3.2 (𝐶

= 3.2). The deep water drag coefficient, 𝐶
n
, can be obtained by;
𝐶
n
=0.22 √𝑥 4-20
𝑋 =𝐿
~
Q
𝐴
t
𝐵 𝑉
4-21
Where, 𝑋 = dimensionless ship parameter; 𝐴
t
= immersed cross-sectional area of the ship
midsection (𝑚
Q
); B = beam of the vessel; V = submerged volume of the ship (𝑚
T
).
Unlike transverse loads, the order of magnitude of the skin friction competes with the
longitudinal form drag for streamlined vessels. Therefore, both types of drag forces have to be
considered while calculating the total static longitudinal current forces. Neglecting the propeller
forces, the longitudinal drag is given by;
𝐹
Jk
=𝐹
Jëìí
+𝐹
Jëîï
4-22
Where,

78
𝐹
Jëìí
=0.5 𝜌
~
𝐵 𝑇 𝐶
Jkg
𝑉
k
|𝑉
k
| 𝑐𝑜𝑠𝜃
k
4-23
𝐹
Jëîï
=0.5 𝜌
~
𝐵 𝑆 𝐶
Jk°
𝑉
k
|𝑉
k
| 𝑐𝑜𝑠𝜃
k
4-24
In Eqns. 4-23 and 4-24; 𝐶
Jkg
is longitudinal current form drag coefficient and equals to 1 (𝐶
Jkg
=
1); and 𝐶
Jk°
and 𝑆 are the longitudinal skin friction coefficient and wetted surface area (𝑚
Q
)
respectively. Here these parameters are estimated using;
𝑆 = 1.7 𝑇 𝐿
~
+
𝐷
𝑇𝛾
~
4-25
𝐶
Jk°
=0.075/[(log
n
𝑅
ò
)−2 ]
Q
4-26
Where, 𝛾
~
is specific weight of water and 𝑅
ò
is the Reynolds number, and for vessels it is
computed using Eqn. 4-27;
𝑅
ò
=ó
𝑉
k
𝐿
~
cos (𝜃
k
)
𝜈
ó 4-27
Finally, the yaw moment 𝑁
~°ÝÂ
, can be calculated as the sum of the torque created due to
eccentricity (𝑥
ô
,𝑦′) of the lateral and longitudinal drag forces:
𝑁
~°ÝÂ
= 𝐹
Jk
𝑦′+𝐹
Lk
𝑥′ 4-28
4.3 Collision Model
The two main problems to address in a collision model are 1) the detection of the collisions
and calculating the velocity changes to break the colliding contact and 2) computing the necessary
contact forces that will lead to changes in velocities to avoid interpenetration. The details of the

79
collision detection algorithm and the collision solver that tackle these problems in a physically
correct manner will be explained in this chapter.
4.3.1 Collision Detection
The literature in collision detection is quite rich due to significant contributions of computer
graphics and robotics researchers. One of the generally accepted and applied methods is called the
Separation Axis Theorem (Baraff, 1989; Gottschalk et. al., 1996). This technique is generic for
rigid bodies defined as convex polyhedral, and particularly useful if collisions between objects
with arbitrary shapes are expected in a given model. However, in this study, particularly in large
scale simulations, the vessels are approximated as ellipses for certain practical reasons that are
described in section 4.4. This approximation in vessel shapes allowed the development of a simpler
and faster but non-generic algorithm for identifying collisions.
The collision detection algorithm starts with keeping track of the minimum distance between
the objects at every time step during the simulation. The minimum distance between two ellipses
consists of 𝑛 points such as 𝐸

= {𝑘

,𝑘
Q
,𝑘
T
,…,𝑘

} and 𝐸
Q
= {𝑙

,𝑙
Q
,𝑙
T
,…,𝑙

}, where each point
is defined in cartesian coordinates as 𝑘
d
(𝑥
d
,𝑦
d
) and 𝑙
e
(𝑥
e
,𝑦
e
), and can be defined as;
𝑑
td
(𝐸

,𝐸
Q
)=min{𝑑𝑘
d
,𝑙
e
}, 𝑖,𝑗 =1,2,…,𝑛
Where 𝑑𝑘
d
,𝑙
e
is the euclidean distance between 𝑘
d
and 𝑙
e
. Also, it must be noted that the
cardinality between two ellipses is not a requirement, yet it is assumed here such that to simplify
the notation. The collision detection algorithm works in the following steps;
- The minimum distances between all objects are tracked at all time steps.
- If the minimum distance between any two bodies (𝑑
td
(𝐸

,𝐸
Q
)) is less than a given
tolerance, (𝜖), then it means that the bodies are in colliding contact.

80
- If bodies interpenetrate at time 𝑡
n
+Δ𝑡 (Figure 4-12b), then the time of collision (𝑡
k
) is
approximated by evaluating 𝑣
ÅÂ
(𝑡
n
)/[𝑑
t d
(𝑡
n
)−𝜖]. Thus, the time of collision (𝑡
k
)
should satisfy 𝑡
s
< 𝑡
k
≤𝑡
s
+Δ𝑡.
- Once all the collisions are detected, they are solved in the order of occurrence and the post
collision velocities are calculated as the details will be discussed in the next section.

Figure 4-12: a) Ellipses 𝐸

and 𝐸
Q
, and the minimum distance, 𝑑
td
between them at time 𝑡
n
; b)
Two bodies collide and interpenetration occurs at time 𝑡
n
+Δ𝑡; c) The position of the bodies at
exact time of collision, 𝑡
k
, within given distance tolerance of 𝜖.

4.3.2 Collision Solver
Once all the collisions are detected at a given time step and the contact point on each body is
identified, the resulting impact force and the torque as well as the velocities of the bodies after the
collision is determined. One common approach for calculating the impact forces is the so-called
“penalty” method and is based on the introduction of the arbitrary forces (e.g. spring forces) to
separate the colliding bodies (Terzopolous et al., 1987). Although this method can generate
realistic animations, it is mostly preferred when the accuracy of the physics is less important. The
E
1
E
2
!
"#$
E
1
E
2
E
1
E
2
%
(') ()) (*)

81
same problem can be tackled analytically using the conservation of momentum principle as will
be discussed here. This method has been employed in computer graphics problems widely and
proven to provide physically correct contact forces (Hahn, 1988; Moore and Wilhelms, 1988;
Baraff, 1989).

Figure 4-13: a) Body A and body B approaching each other with prospective contact points 𝑝
°
(𝑡)
and 𝑝
°
(𝑡); b) Two bodies are in contact at 𝑡 =𝑡
n
, 𝑝
°
(𝑡
n
) = 𝑝
g
(𝑡
n
).

If the two bodies, body A and body B, are in contact at point 𝑝, then a given point on body A
denoted by 𝑝
°
(𝑡) should satisfy 𝑝
°
(𝑡
n
) =𝑝 at the time of the contact, 𝑡
n
. Similarly, point 𝑝
g
(𝑡)
on body B will coincide 𝑝
g
(𝑡
n
)= 𝑝 at time 𝑡
n
(Figure 4-13). While, at time 𝑡
n
, 𝑝
°
(𝑡) and 𝑝
g
(𝑡)
will overlap, the velocities of these two points might be very different. By evaluating the velocities
of these points, which will be denoted by 𝑝
°
̇ (𝑡
n
) and 𝑝
g
̇ (𝑡
n
) respectively, we can determine
whether these two bodies are approaching, separating or colliding each other. The point velocities
are given as:
𝑝
°
̇ (𝑡
n
) =𝑣
°
(𝑡
n
)+𝜔
°
(𝑡
n
) 𝕩 (𝑝
°
(𝑡
n
)−𝑥
°
(𝑡
n
)) 4-29
B
A
!
"
($)
!
&
($)
B
A
!
"
($
'
)
!
&
($
'
)
(") (&)

82
𝑝
g
̇ (𝑡
n
) =𝑣
g
(𝑡
n
)+𝜔
g
(𝑡
n
) 𝕩 (𝑝
g
(𝑡
n
)−𝑥
g
(𝑡
n
)) 4-30
Where, 𝑣
°
,𝑣
g
,𝜔
°
and 𝜔
g
are the linear and angular velocities and 𝑥
°
(𝑡
n
) and 𝑥
g
(𝑡
n
) are the
coordinates of the center of masses of the bodies at the time of the contact. Subsequently, the
relative velocity, 𝑣
ÅÂ
, of the contact points with respect to each other in the direction of contact
can be defined by Eqn. 4-31.
𝑣
ÅÂ
=𝑝
°
̇ (𝑡
n
)−𝑝
g
̇ (𝑡
n
)⋅𝑛 ø⃗(𝑡
n
)
4-31
In the above equation, 𝑛 ø⃗(𝑡
n
) denotes the normal surface vector, which provides the direction of
𝑣
ÅÂ
at the contact point. If 𝑣
ÅÂ
is positive, it means that 𝑣
ÅÂ
and 𝑛 ø⃗(𝑡
n
) are in the same direction
and the bodies are moving away from each other. Then the contact point will diminish at time 𝑡
n
+
Δ𝑡. Yet, we are concerned about the case when 𝑣
ÅÂ
<0, which means that the directions of the
𝑣
ÅÂ
and 𝑛 ø⃗(𝑡
n
) are opposite, and the two bodies are in contact at point 𝑝. To avoid an
interpenetration, the velocities of body A and body B supposed to change instantaneously. To
create such a discontinuity in body velocities, a contact force of 𝐹 should be exerted on the bodies
for a very short instant of time, Δ𝑡, which leads to a new quantity called the impulse ( 𝐽):
𝐽 = 𝐹Δ𝑡 4-32
The impulse is also a vector quantity with the units of mass times acceleration. Applying an
impulse will generate the required instantaneous adjustment in the velocity of a body. Then the
effect of 𝐽 on the linear velocity (𝑣) of a given rigid body with the mass, 𝑀, can be driven simply
by using the relation Δ𝑣 = 𝐽/𝑀. Additionally, the impulse acting at contact point 𝑝 also creates an
impulsive torque of 𝜏
ü
computed using Eqn. 4-33:

83
𝜏
ü
=𝑝
k
−𝑥
n
(𝑡
n
) 𝕩 𝐽
4-33
Likewise, the change in the angular velocity (Δ𝜔) of the body with the yaw moment of inertia
of 𝐼
ßß
due to the impulsive torque is given by Δ𝜔 =𝜏
ü
/𝐼
ßß
. Therefore if a collision takes place
between two bodies, an impulse and a resulting impulse torque will be applied to prevent them
interpenetrating.

Figure 4-14: a) The relative velocity vector at the time of the contact; b) The impulse due to
collision acting on each body; c) Post-impulse relative velocity vector.

Since the tangential or friction forces are neglected here, then the direction of the impulse will be
in the normal direction, 𝑛 ø⃗(𝑡
n
) (Figure 4-14a). Thus, the impulse 𝐽 can be written as;
𝐽 =𝑗 𝑛 ø⃗(𝑡
n
) 4-34
Where 𝑗 is the scalar to be calculated and gives the magnitude of the impulse. The sign convention
adopted here lets the impulse 𝐽 acts positively on body A, that is, A is subject to an impulse of
(+ 𝑗 𝑛 ø⃗(𝑡
n
)), while an equal but opposite impulse (− 𝑗 𝑛 ø⃗(𝑡
n
)) is exerted on body B (Figure 4-14b).
To derive the 𝑗 by using an empirical law for collisions, the following relations need to be identified
first:
A
B
! ($
%
)
'
(
̇ *
$
%
−'
,
̇ *
$
%
A
B
! ($
%
)
'
(
̇ *
$
%
−'
,
̇ *
$
%
A
B
!" (%
&
)
−!" (%
&
)
(𝒂) (𝒃) (𝒄)

84
𝑝
°
Q ̇ (𝑡
n
) = 𝑣
°
Q
(𝑡
n
)+ 𝜔
°
Q
(𝑡
n
) 𝕩 𝑟
°

4-35
Along with,
𝑣
°
Q
(𝑡
n
) = 𝑣
°

(𝑡
n
)+
𝑗𝑛 ø⃗(𝑡
n
)
𝑀
°

4-36
𝜔
°
Q
(𝑡
n
) =𝜔
°

(𝑡
n
)+
𝑟
°
𝕩 𝑗𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°

4-37
Where 𝑝
°
̇
(𝑡
n
) and 𝑝
°
̇ Q
(𝑡
n
) denote the velocity of the body A before and after the impulse 𝐽 is
applied, respectively, and 𝑟
°
=𝑝
k
−𝑥
°
(𝑡
n
) is the displacement distance. Also, the pre-collision
velocities of body A are denoted with 𝑣
°

(𝑡
n
) and 𝜔
°

(𝑡
n
), then the post-collision velocities are
shown by 𝑣
°
Q
(𝑡
n
) and 𝜔
°
Q
(𝑡
n
). Combining Eqns. 4-35 to 4-37 yields the following;
𝑝
°
̇ Q
(𝑡
n
) =𝑣
°

(𝑡
n
)+
𝑗𝑛 ø⃗(𝑡
n
)
𝑀
°
+𝜔
°

(𝑡
n
)+
𝑟
°
𝕩 𝑗𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°
𝑥 𝑟
°

4-38
𝑝
°
̇ Q
(𝑡
n
)= 𝑣
°

(𝑡
n
)+ 𝜔
°

(𝑡
n
) 𝑥 𝑟
°
+
𝑗𝑛 ø⃗(𝑡
n
)
𝑀
°
+
𝑟
°
𝕩 𝑗𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°
𝑥 𝑟
°

4-39
If we substitute Eqns. 4-29 and 4-30 into the above expression, we have:
𝑝
°
̇ Q
(𝑡
n
) = 𝑝
°
̇
(𝑡
n
)+𝑗
𝑛 ø⃗(𝑡
n
)
𝑀
°
+
𝑟
°
𝑥 𝑗𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°
𝕩 𝑟
°

4-40
𝑝
g
̇ Q
(𝑡
n
) = 𝑝
g
̇
(𝑡
n
)− 𝑗
𝑛 ø⃗(𝑡
n
)
𝑀
g
+
𝑟
g
𝑥 𝑗𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
g
𝕩 𝑟
g

4-41

85
𝑝
°
̇ Q
(𝑡
n
)− 𝑝
g
̇ Q
(𝑡
n
)= 𝑝
°
̇
(𝑡
n
)− 𝑝
g
̇
(𝑡
n
)+ 𝑗
ø⃗(I
ÿ
)
í
!
+
ø⃗(I
ÿ
)
í
"
+
+
Å
!
J e ø⃗(I
ÿ
)
(î
##
)
!
𝕩 𝑟
°
+
Å
"
J e ø⃗(I
ÿ
)
(î
##
)
"
𝕩 𝑟
g

4-42
Following Eqn. 4-31, we can also introduce 𝑣
ÅÂ

and 𝑣
ÅÂ
Q
, the pre and post-impact relative
velocities in the normal direction;
𝑣
ÅÂ

=𝑝
°
̇
(𝑡
n
)− 𝑝
g
̇
(𝑡
n
)⋅ 𝑛 ø⃗(𝑡
n
)
4-43
𝑣
ÅÂ
Q
= 𝑝
°
̇ Q
(𝑡
n
)− 𝑝
g
̇ Q
(𝑡
n
)⋅ 𝑛 ø⃗(𝑡
n
)
4-44
𝑣
ÅÂ

= −𝜖𝑣
ÅÂ
Q
4-45
Where 𝜖 is called the coefficient of restitution and ranges from 0 to 1. The upper limit, 𝜖 =1,
yields 𝑣
ÅÂ

= −𝑣
ÅÂ
Q
which means the entire kinetic energy of the system has been preserved and
the collision is perfectly elastic. If 𝜖 =0, then 𝑣
ÅÂ
Q
= 0, indicating a fully inelastic collision in
which bodies stuck each other after the collision.
If we rewrite the Eqn. 4-42 in terms of 𝑣
ÅÂ

, 𝑣
ÅÂ
Q
and 𝑗, and rearrange, then the expression to
compute the 𝑗 becomes;
𝑗 = −
−(1+𝜖)𝑣
ÅÂ

\
1
𝑀
°
+
1
𝑀
g
+𝑛 ø⃗(𝑡
n
) ⋅ \
𝑟
°
𝑥 𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°
_ 𝑥 𝑟
°
+𝑛 ø⃗(𝑡
n
)⋅\
𝑟
g
𝑥 𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
g
_𝑥 𝑟
g
_

4-46
In which 𝑀
°
,𝑀
g
are the masses, and (𝐼
ßß
)
°
and (𝐼
ßß
)
g
are the moment of inertias for the yaw
moment of the two colliding bodies. In the case where the body collides in to a non-movable body
like a berth, wharf or another port infrastructure, the terms with subscript 𝑏 will vanish in Eqn. 4-

86
46 since the mass and the inertia of an object that cannot be moved is infinitely large relative to
the colliding body. Therefore, for this type of collision, Eqn. 4-46 simplifies in to Eqn. 4-47:
𝑗 = −
−(1+𝜖)𝑣
ÅÂ

\
1
𝑀
°
+𝑛 ø⃗(𝑡
n
)⋅\
𝑟
°
𝑥 𝑛 ø⃗(𝑡
n
)
(𝐼
ßß
)
°
_ 𝑥 𝑟
°
_

4-47
Once the impact force, 𝑗, is determined, the post-collision velocities, 𝑢,𝑣 and 𝜔, for each body
involved in the collision can be calculated using Eqn. 4-36 and 4-37.
To validate the ability of the collision model, various benchmark experiments were produced
and example snapshots of the four selected cases are shown in Figure 4-15 to Figure 4-18. It was
assumed that the collisions were perfectly elastic in all tests. The first case, given in Figure 4-15,
was an eccentric collision that included two bodies with equal masses approaching to each other
with velocities of equal magnitude of 0.1 𝑚/𝑠.

Figure 4-15: Time lapse of the eccentric collision takes place between two bodies with equal
masses.
U
-U

87

Figure 4-16: Time lapse of the collision takes place between two bodies with different masses.

After the collision they started to spin in the same direction and moved apart from each other in
both 𝑥 and 𝑦 directions.
The last two collision tests involved three bodies of equal masses, where two of them were
initially at rest while the third one approached them with a velocity of −0.1 𝑚/𝑠. In the case shown
in Figure 4-17, the center of gravities of the bodies were aligned, and the two resting bodies were
in contact with each other. Once the moving body hit the middle one, the energy was transferred
o the outer body and it separated from the other two with a velocity of −0.1 𝑚/𝑠.
U
-U

88

Figure 4-17: Time-lapse of the three-body collision with equal masses, and their center of gravities
are on the same line. The body on the left moves with the initial velocity of −0.1 𝑚/𝑠, while the
other two are initially at rest.

Figure 4-18: Time-lapse of a three-body collision with equal masses. The body on the left moves
with the initial velocity of −0.1 𝑚/𝑠, while the other two are initially at rest.
-U
-U

89
In the final experiment presented in Figure 4-18, three bodies placed arbitrarily, whilst only
the leftmost object had an initial velocity of −0.1 𝑚/𝑠. Multiple collisions take place between
these bodies and their calculated positions at different instants are demonstrated in the figure. All
these tests provided realistic results and verified the validity of the collision model.
4.4 Large Scale Numerical Experiments
The developed model was tested in two locations for two different tsunamis. In the first test,
the effects of a hypothetical tsunami arising from Alaska-Aleutians Subduction Zone was
investigated for the ships moored in Port of Long Beach. The purpose of this analysis was to see
if the presence of vessels within the port changes the flow field created by the tsunami, and if so,
does the variation intensify as the number of the ships in and around the terminal increases? Then,
another experiment was conducted based on the post-tsunami data collected after the 2011 Tohoku
Tsunami in Port of Ishinomaki, Japan. Based on a personal communication with Japanese
researches in April 2018, it has been found out that a container ship broke free off its’ moorings
and was pulled out of the terminal by the currents during the event. A realization of this event was
made using the developed transport model based on the best estimates of the vessel’s draft, location
and release time at the time of the incident. However, the values of these parameters cannot be
known precisely before the incident. Therefore, a sensitivity analysis was also performed based on
this case to see the effects of the initial choice of these parameters on the predicted behavior of the
vessel.
In these experiments, the ships were approximated as ellipses, and the pressure disturbances
therefore were defined as elliptic gaussian functions given by Eqn. 4-48. The coefficients 𝑐

,𝑐
Q
and
𝑐
T
are the shape parameters that define the width and the length of the gaussian, as well as the
ship’s alignment with respect to earth fixed coordinate system based on the angle of rotation, 𝜃.

90
𝑍(𝑥,𝑦) = 𝐴∗exp[− (𝑐

(𝑋−𝑥0)
Q
+ 2𝑐
Q
(𝑋−𝑥0)(𝑌−𝑦0)
+ 𝑐
T
(𝑌−𝑦0)
Q
)]
Where;
4-48
𝑐

=
𝑐𝑜𝑠(𝜃)
Q
(2𝜎
J
Q
)
+
𝑠𝑖𝑛(𝜃)
Q
(2𝜎
L
Q
)
, 𝑐
Q
=−
sin(2𝜃)
(4𝜎
J
Q
)
+
sin(2𝜃)
4𝜎
L
Q

, 𝑐
T
=
𝑠𝑖𝑛(𝜃)
Q
(2𝜎
J
Q
)
+
𝑐𝑜𝑠(𝜃)
Q
(2𝜎
L
Q
)

Employment of the elliptic gaussian function as the pressure disturbance to define the vessels
has numerous advantages with some tolerable drawbacks, and all of these will be discussed here.
First of all, since MOST operates on a structured grid system, it has been seen in preliminary tests
that discontinuities can occur in the motion of the ships when they were expressed with discrete
functions (e.g. rectangular distribution), or with functions that have steep slopes on the edges (e.g.
top-hat function). Due to discontinuity in motion, abrupt increases in pressure at the grid cells near
the ship boundary were observed. The smoothness of the gaussian helped the neighboring grid
cells to experience pressure changes gradually as the ship translates. The other advantage of using
an elliptic gaussian function is in the transformation between body-fixed and global coordinate
systems. Likewise, with other types of functions, for individual vessels the pressure field had to
be interpolated from ship-fixed coordinates to global coordinate systems. The interpolation of the
pressure field at each time step caused a significant and unnecessary overhead; the system is
solved in a computationally cheaper way with elliptic gaussian functions since the rotation terms
are already accommodated in the shape parameters. Finally, if ships are assumed to have elliptic
shapes, then an efficient collision model that is described comprehensively in the previous section
can be used. The only shortcoming of using gaussian function as pressure disturbance is that the
draft of the ship is maximum at the center and decreases along its width and length towards the

91
edges, whereas the draft generally is distributed more uniform along the hull in actual ships.
Nonetheless, this issue can be tolerated considering the benefits acquired.

Figure 4-19: An example of the elliptic gaussian function employed, representing a vessel with the
length of 250 meters and width of 30 meters. a) Shows the perspective-view in 3D, b) Side-view
of the pressure distribution

4.4.1 Hypothetical Tsunami Scenario in Port of Long Beach
The Alaska-Aleutians Subduction Zone (AASZ) is a active and one of the most threating far-
field tsunami source regions for California. In recent history, there had been three remarkable
tsunamigenic earthquakes: the 8.2 Mw 1946 Aleutian Islands (Johnson and Satake 1997), 8.6 Mw
1957 Andreanof Islands (Johnson et al. 1994) and 9.2 Mw 1964 Prince William Sound (Kanamori
a)
b)

92
1970) earthquakes. Each of these earthquakes were followed by tsunamis in the Pacific which all
impacted California.
The Good Friday 1964 earthquake caused the most destructive tsunami both in the near field
and in the Pacific due to the enormous amount of energy released. Northern California was hit the
worst; the run-up reached 6.1 m in Crescent City and caused significant inundation. Moving south
along the California coast the wave height decreased to 3 m in San Francisco Bay and 3.8 m in
Santa Cruz (Wilson and Torum 1968). Across the state, 12 fatalities occurred in addition to 17
million dollars of damage (in 1964 USD)
The source characteristics of the hypothetical tsunami scenario considered in this study are
based on the 1964 earthquake, but it is located to the west of the 1964 rupture to increase the
energy beamed towards California. The estimated rupture area is 700 km long and 100 km wide
with an average slip of 25 m. Overall, the source corresponds to a magnitude 9.2 earthquake. The
details of the earthquake source parameters of this event were given by Barberopoulou et al.
(2011).
Based on this tsunami scenario, two different simulations are set-up focusing the Port of Long
Beach. In the first test, only one ship was placed within the port, while the second case had eight
identical ships located in and around the port (Figure 4-20). The parameters of the ships considered
in this part of the analysis is given in Table 4-1. Additionally, the results from these two
simulations were compared with a baseline simulation where no ships were present to see the
extent of the variations in the flow field induced by the vessels. A nested grid system similar to
the ones described in previous chapters was established here to propagate the tsunami waves from
Alaska to the Port of Long Beach. Here the grid size (Δ𝑥) was set to 5 meters instead of 10 meters

93
in the innermost computational grid to better resolve the vessels. In this final level of the nested
grids, the simulations were run for 10 hours of tsunami action.
Table 4-1: The main parameters of the vessels used in Port of Long Beach simulations
Length Meters 300
Beam Meters 60
Draft Meters 10
Deadweight Tons 150,000
Displacement Cubic meters 120,000

Figure 4-20: Southeast basin of the Port of Long Beach. a) Initial location of the only ship in test
case one; b) Initial locations of all eight identical ships are shown.

The results of these two sets of simulations were compared in Figure 4-21 with the plots
showing the current speeds estimated 100 minutes after the first arrival of the tsunami to the port.
When only one ship was included, the difference in the flow fields compared to the baseline
a) b)

94
simulation were local, and relatively weak. Whilst the larger eddy generated at the center of the
terminal in the baseline simulation also appeared in the one-ship case with a slight shift in south,
the variations outside the terminal are negligible. Contrary, in the multi-ship case the disparities in
the predicted flow fields became more profound. Both in and out of the terminal, the flow patterns
were altered drastically. As it has been put forward earlier, the presence of the vessels has a
considerable impact on the flow fields both in the vicinity of the terminal as can be seen from the
plots given in the left panel of Figure 4-21, which intensifies with the number of ships included in
the analysis.

Figure 4-21: Current speeds extracted 100 mins after the arrival of the first wave. Right panel
shows the current field from the baseline simulation; middle panel are the current speeds when
vessels are included; right panel absolute difference of current speeds predicted in both cases. The
figures on top are from one-ship simulation, while the bottom plots are from the multi-ship case.

95
4.4.2 Incident in Ishinomaki Harbor During 2011 Tohoku Tsunami
Like the rest of the Tohoku coast, the Ishinomaki region in Miyagi prefecture was hit hard by
the 2011 Tohoku Tsunami. In Ishinomaki city alone, ~4,000 casualties were reported, ~20,000
houses were completely destroyed, and entire port facilities were flooded (Takagi and Bricker,
2014).
The C.S Viceroy, a 177-meter-long bulk carrier was moored inside the Ishinomaki Harbor at
the time of the tsunami. A personal communication with Dr. Tatsuya Asai revealed that
approximately 66 minutes after the earthquake, around 16:00 JST, the ship broke all its moorings
and started to drift freely inside the terminal until it was carried out of the port by the tsunami. The
rough location of the ship at the time of the incident and the path it followed are shown in Figure
4-22.

Figure 4-22: Aerial photo of the Ishinomaki Port. The red cross shows the initial location of the
C.S. Victory. The red line indicates the approximate path of the ship after it broke free.

A realization of this incident is depicted here, using the dimensions and the physical properties
of the vessel that are listed in Table 4-2. Since the precise draft of the ship at the time of the incident

96
is unknown, in this analysis it was taken as 7.5 meters; the midpoint of the minimum and maximum
possible draft. The tsunami waves were brought to Ishinomaki Port through a nested grid system,
where the final grid has the horizontal resolution of 5 meters. The duration of the simulations in
this innermost grid was 5 hours. Initially, the C.S Viceroy was placed to the mooring spot in the
terminal where it was thought to be moored when the tsunami had arrived. Since the mooring lines
were not included in this study, the vessel’s motion had been restricted in the model for the first
66 minutes of the simulation, then it was released and let drift freely with the tsunami currents
until the end of the run.
Table 4-2: The main parameters of C.S Victory
Length Meters 177
Beam Meters 26
Draft Meters 6.1 – 9.0
Deadweight Tons 32,500
Displacement Cubic meters 20,000

Snapshots of the simulation showing the motion of the vessel as well as the tsunami current
fields, taken at the different times of the simulation after the ship was released are illustrated in
Figure 4-23. The figure shows that the speed of the tsunami currents entering to the terminal where
the ship was residing about ~3 𝑚/𝑠 around 𝑡 =60 mins. This was approximately the time when
the ship broke its mooring lines. A large eddy formed near the ship shortly after which pulled it
away from the berth and carried it south near the terminal entrance. At around 𝑡 =90 mins, the
next tsunami surge has arrived and pushed the ship back in. After drifting and spinning with the

97
flow in the port until 𝑡 ≈ 150 mins, finally it was taken out of the harbor and brought to the
southern boundary of the computational domain.

Figure 4-23: Snapshots of the current field and the vessel’s location at different instants.

(m/s)

98
Moreover, in Figure 4-24, the entire path of the vessel predicted by the model compared with the
data acquired from Dr. Asai. According to this figure, the model arguably delivered a good
realization of the incident based on the initial location, release time and draft assignments.
However, as briefly discussed earlier, these parameters were in fact unknown a priori. Also,
during the preliminary analysis, it has been noticed that if these parameters are changed in the
beginning of the simulation, then the ship can follow a different path. Therefore, based on this
observation, a sensitivity analysis was performed to better understand the role of the initial choice
of these “free” parameters in model outputs.

Figure 4-24: The comparison of the path of the vessel predicted by the model for the draft of 7.5
meters and release time of 66 mins, with the data provided by Dr. Asai.

Model
Data

99
The release time, draft and the initial location of the vessel in addition to the drag coefficient
for the lateral current loading (𝐶
&L
) were selected for the sensitivity analysis. While, in general,
the release time depends on the resultant forces on the hull along with the capacity of the vessel’s
mooring lines, the draft is determined by the total load (cargo+fuel). Similarly, the exact location
of the vessel is an unknown, and moving the center of mass even by 1-2 grid points corresponds
to 5-10 meters of physical distance in the model, which can also lead to significant changes in
exerted current forces. Although the drag coefficient (𝐶
&L
) is calculated empirically by Eqn. 4-19,
which was derived by curve fitting to experimental date and aimed to be used by engineers. Yet,
it can still vary depending on the hull shape and the flow conditions, and therefore was included
in to this sensitivity analysis. The release time varied between 55 mins and 95 mins with one-
minute increments, while the draft was increased by 0.1 meters at each trial from 6.1 to 9 meters.
The initial center of mass moved 20 meters laterally and longitudinally, which adds up to 16 trials,
and the drag coefficient ranged from 1 to 2.5.
The results of the sensitivity analysis are summarized in Figure 4-25 and Figure 4-26: The
points show where the vessel is ended up after 5 hours long simulations. Different colors represent
the different parameters included in the sensitivity analysis., where all the paths calculated by the
model for different initial values of different parameters, and the final destinations of the vessel
after 5 hours long simulations are shown respectively. Both figures tell that the model’s outputs
are very sensitive to the initial choice of parameters and even small changes in the release time,
draft, initial location and the drag coefficient lead to remarkably different trajectories. These results
indicate that the process is indeed random as no apparent trends or patterns can be identified from
the figures. Although the model is capable of predicting the path of vessels for a given set of input
parameter, unfortunately the randomness associated with the selection of the initial conditions

100
implies that a pure deterministic approach like the one presented here is of limited use due to
chaotic behavior of the system.

Figure 4-25: The paths of the vessel estimated by the model for different values of a) Release time;
b) Vessel draft; c) Initial location; d) Drag coefficient for lateral loading.
a)
b)
c)
d)

101

Figure 4-26: The points show where the vessel is ended up after 5 hours long simulations. Different
colors represent the different parameters included in the sensitivity analysis.
4.5 Conclusions
In this chapter a large vessel transport model was presented which was developed to predict
the motion of the large vessel within the ports under the influence of the strong currents generated
by tsunamis. The transport module and the hydrodynamic solver are fully coupled, which allows
flow and vessels to interact with each other. This aspect was ignored in earlier attempts and indeed
is the key feature of this study. Considering the ratio of the draft of the large vessels (~8-10 meters)
to the typical depths in ports (~10-15 meters), the interaction between the flow and the vessels gets
amplified. Therefore, for accurate representation of the flow field in the presence of the vessels
consideration of their interaction is crucial.

102
The vessels were included into the hydrodynamic computations as free surface pressure
distributions which also accounted for the form drag that the flow experiences due to vessels.
Additionally, the bottom stress term was modified locally at the grid points where the pressure
gradient was non-zero to include effect of the skin friction. Then, the capability of the modified
hydrodynamic model was assessed first using the data provided in Ertekin et.al., (1986) where the
waves generated by the moving pressure disturbances had been investigated. It has been found that
the model can sufficiently capture the amplitude and the phase of the leading wave for the pressure
disturbances moving at the critical speed. Furthermore, the model was tested for the case when the
flow passes a floating object and the results were compared to the numerical benchmark test
created in OpenFOAM. With the modifications made, MOST captured the generation, separation
and the advection of the eddies generated in the lee of the floating block quite well when compared
to the results obtained in OpenFOAM.
The transport model was based on the linear equations of motion and only three degrees of
freedom are considered; surge, sway and the yaw. Added mass and the first order damping terms
were included as the restoring hull forces, but the forces imposed by the ship’s controls, thrusters
and rudders were not. The reasoning behind this is the reports from past incidents emphasizing
that once the vessel is picked-up by the flow, the ship’s controls are of limited assistance in taking
back the ship under control. Thus, the forces due to currents were the only sources of the external
loading. Moreover, the present model was also equipped with a built-in collision solver founded
on the concept of the conservation of momentum and impulse. It can detect and evaluate the
collisions between two or more vessels as well as the collisions that take place between the vessels
and the port structures.

103
In the first test, the effects of a hypothetical tsunami arising from Alaska-Aleutians Subduction
Zone was investigated for ships moored in the Port of Long Beach, which is the most important
commercial hub in the US west coast with the total annual value of the export and import cargo of
$100 billion. Two different setups have been investigated here; one with a single ship moored
within the port, while eight ships were placed in and around the port in the second. Then, the results
from these two simulations were compared with the baseline simulation where no vessels were
included. This analysis revealed that the flow in and around the vessels is indeed affected by the
vessels and this gets intensified as the number of the vessels increases, which proves the initial
hypothesis of this research.
The last large-scale test was based on an incident reported in Ishinomaki Port, Japan that
happened during the 2011 Tohoku Tsunami. In this incident a ~180-meter-long bulk carrier broke
its mooring approximately 65 mins after the first arrival of the tsunami and was taken out of the
port after spinning and drifting freely inside the port. First, a realization of this event has been done
using the developed transport model, and then, a sensitivity analysis has been performed based on
this case to see the effects of the parameters that cannot be known precisely before the incident
like the draft of the vessel, or the exact time of the mooring line failure. The results from these
analyses showed that the model is capable of estimating the ship’s path for certain combinations
of input parameters, but at the same time the model results are very sensitive to the initial choice
of the input parameters, and the ultimate predicted path of the ship could entirely be different.

104
5. Conclusions and Future Work Ideas
5.1 Conclusions
The main objective of this research was to understand the transport potential of the currents
generated by tsunamis in the nearshore, particularly ports and harbors. Before investigating the
details of transport phenomena, a sensitivity analysis was conducted to understand the effects of
tides and wave directionality on the localized tsunami-induced currents. A main objective of this
research component was to see if the inclusion of tides in tsunami simulations has any meaningful
effect on the tsunami-generated currents. It has been found that the tide-induced variability in the
tsunami currents ranged from 10% to 25%, which is a variability lower than the typical forecast
accuracy of the current models. Then, role of the source location on the maximum nearshore
current speed predictions was assessed to see if the decisions about maritime advisories or
warnings might be made based only on the local wave height forecast. The results suggested that
the source region does not play a significant role in the prediction of maximum speeds and that
spatial variability in maximum currents is dominated by the presence of eddies, for which
deterministic simulation has limited use.
Next, a sediment transport model has been developed where the erosion and deposition rates
were calculated using empirical relations, whereas the transportation of entrained particles is
estimated using an advection-diffusion scalar transport equation, evaluated numerically. The
developed transport model performed particularly well in large scale applications and provided
good estimates of expected morphological changes during tsunamis within the harbors. The only
drawback of the present tool, which is common in all sediment transport models, is that the initial
assignment of the parameters either regarding the sediment properties or in deposition and erosion

105
formulae influences the results notably. Therefore, an initial calibration is needed, if available,
based on the data collected after past events before using the model for future planning purposes.
Finally, a large vessel transport model was presented that was shown to be able to solve the
interaction between the flow and the vessel reasonably. The integrated collision model can
evaluate the colliding contact between multiple vessels or a vessel and port structures in a
physically correct manner as long as the vessels have elliptical shapes. In general, the developed
model can provide reasonable predictions on the ship’s path, yet a sensitivity analysis conducted
in Ishinomaki Port revealed that these predictions are extremely sensitive to the initial conditions.
This finding indicates that the process in fact is chaotic and requires further assessment.
5.2 Suggested Future Works
The developed sediment model can further be improved by replacing the empirical formulae
of erosion and deposition rates with the local gradients of sediment flux. Calculating the depth
changes based on the gradients of the sediment fluxes rather than fully empirical formulations that
treats erosion and deposition processes separately, can provide better model estimates. An example
of this approach is given in McCall et. al., (2010). Additionally, as discussed in Chapter 3, an NSW
model is not fully capable of resolving complex flow fields in the nearshore created by the
tsunamis. Therefore, coupling the sediment transport model to a higher order numerical solver
(e.g. Boussinesq-type) suited with better physics can represent the flow conditions, and eventually
lead to more accurate predictions by the sediment transport model.
The large vessel transport model can as well be improved and extended in several ways. First
of all, the collision model can be modified so that it will be able to detect and evaluate the collisions
between not only ellipses but also any convex polyhedral. Once it has this ability it can serve well

106
in a debris transport model, in which objects with random shapes and sizes are expected to be
introduced. Additionally, as the results of the sensitivity analysis presented in Chapter 4 suggests,
the large vessel transport process is chaotic by nature. Therefore, instead of a pure deterministic
approach like the one presented in this study, hybrid models that combine deterministic methods
with stochastic or statistical techniques should be employed to come up with a model that can be
used as a guidance tool in the future.

107
References
Adrian, R. J., Christensen, K. T., & Liu, Z.-C. (2000). Analysis and interpretation of instantaneous
turbulent velocity fields. Experiments in fluids, 29(3), 275–290.
Apotsos, A., Gelfenbaum, G., & Jaffe, B. (2011). Process-based modeling of tsunami inundation
and sediment transport. Journal of Geophysical Research: Earth Surface, 116(1), 1–20.
https://doi.org/10.1029/2010JF001797
Baraff, D. (1989). Analytical methods for dynamic simulation of non-penetrating rigid bodies.
No ACM SIGGRAPH Computer Graphics (Sēj. 23, lpp. 223–232). ACM.
Barberopoulou, A., Borrero, J. C., Uslu, B., Legg, M. R., & Synolakis, C. E. (2011). A second
generation of tsunami inundation maps for the State of California. Pure and applied
geophysics, 168(11), 2133–2146.
Borrero, J. C., Lynett, P. J., Kalligeris, N., & Borrero, J. C. (2015). Tsunami currents in ports
Subject Areas Author for correspondence :
Borrero, J. C., Bell, R., Csato, C., DeLange, W., Goring, D., Dougal Greer, S., … Power, W.
(2013). 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.
https://doi.org/10.1007/s00024-012-0492-6
Borrero, J. C., Goring, D. G., Greer, S. D., & Power, W. L. (2015). Far-Field Tsunami Hazard in
New Zealand Ports. Pure and Applied Geophysics, 172(3–4), 731–756.
https://doi.org/10.1007/s00024-014-0987-4
Cao, Z., Pender, G., Wallis, S., & Carling, P. (2004). Sediment Bed, (July), 689–703.
Cao, Z., Pender, G., Wallis, S., Carling, P., Simpson, G., Castelltort, S., … Wu, T.-R. Y.-T. Y.
(2011). Sedimentologic and geomorphologic tsunami imprints worldwide - A review.
Coastal Engineering, 22(4), 83–92. https://doi.org/10.1016/S0012-8252(03)00018-7
David, C. G., Roeber, V., Goseberg, N., & Schlurmann, T. (2017). Generation and propagation
of ship-borne waves-Solutions from a Boussinesq-type model. Coastal Engineering, 127,
170–187.
Dawson, A. G., & Shaozhong, S. (2000). Tsunami Deposits. Pure and Applied Geophysics,
157(6–8), 875–897. https://doi.org/10.1007/s000240050010
Dengler, L., Uslu, B., Barberopoulou, a., Borrero, J., & Synolakis, C. (2008). The Vulnerability
of Crescent City, California, to Tsunamis Generated by Earthquakes in the Kuril Islands
Region of the Northwestern Pacific. Seismological Research Letters, 79(5), 608–619.
https://doi.org/10.1785/gssrl.79.5.608

108
Ersan, D. B., & Beji, S. (2013). Numerical simulation of waves generated by a moving pressure
field. Ocean Engineering, 59, 231–239.
Fennema, R. J., & Chaudhry, M. H. (1990). Explicit methods for 2-D transient free surface
flows. Journal of Hydraulic Engineering, 116(8), 1013–1034.
https://doi.org/10.1061/(ASCE)0733-9429(1990)116:8(1013)
Fossen, T. I. (1994). Guidance and control of ocean vehicles (Sēj. 199). Wiley New York.
Fritz, H. M., Phillips, D. A., Okayasu, A., Shimozono, T., Liu, H., Mohammed, F., … Takahashi,
T. (2012). The 2011 Japan tsunami current velocity measurements from survivor videos at
Kesennuma Bay using LiDAR. Geophysical Research Letters, 39(2), 1–6.
https://doi.org/10.1029/2011GL050686
Goto, K., Takahashi, J., Oie, T., & Imamura, F. (2011). Remarkable bathymetric change in the
nearshore zone by the 2004 Indian Ocean tsunami: Kirinda Harbor, Sri Lanka.
Geomorphology, 127(1–2), 107–116. https://doi.org/10.1016/j.geomorph.2010.12.011
Hahn, J. K. (1988). Realistic animation of rigid bodies. No Acm Siggraph Computer Graphics
(Sēj. 22, lpp. 299–308). ACM.
Higuera, P., Lara, J. L., & Losada, I. J. (2014). Three-dimensional interaction of waves and
porous coastal structures using OpenFOAM®. Part I: Formulation and validation. Coastal
Engineering, 83, 243–258.
Johnson, J. M., & Satake, K. (1997). Estimation of seismic moment and slip distribution of the
April 1, 1946, Aleutian tsunami earthquake. Journal of Geophysical Research: Solid Earth,
102(B6), 11765–11774.
Johnson, J. M., Tanioka, Y., Ruff, L. J., Satake, K., Kanamori, H., & Sykes, L. R. (1994). The
1957 great Aleutian earthquake. No Shallow Subduction Zones: Seismicity, Mechanics and
Seismic Potential (lpp. 3–28). Springer.
Kamrin, K., Rycroft, C. H., & Nave, J.-C. (2012). Reference map technique for finite-strain
elasticity and fluid–solid interaction. Journal of the Mechanics and Physics of Solids,
60(11), 1952–1969.
Kanamori, H. (1970). The Alaska earthquake of 1964: Radiation of long-period surface waves
and source mechanism. Journal of Geophysical Research, 75(26), 5029–5040.
Keen, A. S., Lynett, P. J., Eskijian, M. L., Ayca, A., & Wilson, R. (2017). Monte Carlo–based
approach to estimating fragility curves of floating docks for small craft marinas. Journal of
Waterway, Port, Coastal, and Ocean Engineering, 143(4), 4017004.
Kim, D. H., & Lynett, P. J. (2011). Turbulent mixing and passive scalar transport in shallow
flows. Physics of Fluids, 23(1). https://doi.org/10.1063/1.3541840

109
Kim, D.-H., Cho, Y.-S., & Kim, H.-J. (2008). Well-balanced scheme between flux and source
terms for computation of shallow-water equations over irregular bathymetry. Journal of
engineering mechanics, 134(April), 277–290. https://doi.org/10.1061/(ASCE)0733-
9399(2008)134:4(277)
Kim, D.-H., & Lee, S. O. (2011). Stable numerical model for transcritical flow and sediment
transport on uneven bathymetry. Journal of Hydraulic Engineering, 138(1), 46–56.
https://doi.org/10.1061/(ASCE)HY.1943-7900.0000473
Lee, S.-J., Yates, G. T., & Wu, T. Y. (1989). Experiments and analyses of upstream-advancing
solitary waves generated by moving disturbances. Journal of Fluid Mechanics, 199, 569–
593.
Liu, P. L.-F., & Wu, T.-R. (2004). Waves generated by moving pressure disturbances in
rectangular and trapezoidal channels. Journal of Hydraulic Research, 42(2), 163–171.
Lynett, P. J., Borrero, J. C., Weiss, R., Son, S., Greer, D., & Renteria, W. (2012). Observations
and modeling of tsunami-induced currents in ports and harbors. Earth and Planetary
Science Letters, 327–328, 68–74. https://doi.org/10.1016/j.epsl.2012.02.002
Lynett, P. J., Borrero, J., Son, S., Wilson, R., & Miller, K. (2014). Assessment of the tsunami-
induced current hazard. Geophysical Research Letters, 41(6), 2048–2055.
https://doi.org/10.1002/2013GL058680
Merritt, W. S., Letcher, R. A., & Jakeman, A. J. (2003). A review of erosion and sediment
transport models. Environmental Modelling and Software, 18(8–9), 761–799.
https://doi.org/10.1016/S1364-8152(03)00078-1
Moore, M., & Wilhelms, J. (1988). Collision detection and response for computer animation. No
ACM Siggraph Computer Graphics (Sēj. 22, lpp. 289–298). ACM.
Morgan, M. J. (1978). Dynamic positioning of offshore vessels.
Morrison, F. (2010). Data correlation for drag coefficient for sphere. Michigan Technology
University, Houghton, MI, 6(April), 1–2. Iegūts no
http://www.chem.mtu.edu/~fmorriso/DataCorrelationForSphereDrag2013.pdf
Okal, E. A., Fritz, H. M., Raad, P. E., Synolakis, C., Al-Shijbi, Y., & Al-Saifi, M. (2006). Oman
field survey after the December 2004 Indian Ocean tsunami. Earthquake Spectra,
22(SUPPL. 3). https://doi.org/10.1193/1.2202647
Okal, E. A., Fritz, H. M., Raveloson, R., Joelson, G., Pančošková, P., & Rambolamanana, G.
(2006). Madagascar field survey after the December 2004 Indian Ocean tsunami.
Earthquake Spectra, 22(SUPPL. 3). https://doi.org/10.1193/1.2202646

110
Okal, E. A., Sladen, A., & Okal, E. A. S. (2006). Rodrigues, Mauritius, and R??union Islands
field survey after the December 2004 Indian Ocean tsunami. Earthquake Spectra,
22(SUPPL. 3). https://doi.org/10.1193/1.2209190
Papanicolaou, A. N. T., Krallis, G., & Edinger, J. (2008). Sediment Transport Modeling Review
— Current and, 134(January), 1–14.
Paris, R., Wassmer, P., Sartohadi, J., Lavigne, F., Barthomeuf, B., Desgages, E., … Gomez, C.
(2009). Tsunamis as geomorphic crises: Lessons from the December 26, 2004 tsunami in
Lhok Nga, West Banda Aceh (Sumatra, Indonesia). Geomorphology, 104(1–2), 59–72.
https://doi.org/10.1016/j.geomorph.2008.05.040
Peskin, C. S. (2002). The immersed boundary method. Acta Numerica, 11, 479–517.
Platt, J. C., & Barr, A. H. (1988). Constraints methods for flexible models. ACM SIGGRAPH
Computer Graphics, 22(4), 279–288.
Postacchini, M., Othman, I. K., Brocchini, M., & Baldock, T. E. (2014). Sediment transport and
morphodynamics generated by a dam-break swash uprush: Coupled vs uncoupled modeling.
Coastal Engineering, 89, 99–105. https://doi.org/10.1016/j.coastaleng.2014.04.003
Rakha, K., Deigaard, R., & Broker, I. (1997). A phase-resolving cross-shore sediment transport
model for beach profile evolution. Coastal Eng., 31, 231–261.
https://doi.org/10.1016/S0378-3839(97)00008-2
Scheffers, A., & Kelletat, D. (2003). Sedimentologic and geomorphologic tsunami imprints
worldwide - A review. Earth-Science Reviews, 63(1–2), 83–92.
https://doi.org/10.1016/S0012-8252(03)00018-7
Scheffers, A., & Kelletat, D. (2003). Sedimentologic and geomorphologic tsunami imprints
worldwide - A review. Earth-Science Reviews, 63(1–2), 83–92.
https://doi.org/10.1016/S0012-8252(03)00018-7
Simpson, G., & Castelltort, S. (2006). Coupled model of surface water flow, sediment transport
and morphological evolution. Computers & Geosciences, 32(10), 1600–1614.
https://doi.org/10.1016/j.cageo.2006.02.020
Sorensen, R. M. (1973). Ship-generated waves. Adv. Hydroscience, 9, 49–83.
Steefel, C.I., Van Cappellan, P. (1998). Reactive transport modeling of natural systems. Journal
of hydrology, 209(1–4), 1–7.
Sugawara, D., Goto, K., & Jaffe, B. E. (2014). Numerical models of tsunami sediment transport -
Current understanding and future directions. Marine Geology, 352, 295–320.
https://doi.org/10.1016/j.margeo.2014.02.007

111
Takagi, H., & Bricker, J. D. (2014). Assessment of the effectiveness of general breakwaters in
reducing tsunami inundation in Ishinomaki. Coastal Engineering Journal, 56(4), 1450011–
1450018.
Takahashi, T., Shuto, N., Imamura, F., & Asai, D. (2001). Modeling sediment transport due to
tsunamis with exchange rate between bed load layer and suspended load layer. In Coastal
Engineering 2000, 1508-1519.
Terzopoulos, D., Platt, J., Barr, A., & Fleischer, K. (1987). Elastically deformable models. ACM
Siggraph Computer Graphics, 21(4), 205–214.
Titov, V., Kânoğlu, U., & Synolakis, C. (2016). Development of MOST for Real-Time Tsunami
Forecasting. Journal of Waterway, Port, Coastal, and Ocean Engineering, 142(6),
03116004. https://doi.org/10.1061/(ASCE)WW.1943-5460.0000357
Tomita, T., Imamura, F., Arikawa, T., Yasuda, Y., & Kawata, Y. (2006). Damage Caused by the
2004 Indian Ocean Tsunami on The Southwestern Coast of Sri Lanka. Coastal Engineering
Journal, 48(02), 99–116. https://doi.org/10.1142/S0578563406001362
U.S. Army Corps of Engineers (2005). Design: Moorings, unified facilities criteria.” UFC 4-159-
3, Washington, DC.
van Rijn, L. C., Walstra, D.-J. R., & van Ormondt, M. (2007). Unified View of Sediment
Transport by Currents and Waves. IV: Application of Morphodynamic Model. Journal of
Hydraulic Engineering, 133(7), 776–793. https://doi.org/10.1061/(ASCE)0733-
9429(2007)133:7(776)
Vennell, R. (2010). Resonance and trapping of topographic transient ocean waves generated by a
moving atmospheric disturbance. Journal of Fluid Mechanics, 650, 427–442.
Weiss, R., & Bourgeois, J. (2012). Reducing Tsunami Risk, 336(June), 1117–1118.
https://doi.org/10.1126/science.1221452
Wilson, B. W., & Torum, A. (1968). The tsunami of the Alaskan Earthquake, 1964: Engineering
evaluation. COASTAL ENGINEERING RESEARCH CENTER VICKSBURG MS.
Wilson, R. I., Admire, A. R., Borrero, J. C., Dengler, L. A., Legg, M. R., Lynett, P., …
Whitmore, P. M. (2013). Observations and Impacts from the 2010 Chilean and 2011
Japanese Tsunamis in California (USA). Pure and Applied Geophysics, 170(6–8), 1127–
1147. https://doi.org/10.1007/s00024-012-0527-z
Wilson, R., Davenport, C., & Jaffe, B. (2012). Sediment scour and deposition within harbors in
California (USA), caused by the March 11, 2011 Tohoku-oki tsunami. Sedimentary
Geology, 282, 228–240. https://doi.org/10.1016/j.sedgeo.2012.06.001

112
Wu, T. Y.-T. (1987). Generation of upstream advancing solitons by moving disturbances.
Journal of fluid mechanics, 184, 75–99.
Wu, W., & Wang, S. (2007). One-Dimensional Modeling of Dam-Break Flow over Movable
Beds. Journal of Hydraulic Engineering, 133(1), 48–58.
https://doi.org/10.1061/(ASCE)0733-9429(2007)133:1(48)
Xia, J., Lin, B., Falconer, R. A., & Wang, G. (2010). Modelling dam-break flows over mobile
beds using a 2D coupled approach. Advances in Water Resources, 33(2), 171–183.
https://doi.org/10.1016/j.advwatres.2009.11.004
Yamamoto, S., & Daiguji, H. (1993). Higher-order-accurate upwind schemes for solving the
compressible Euler and Navier-Stokes equations. Computers and Fluids, 22(2–3), 259–270.
https://doi.org/10.1016/0045-7930(93)90058-H

113
Appendix
This appendix includes the additional details about the grid setup used in the simulations, as
well as the source parameters of the hypothetical tsunami sources that were used to analyze the
effect of wave directionality on nearshore currents. With the help of supporting information
provided here along with the paper should be enough for other researchers to be able to reproduce
the results.
Figures A1 to A2 plot the nested bathymetry system used in all locations; Crescent City, Half
Moon Bay and San Diego Bay. Figures A4 to S7 show the time series of currents for tide-only,
tsunami-only and four tsunami+tide simulations extracted from Crescent City Harbor. The
locations of gauges where the time series were extracted at can be seen in Figure A3. Additionally,
Figures A8 and A9 show the comparison of modeled free surface elevations of tsunami only and
detided tide+tsunami simulations to measured data by the tide station in Crescent City during 2011
Tohoku Tsunami. Figure A10 illustrates the direction of the waves propagating towards the Pillar
Point Harbor in the open ocean, and their directions when entering to the innermost grid.
Finally, Tables A1 to A3 give the source parameters of all three hypothetical scenarios
extensively.

114

Figure A1. The cross Pacific propagation grid, and first level of nested bathymetries used in the
simulations.

115

Figure A2: Nested grid setup for, a) Crescent City, b) San Diego Bay and c) Half Moon Bay.
Footprints of each nested grid are identified with red boxes. Red plus signs in bottom figures mark
the tide stations.
a) b)
c)

116

Figure A3: Locations of synthetic gauges, at which the current time series shown in Figure S4-S7
were extracted.

Figure A4: Time series of current speeds extracted from a) Tide only simulation, b) tsunami only
simulation and c) Tide+Tsunami simulation for max tide case, at the locations shown in Figure
A3.
G4
G1
G2
G3
a)
b)
c)

117

Figure A5: Time series of current speeds extracted from a) Tide only simulation, b) tsunami only
simulation and c) Tide+Tsunami simulation for min tide case, at the locations shown in Figure
A3.

Figure A6: Time series of current speeds extracted from a) Tide only simulation, b) tsunami only
simulation and c) Tide+Tsunami simulation for mid-high tide case, at the locations shown in
Figure A3.

a)
b)
c)
a)
b)
c)

118

Figure S7. Time series of current speeds extracted from a) Tide only simulation, b) tsunami only
simulation and c) Tide+Tsunami simulation for mid-low tide case, at the locations shown in Figure
A3.

Figure S8. Comparison of the free surface elevations in Crescent City of Tide Station (blue),
tsunami only MOST simulation (orange) and tsunami+tide MOST simulation (yellow) for the first
5 hours of wave actions a) max tide case and b) min tide case. Location of tide station is shown in
Figure A2a.
a)
b)
c)

119

Figure A9. Comparison of the free surface elevations in Crescent City of Tide Station (blue),
tsunami only MOST simulation (orange) and tsunami+tide MOST simulation (yellow) for the first
5 hours of wave actions a) mid-high tide case and b)mid-low tide case. Location of tide station is
shown in Figure A2a.

120

Figure A10. Direction of incoming tsunamis (left panel), and near the study harbor (Pillar Point
Harbor) for a)Alaska, b)Chile and c)Mariana scenarios

121
Table A1: Source parameters for Alaska Scenario
Scenario Lon (
o
E) Lat (
o

N)
Slip
(m)
Strike (
o
) Dip
(
o
)
Depth
(km)
Length
(km)
Width
(km)
Rake
(
o
)
Alaska
Subfault 1 202.2610 55.1330 3.75 247 15 17.94 100 50 90
Subfault 2 201.8797 55.4908 3.75 246.2137 15 30.88 100 50 90
Subfault 3 203.6040 55.5090 3.75 240 15 17.94 100 50 90
Subfault 4 203.1521 55.8534 3.75 240.4869 15 30.88 100 50 90
Subfault 5 204.8950 55.9700 3.75 236 15 17.94 100 50 90
Subfault 6 204.4315 56.3016 3.75 235.669 15 30.88 100 50 90
Subfault 7 206.2080 56.4730 3.75 236 15 17.94 100 50 90
Subfault 8 205.7880 56.8210 3.75 235.4756 15 30.88 100 50 90
Subfault 9 207.5370 56.9750 3.75 236 15 17.94 100 50 90
Subfault 10 207.1091 57.3227 3.75 235.4119 15 30.88 100 50 90
Subfault 11 208.9371 57.5124 3.75 236 15 17.94 100 50 90
Subfault 12 208.4198 57.8213 3.75 234.6891 15 30.88 100 50 90

Table A2: Source parameters for Chile Scenario
Scenario Lon (
o
E) Lat (
o
N) Slip
(m)
Strike
(
o
)
Dip
(
o
)
Depth
(km)
Length
(km)
Width
(km)
Rake
(
o
)
Chile
Subfault 1 289.1231 -27.7826 23 14.16 14.54 11.96 100 50 90
Subfault 2 288.6469 -27.6762 23 14.16 8 5 100 50 90
Subfault 3 288.8943 -28.6409 23 13.19 14.63 11.96 100 50 90
Subfault 4 288.4124 -28.5417 23 13.19 8 5 100 50 90
Subfault 5 288.7113 -29.4680 23 9.68 14.72 11.96 100 50 90
Subfault 6 288.2196 -29.3950 23 9.68 8 5 100 50 90
Subfault 7 288.5944 -30.2923 23 5.36 14.81 11.96 100 50 90
Subfault 8 288.0938 -30.2517 23 5.36 8 5 100 50 90
Subfault 9 288.5223 -31.1639 23 3.8 14.9 11.96 100 50 90
Subfault 10 288.0163 -31.1351 23 3.8 8 5 100 50 90
Subfault 11 288.4748 -32.0416 23 2.55 15 11.96 100 50 90
Subfault 12 287.9635 -32.0223 23 2.55 8 5 100 50 90
Subfault 13 288.3901 -33.0041 23 7.01 15 11.96 100 50 90
Subfault 14 287.8768 -32.9512 23 7.01 8 5 100 50 90

122
Table A3: Source parameters for Mariana Scenario
Scenario Lon (
o
E) Lat (
o
N) Slip
(m)
Strike
(
o
)
Dip
(
o
)
Depth
(km)
Length
(km)
Width
(km)
Rake
(
o
)
Mariana
Subfault 1 146.6689 19.3123 8 164.48 25.07 21.41 100 50 90
Subfault 2 147.0846 19.4212 8 164.48 19.16 5 100 50 90
Subfault 3 146.9297 18.5663 8 172.07 22 22.1 100 50 90
Subfault 4 147.3650 18.6238 8 172.07 20 5 100 50 90
Subfault 5 146.9495 17.7148 8 175.11 22.06 22.04 100 50 90
Subfault 6 147.3850 17.7503 8 175.11 19.93 5 100 50 90
Subfault 7 146.9447 16.8869 8 180 25.51 18.61 100 50 90
Subfault 8 147.3683 16.8869 8 180 15.79 5 100 50 90
Subfault 9 146.8626 16.0669 8 185.18 27.39 18.41 100 50 90
Subfault 10 147.2758 16.0309 8 185.18 15.56 5 100 50 90
Subfault 11 146.7068 15.3883 8 199.05 28.12 20.91 100 50 90
Subfault 12 147.0949 15.2590 8 199.05 18.56 5 100 50 90
Subfault 13 146.4717 14.6025 8 204.35 29.6 26.27 100 50 90
Subfault 14 146.8391 14.4415 8 204.35 25.18 5 100 50 90
Subfault 15 146.1678 13.9485 8 217.45 32.04 26.79 100 50 90
Subfault 16 146.4789 13.7170 8 217.45 25.84 5 100 50 90

Abstract (if available)

Abstract The purpose of this research is to study the transport processes induced by the tsunami-generated currents in the nearshore. By this research, it is intended to significantly advance our ability to provide quantitative guidance to ports and harbors in the development of resiliency and recovery plans. ❧ First, the effects of dynamic tides and wave directivity on maxima of tsunami-induced currents were investigated. The goal of this part was to construct a basis for understanding the dynamic effects of tides and wave directivity on current-based tsunami hazards in a coastal zone through the application of numerical simulation tools for hazard mapping purposes. The effects of both were quantified through sensitivity analysis. It was concluded that the tide induced variability in maximum tsunami currents can be in the order of 25%, while wave directivity has negligible impact. ❧ Then, an integrated numerical tool was developed to predict the changes in the bottom topography inside the harbor basins due to sediment transport driven by the tsunami-induced currents. The recent tsunamis that were affected California also led to led to areas of extreme scour as well as areas of large deposition caused by the strong currents. Thus, the main goal here was to be able to provide quantitative guidance to the authorities and/or to stakeholders to develop counter measures against this phenomenon through numerical modeling. ❧ Finally, the motion of large vessels in ports under the effect of the tsunami-induced currents was investigated in this thesis. The recent observations revealed that the maneuvering capabilities of vessels can become very limited when they are captured by strong currents, and within eddies in particular. The vessels that break free and drift-off in and around the port during tsunamis generally exhibits random behaviors, and thus poses a great risk of hazard to the port and its vicinity. The numerical tool was designed to allow for the interaction between flow and the vessel. Since the draft of the vessel in shallow waters reduces to the effective flow depth considerably, this interaction gets intensified when ships are in restricted waters like ports. It was shown that the developed tool is be able to solve the interaction between the flow and the vessel reasonably The integrated collision model can evaluate the colliding contact between multiple vessels or a vessel and port structures in a physically correct manner as long as the vessels have elliptical shapes. In general, the developed model can provide reasonable predictions on the ship’s path, yet a sensitivity analysis conducted in Ishinomaki Port revealed that these predictions are extremely sensitive to the initial conditions. This finding indicates that the process in fact is chaotic and requires further assessment.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Wave induced hydrodynamic complexity and transport in the nearshore

PDF

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

PDF

Numerical analysis of harbor oscillation under effect of fluctuating tidal level and varying harbor layout

PDF

Tsunami-induced turbulent coherent structures

PDF

Computer modeling for wave oscillation problems in harbors and coastal regions

PDF

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

PDF

Investigating the debris motion during extreme coastal events: experimental and numerical study

PDF

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

PDF

Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves

PDF

Immersive computing for coastal engineering

PDF

On the simulation of stratified turbulent flows

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

A subgrid-scale model for large-eddy simulation based on the physics of interscale energy transfer in turbulence

PDF

Accuracy and feasibility of combustion studies under engine relevant conditions

PDF

Modeling and simulation of complex recovery processes

PDF

The role of rigid foundation assumption in two-dimensional soil-structure interaction (SSI)

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Wave method for structural system identification and health monitoring of buildings based on layered shear beam model

PDF

Understanding human-building-emergency interactions in the built environment

PDF

On the transport dynamics of miscible fluids in porous media under different sources of disorder

Asset Metadata

Creator Ayca, Aykut (author)

Core Title Understanding the transport potential of nearshore tsunami currents

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Publication Date 10/11/2019

Defense Date 08/24/2018

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

Tag currents,debris transport,harbors,large vessel transport,numerical modeling,OAI-PMH Harvest,Ports,sediment transport,source directivity,tide interaction,tsunami,wave modeling

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Lynett, Patrick Joseph (committee chair), Domaradzki, Julian Andrzej (committee member), Lee, Jiin-Jen (committee member)

Creator Email ayca@usc.edu,aykutayca@yahoo.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-91200

Unique identifier UC11675425

Identifier etd-AycaAykut-6820.pdf (filename),usctheses-c89-91200 (legacy record id)

Legacy Identifier etd-AycaAykut-6820.pdf

Dmrecord 91200

Document Type Dissertation

Format application/pdf (imt)

Rights Ayca, Aykut

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