A Tsunami in Lisbon
Where to run?
Daniel André Silva Conde
Dissertação para obtenção do Grau de Mestre em
Engenharia Civil
Júri
Presidente:
Prof. Doutor António Jorge Silva Guerreiro Monteiro
Orientador:
Prof. Doutor Rui Miguel Lage Ferreira
Co-Orientador:
Prof. Doutor Carlos Alberto Ferreira de Sousa Oliveira
Vogal:
Prof. Doutora Maria Ana Viana Baptista
Vogal:
Engenheira Maria João Martins Telhado
Outubro 2012
To my parents, Armando e Guilhermina. Thank you.
Aos meus pais, Armando e Guilhermina. Obrigado.
Acknowledgments
This dissertation was developed at Instituto Superior Técnico - Technical University of Lisbon, Portugal,
under the guidance of Professors Rui M. L. Ferrira and C. Sousa Oliveira. During this time, part of the
work was supported by CEHIDRO - Centro de Estudos de Hidrossistemas - with two research initiation
scholarships, also under the guidance of Professor Rui M. L. Ferreira.
For the great opportunities that were given to me with no hesitation and for all the support, knowledge
and motivation transmitted throughout this year, my sincere gratefulness to Professor Rui M. L. Ferreira.
And I really do mean sincere.
My many thanks to Ricardo for all the support and valuable help with STAV-2D. And of course to Edgar,
for keeping the morale high when odds weren’t.
Also, to Professora Maria Ana Baptista, my sincere thanks for all the kindness and availability. Without
her contribution this work wouldn’t be possible. To Professor C. Sousa Oliveira for the thoughtful insights
about various topics on seismicity and tsunamis.
To my dear colleagues, specially André Marques and Pedro Pinotes, for making those hardworking hours
as enjoyable as they could have been. Also to my fellow colleagues and friends at the Residência Engo
Duarte Pacheco for providing the most hilarious 5 years that one could have had.
My non-translatable-to-words gratitude to my parents. To my mother, Mira, and my father, Armando,
for the constant support, guidance and love throughout my entire existence.
To my grandfather, Joaquim, for the sensitive human being that he is and inspires to be. To my
grandmother, Celeste, for being the craziest grandmother on this earth. And I mean this as the truest
compliment. To my grandmother, Alcínia, for all the support and welcoming affections every single time
she laid eyes on me.
To my godparents, Fátima and Sérgio, for all the support, affection and enjoyable moments that we have
had together so far. Also, my affectionate gratitude to my older "godfather", Sérgio, whose memory will
last forever.
A very special place for the younger family members. To Zi and Lino for having the most inspiring
relation that can be seen between two brothers. Everyone feels better if you’re around. Also, to the rising
star in the game (football!), Miguel for the sincere feelings and affections. To Catarina for the special
i
connection we share. No matter how long or how far I am away from all of you, I will never cease to have
four brothers.
And that’s about it. Oh wait... (ahah!) my loving gratitude to my girlfriend, Dora, for making the last
5 years probably the best ones that I have had so far. Thank you for fulfilling a long lasting void.
It would take forever to thank everyone I owe some part of myself or my life. To all those I have not
mentioned, my apologies and sincere thanks.
ii
Resumo
Uma revisão recente do histórico de tsunamis em Portugal mostrou que o estuário do Tejo foi afetado
por vários tsunamis catastróficos ao longo dos últimos dois milénios. Dada esta herança histórica, e a
crescente consciencialização devida a fenómenos recentes com mediatização mundial, o presente trabalho
procura fornecer informação relevante quanto à exposição das zonas ribeirinhas do estuário do Tejo a um
tsunami semelhante ao ocorrido em 1 de Novembro de 1755.
Novas abordagens, teóricas e metodológicas, foram usadas para modelar um cenário de Tsunami, com
especial enfase para a propagação em terra. Um dos aspetos mais relevantes desta fase de um tsunami
é a sua capacidade de incorporar detritos, sejam eles de origem natural, como sedimentos ou vegetação,
ou artificial, como ruínas de infraestrutura construída ou outros equipamentos do quotidiano humano.
Um modelo conceptual específico para escoamentos deste tipo - debris flows - foi revisto e ampliado com
um flexível esquema de condições de fronteira. O modelo numérico contempla as mais recentes técnicas de
discretização matemática e uma implementação computacional eficiente. Os seus fundamentos teóricos
são apresentados a um nível introdutório enquanto que as novas funcionalidades são exaustivamente
explicadas, demonstradas e discutidas.
Um esforço muito considerável foi direcionado para a discretização da cidade de Lisboa. O ambiente construído deixa de ser um simples pano de fundo para visualizar resultados e passa a definir, explicitamente,
a superfície do terreno. Todas as estruturas relevantes foram adequadamente modeladas com modelos
digitais de terreno e malhas computacionais muitíssimo finas para a escala do problema em mão.
Os resultados são apresentados em vários níveis diferentes, desde perspetivas abrangentes a análises
detalhadas ao nível de cada arruamento. Nestas últimas, é fornecida informação pormenorizada quanto
aos tempos de chegada do tsunami, as suas propriedades hidrodinâmicas e os padrões de deposição de
detritos subsequentes à inundação.
Este estudo demonstrou que vários locais na frente ribeirinha do estuário do Tejo são especialmente
preocupantes no que a um cenário de Tsunami diz respeito. Este documento fornece informação relevante
para apoiar o planeamento das respostas de socorro e eventual evacuação da cidade, constituindo mais
um passo na promoção da segurança dos cidadãos de Lisboa.
Palavras-Chave:
Tsunamis, modelação matemática, estuário do Tejo, Lisboa
iii
iv
Abstract
A recent revision of the catalogue of tsunamis in Portugal has shown that the Tagus estuary has been
affected by catastrophic tsunamis numerous times over the past two millennia. Provided this historical
heritage, and the growing awareness due to recent events worldwide, the present work aims at providing
relevant data on the exposure of the Tagus estuary waterfront to a tsunami similar to the one occurred
on the 1st November of 1755.
New theoretical and methodological approaches to tsunami modeling were employed, with special emphasis given to propagation over dry land. One of the most relevant features of this stage of the tsunami is its
ability to incorporate debris, either natural sediment incorporated from the bottom boundary or remains
of human built environment. The increased mass and momentum of the run-up can inflict tremendous
damage and, regrettably, severe human losses.
A conceptual model specifically suited for describing debris flows is reviewed and extended with flexible
boundary conditions. The employed numerical tool features the most up to date mathematical discretization techniques and an efficient computational performance. The theoretical fundamentals of the
numerical model are presented at an introductory level while the new features implemented are given a
thorough explanation and discussion.
A very considerable effort was directed to Lisbon’s topology discretization. The built environment is no
longer seen as a simple background layer for result visualization, it actually explicitly defines the terrain
surface. All relevant structures are adequately modeled with very fine elevation models and computational
meshes, provided the scale of the problem at hand.
The results are provided at different levels, ranging from generic overviews to street-level analysis on
designated areas. At these locations detailed data on wave arrival times, hydrodynamic features of the
advancing inundation front and urban debris scattering patterns are provided.
This study has shown that several locations on the Tagus estuary waterfront are especially worrisome
in what respects to a tsunami impact scenario. The present document provides relevant information to
support the design of emergency response plans, comprising one further step in promoting the safety of
Lisbon’s citizens.
Keywords:
Tsunamis, mathematical modelling, Tagus estuary, Lisbon.
v
vi
Notation
A
Cell area
[L2 ]
C
Sediment concentration
[−]
CL
Homogenous, depth-averaged sediment concentration in layer L
[−]
Cf
Friction coefficient
[−]
c
Shallow water wave velocity
[ms−1 ]
c̃ik
Approximate c in k edge
[ms−1 ]
ds
Reference sediment diameter
[m]
(n)
ẽik
n eigenvector
[ms−1 ]
E
Flux vector
fsij
F~
Stress tensor
[P a]
Generic force
[N ]
F
Flux vector in x
G
Flux vector in y
g
Gravitic acceleration
H
Source terms vector
h
Fluid height
[m]
H
Average depth
[m]
hL
Thickness of layer L
[m]
hL
Fluid height on the left side of a shock
[m]
hR
Fluid height on the right side of a shock
[m]
~n
Unit normal to a plane
[m]
Ks
Manning-Strickler coefficient
[m1/3 s−1 ]
L
Wavelength (KdV notation)
[−]
Ll
Lower interface of a given layer
[−]
Lu
Upper interface of a given layer
[−]
p
Bed porosity
[−]
p
Hydrostatic Pressure
[P a]
PL
Depth-averaged hydrostatic pressure in Layer L
[P a]
qs
Solid discharge
[m3 s−1 ]
qs∗
Solid discharge capacity
[m3 s−1 ]
R
Friction source term vector
s
Specific sediment gravity
[−]
S
Shock speed
[ms−1 ]
[m]
[ms−2 ]
tik
Unit tangent to k edge, in natural rotation from nik
T
Bottom slope source terms numerical flux matrix
Tij
Depth integrated turbulent tensions tensor
[P a]
uφ
Velocity associated with the vertical mass flux
[ms−1 ]
uI
Interface velocity
[ms−1 ]
uL
Velocity on the left side of a shock
[m]
uR
Velocity on the right side of a shock
[m]
u∗
Friction velocity
[ms−1 ]
UiL
Depth averaged velocity in layer L, in the xi direction
[ms−1 ]
UL
Depth averaged velocity in layer L, in the x direction
[ms−1 ]
U
Independent variables vector
vii
V
Primitive variables vector
VL
Depth averaged velocity in layer L, in the y direction
[ms−1 ]
wi
Weight factor for cell i
[m]
ws
Sediment settling velocity
[ms−1 ]
Zb
Bed elevation
[m]
(n)
αik
(n)
βik
Wave strengths in k edge
[m]
Bottom source flux coefficient in k edge
[m]
δij
Kronecker delta
[−]
Depth-averaged turbulent kinetic energy rate of dissipation
[m2 s−3 ]
ΦLu
Net vertical flux across the lower interface
[ms−1 ]
ΦLu
Net vertical flux across the upper interface
[ms−1 ]
Φa,b
Net vertical flux from layer a to layer b
[ms−1 ]
φ
Generic function
[m3 s−1 ]
η
Free-surface elevation (KdV notation)
[−]
θ
Shields parameter
[−]
κ
Depth-averaged turbulent kinetic energy
[m2 s−2 ]
(n)
λ̃ik
n eigenvalue of
[ms−1 ]
λt
Wavelength
[m]
λL
Characteristics speed on the left side of a shock
[m]
λR
Characteristics speed on the right side of a shock
[m]
Λ
Adaptation length
[m]
µ
Dynamic viscosity of the fluid
[P a.s]
ν
Poisson coefficient
[−]
νT
Turbulent viscosity
[P a.s]
ρ
Mixture density
[kgm−3 ]
ρL
Depth-averaged density of the mixture on layer L
[kgm−3 ]
ρ(w)
Clean water density
[kgm−3 ]
sigmaij
Stress tensor
[P a]
τij
Turbulent stress tensor
[P a]
τb
Bed shear stress
[P a]
τy
Yield stress
[P a]
τv
Viscous stress
[P a]
τt
Turbulent stress
[P a]
viii
Acronyms
BC
Boundary Conditions
BKBC
Bottom Kinematic Boundary Condition
CFD
Computational Fluid Dynamics
CFL
Courant-Friedrichs-Lewy
CPU
Central Processing Unit
DEM
Digital Elevation Model
FEM
Finite Element Method
FDM
Finite Difference Method
FSDBC
Free Surface Dynamic Boundary Condition
FSKBC
Free Surface Kinematic Boundary Condition
FVM
Finite Volume Method
GIS
Geographical Information Systems
IC
Initial Conditions
KdV
Kortweg-de Vries (equations)
PDE
Partial Differential Equation
RH
Rankine-Hugoniot (conditons)
RP
Riemann Problem
STAV-2D
Strong Transients in Alluvial Valleys 2D
ix
x
Contents
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
i
Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vii
Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ix
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xv
1 Introduction
1
1.1
Motivation and framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Proposed objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3
Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2 Tsunamis
5
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.2
What is a tsunami? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.3
Generation mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.4
Tsunami dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.5
Tsunami consequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.5.1
The 1755 Lisbon tsunami . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.5.2
Recent tsunamis in the Indian Ocean and Japan . . . . . . . . . . . . . . . . . . .
13
3 Conceptual Model
3.1
3.2
3.3
Basic equations of fluid and sediment dynamics . . . . . . . . . . . . . . . . . . . . . . . .
17
3.1.1
Material derivative and transport theorem . . . . . . . . . . . . . . . . . . . . . . .
17
3.1.2
Conservation of mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
3.1.3
Conservation of momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
Governing equations for stratified flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.2.1
Conservation of mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.2.2
Conservation of momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
Closure equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
4 Numerical Model
4.1
17
Discretization scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
33
xi
4.2
4.1.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
4.1.2
Finite volume schemes: The case of the 2D shallow-water equations . . . . . . . .
34
4.1.3
Stability region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
4.1.4
Wetting-drying algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
4.1.5
Entropy correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.2.2
Rankine-Hugoniot conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.2.3
Shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
4.2.4
Expansion waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
4.2.5
Analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.2.6
Physical feasibility of the mathematical solutions . . . . . . . . . . . . . . . . . . .
43
4.2.7
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
5 Integration with Geographic Information Systems
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.2
Pre-Processing utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.2.1
Generating boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.2.2
Refining meshes
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
Post-processing utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
5.3.1
51
5.3
Resampler and raster converter . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6 A 1755 tsunami in today’s Tagus estuary topography
6.1
55
Pre-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
6.1.1
Lisbon and Tagus estuary topography . . . . . . . . . . . . . . . . . . . . . . . . .
55
6.1.2
Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
6.1.3
Boundary Conditions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
6.2
Brief framework for result presentation and discussion . . . . . . . . . . . . . . . . . . . .
61
6.3
Overview of a tsunami scenario in Lisbon . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
6.3.1
General description of tsunami propagation . . . . . . . . . . . . . . . . . . . . . .
62
6.3.2
Impact at Carcavelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
6.3.3
Impact at Caxias and Cruz Quebrada . . . . . . . . . . . . . . . . . . . . . . . . .
64
6.3.4
Impact at Belém and Algés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
6.3.5
Impact at Caparica and Trafaria . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
Street-level results for inland propagation . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
6.4.1
Results for Alcântara
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
6.4.2
Results for Downtown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
6.4.3
Results for Almada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
Morphological impacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.5.1
The relevance of flow-morphology interaction . . . . . . . . . . . . . . . . . . . . .
87
6.5.2
Vale do Zebro Royal complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.4
6.5
xii
49
6.5.3
Urban debris washed at Downtown and Alcântara . . . . . . . . . . . . . . . . . .
89
6.5.4
Morphological impacts on other locations . . . . . . . . . . . . . . . . . . . . . . .
93
7 Conclusions and recommendations
95
7.1
Main conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
95
7.2
Directions for future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
96
Bibliography
Appendices
97
101
xiii
xiv
List of Figures
2.1
Oscillatory and translatory waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
Three main stages of tsunami propagation . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3
Schematics of a seismic generation mechanism . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.4
A soliton or solitary wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.5
A regular wave evolving into a bore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.6
1755 copper engraving depicting the 1755 Lisbon earthquake and tsunami (City Museum)
12
2.7
Tide gauge measurements of the 2004 tsunami at several locations (Rabinovich & Thomson,
2007) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.8
2.9
13
Morphological impact of the 2004 Indian Ocean at Lhoknga, northwest coast of Sumatra,
and Phuhket, Thailand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
Pressure gage readings for the 2011 Tohoku tsunami (Maeda et al., 2011) . . . . . . . . .
15
2.10 Morphological impact of the 2011 tsunami at Natori and Fukushima, south of Tohoku
(ABC News, 2012) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3.1
Layer set featured in the conceptual model
. . . . . . . . . . . . . . . . . . . . . . . . . .
22
3.2
A non-material surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.3
Layer set featured in the conceptual model
. . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.1
A shock (compressive) at the boundary, ∂Ω . . . . . . . . . . . . . . . . . . . . . . . . . .
41
4.2
Physical admissibility with respect to the entropy condition . . . . . . . . . . . . . . . . .
44
4.3
Comparison results for the implemented boundary conditions . . . . . . . . . . . . . . . .
45
4.4
Steady conditions at time t = 100 [s] and t = 2000 [s] . . . . . . . . . . . . . . . . . . . .
47
4.5
Steady conditions at time t = 100 s (still the same at t = 2000 s) . . . . . . . . . . . . . .
47
5.1
Pre-processing samples: example of ArcMap isoline defined boundary (left) and conversion
to Gmsh (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
5.2
Example of mesh refinement procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
5.3
Point-in-polygon (PIP) problem examples . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
5.4
Examples of VTK file resampling results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
6.1
DEM Raster with 5x5 [m] resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
6.2
Bathimetry with 10x10 [m] resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
xv
xvi
6.3
Example of DEM built environment enhancement . . . . . . . . . . . . . . . . . . . . . . .
57
6.5
Mean annual flow duration curve (Ómnias - Sacavém hydrometric station) . . . . . . . . .
58
6.4
Warm-up setup and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
6.6
Supplied data for a 1755-like tsunami . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
6.7
Supplied data for a 1755-like tsunami (2) . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
6.8
Relevant locations for result presentation and discussion . . . . . . . . . . . . . . . . . . .
61
6.9
Tsunami propagation overview at the Tagus estuary . . . . . . . . . . . . . . . . . . . . .
62
6.10 Tsunami run-up at Carcavelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
6.11 Tsunami run-up at Caxias and Cruz Quebrada . . . . . . . . . . . . . . . . . . . . . . . .
64
6.12 Tsunami run-up at Belém and Algés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
6.13 Tsunami run-up at Caparica and Trafaria . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
6.14 One-way traffic zones (top) and tsunami approach to Alcântara (bottom; t=10’20”) . . . .
67
6.15 Detailed results for Alcântara, t = 12 minutes 20 seconds . . . . . . . . . . . . . . . . . .
69
6.16 Detailed results for Alcântara, t = 14 minutes 20 seconds . . . . . . . . . . . . . . . . . .
70
6.17 Detailed results for Alcântara, t = 16 minutes 20 seconds . . . . . . . . . . . . . . . . . .
71
6.18 Detailed results for Alcântara, t = 20 minutes 20 seconds . . . . . . . . . . . . . . . . . .
72
6.19 One-way traffic zones (top) and tsunami approach to Downtown (bottom; t=13’20”) . . .
73
6.20 Detailed results for Downtown, t = 15 minutes . . . . . . . . . . . . . . . . . . . . . . . .
75
6.21 Detailed results for Downtown, t = 15 minutes 40 seconds . . . . . . . . . . . . . . . . . .
76
6.22 Detailed results for Downtown, t = 16 minutes 20 seconds . . . . . . . . . . . . . . . . . .
77
6.23 Detailed results for Downtown, t = 17 minutes 20 seconds . . . . . . . . . . . . . . . . . .
78
6.24 Detailed results for Downtown, t = 20 minutes 20 seconds . . . . . . . . . . . . . . . . . .
79
6.25 Detailed results for Cacilhas, t = 14 minutes 30 seconds . . . . . . . . . . . . . . . . . . .
81
6.26 Detailed results for Cacilhas, t = 16 minutes 10 seconds . . . . . . . . . . . . . . . . . . .
82
6.27 Detailed results for Cacilhas, t = 17 minutes 20 seconds . . . . . . . . . . . . . . . . . . .
83
6.28 Detailed results for Cacilhas, t = 22 minutes 50 seconds . . . . . . . . . . . . . . . . . . .
84
6.29 Detailed results for Cacilhas, t = 24 minutes 20 seconds . . . . . . . . . . . . . . . . . . .
85
6.30 Detailed results for Cacilhas, t = 27 minutes 40 seconds . . . . . . . . . . . . . . . . . . .
86
6.31 Tsunami approaching the Royal Complex at Vale do Zebro, t=25’ . . . . . . . . . . . . . .
87
6.32 Tsunami run-up at Coina river (Vale do Zebro) . . . . . . . . . . . . . . . . . . . . . . . .
88
6.33 Tsunami impact at the Royal complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
89
6.34 Idealized setup for urban debris, mainly composed of vehicle-like solid materials . . . . . .
90
6.35 Mobile material distribution along the riverfront . . . . . . . . . . . . . . . . . . . . . . .
90
6.36 Debris scattering patterns in Alcântara at the aftermath of the tsunami . . . . . . . . . .
91
6.37 Debris scattering patterns in Downtoen at the aftermath of the tsunami . . . . . . . . . .
91
6.38 Comparison between debris and clearwater formulations . . . . . . . . . . . . . . . . . . .
92
6.39 Morphological impacts on remaining locations . . . . . . . . . . . . . . . . . . . . . . . . .
93
Chapter 1
Introduction
1.1
Motivation and framework
Often considered as rare events when compared to other natural hazards tsunamis have an enormous
devastating potential. The frequency concept associated to these phenomena is to be understood within
consistent timescales, namely the geological timescale. At this level, tsunamis and earthquakes, as well
as many other geophysical occurrences such as volcanic eruptions or major landslides, become relatively
common events whose predictability, at finer timescales, becomes practically impossible. Other difficulties
arise when studying tsunamis since one of the main uncertainties lies in the potential of an oceanic
earthquake to generate tsunamis. Not all oceanic earthquakes occur in such a way that the seabed gets
vertically displaced and vertical momentum is transferred to the water column.
All these variables highlight a much required awareness of the stochastic nature of tsunamis and other
geophysical phenomena, in general. Furthermore, seismic risk mitigation is presently a standard in most
structural design regulations and public safety procedures. In coastal countries like Portugal, where
tectonic activity is likely to concentrate in oceanic faults, accounting for seismic risks while negletcing
potential tsunami impacts seems, at least, inconsistent and incautious.
The large ratio between seafront extension and continental territory made Portugal hugely dependant
on its maritime Economic Exclusive Zone. Whether it is due to turism, fishing activities or transoceanic
transport, the economic contribution of the seafront to the global wealth of the country has been very
significant and drove most of Portuguese towards the coast. Over centuries, this generalized demographic
behaviour ultimately converged into a country with two noticeably distinct development levels: a relatively unpopulated countryside on interior regions contrasting with densely populated and commercially
active coastal cities. Although the last decades had a significant contribution in fading such dissimilarities, drawn by major communication and technological progresses, the heritage left by recent centuries is
unavoidable: most of the population, infrastructure and valuable historical assets are located in coastal
areas. This emphasizes the relevance of seismic risk mitigation procedures that involve structural damage
as well as inundation extents.
1
1. Introduction
Recent revisions of the Portuguese catalog of tsunamis (Baptista & Miranda, 2009) have also shown that
the seafront has been struck by numerous events over the past two millenia.
Presently, numerous shallow water models are used within the scientific community for tsunami propagation (Imamura et al., 2006). While these include all the relevant components for offshore propagation,
their ability to account for very specific features in waterfront and overland propagation is limited. One
of the most relevant features of tsunami propagation over solid boundaries is its ability to incorporate
debris, either natural sediment incorporated from the bottom or remains of human built environment.
Despite acknowledged, this feature is commonly neglected in most tsunami inland propagation studies,
probably by virtue of a relatively incomplete understanding of the involved phenomena when sediment
transport imply significant changes in the rheological behaviour of the flow.
1.2
Proposed objectives
A computational fluid dynamics (CFD) model based on the finite volume method (FVM) (Canelas,
2010) - STAV2D (Strong Transients over Alluvial Valleys 2D) - featuring enhanced sediment transport
capabilities (Canelas, 2010) and up to date numerical discretization schemes (Murillo, 2006) was employed
to analyze a tsunami scenario in the Tagus estuary. The tsunami parameters resemble the widely known
event that struck Lisbon on the 1st November of 1755.
Not only is the present work based on improved conceptual and numerical models, but also, the physical
environment discretization was taken to a level not known in other studies. The human built environment
was detailed at a very refined scale, featuring buildings, streets and even urban mobile debris.
The main objective of this work lies in providing detailed data on the exposure of designated areas along
the Tagus estuary riverfront. Lisbon’s Downtown, Alcântara and Cacilhas were thoroughly discretized
in order to perform a street level description, a significant improvement over the present knowledge of a
tsunami impact on these locations. Morphological changes and debris scattering were also analyzed by
taking full advantage of the solid transport capabilities of the numerical model.
A set of secondary objectives was also delineated, consisting of computational work on the CFD code:
a new scheme of boundary conditions was implemented and pre/post-processing toolkits were developed
to promote functional integration with commercial Geographical Infromation Systems (GIS).
1.3
Structure
The first chapter is devoted to a brief bibliographical review on the subject of Tsunamis. Particular
emphasis is given to the hydrodynamic concepts and numerical modeling trends. Historical context is
also provided on a very relevant event for the present work: The 1755 Great Lisbon Earthquake and
Tsunami. Recent events that were given great interest in worldwide media and scientific communities
2
1.3. Structure
are also revised.
Theoretical background on the conceptual model that supports the present work and a revision of the employed closure models are provided in Chapter 2, with special care being given to particular mathematical
formulations inherent to multiphase flows.
A brief description of the discretization scheme is presented in the following chapter. The implemented
boundary condition scheme is thoroughly explained and its performance assessed with a set of functional
tests.
The fourth chapter introduces the developed toolkit for pre and post-processing tasks. It is particularly
useful for future reference and also highlights the expected ease of use of these tools.
The last chapter features a comprehensive description of the simulated scenario. All details concerning
the setup and outsourced inputs for these simulations are provided, including boundary and initial conditions and the employed topographic datasets. The results are displayed at different levels: a tsunami
scenario overview on the whole Tagus estuary, a detailed street-by-street analysis on designated areas
and morphological impact assessment.
3
1. Introduction
4
Chapter 2
Tsunamis
2.1
Introduction
Human communities living near the ocean have always experienced the devastating effects of a tsunami.
These waves have surely claimed hundreds of thousand of lives and have caused incalculable material
losses. However, since they are not frequent in industrialized western countries, they have not been given
relevant media attention. The morning of the December 24th changed this forever when an earthquake of
magnitude 9.1 stroke the coast of Sumatra in Indonesia. The tremendous energy released in this oceanic
earthquake generated a vertical displacement of the seabed, rising the sea surface to form a devastating
tsunami that killed 280,000 people and affected the lives of several million others (Helal, 2008).
While people living along the coastline of Sumatra had little time and no warning on the approaching
tsunami, those living along the further coasts of Thailand, Sri Lanka, India and East Africa had plenty
of time to escape to higher grounds. Unaware, they stayed at lower ground, as no tsunami warning
system was implemented in the Indian Ocean at that time. Unfortunately, it has taken a disaster of
such magnitude to awake the scientific communities to this silent but fearful natural hazard. The present
chapter aims to describing some of the most common physical features and modeling trends on tsunamis.
2.2
What is a tsunami?
The term tsunami originates from the Japanese tsu and nami meaning wave and harbor, respectively. It
may have been originally used by fishermen, who would sail out, encountering no visually unusual waves
while out at sea, and later returned to find an overwhelmed hometown.
Tsunamis are waves with large wavelengths that are originated by some kind of sudden disturbance in
standing
1
water and are characterized as shallow-water waves. These are different waves from the ones
most commonly observed on a beach, which are caused by the wind blowing across the ocean’s surface.
1 By standing water one means a water mass whose velocity is much lower than the velocity of a propagating tsunami,
√
which is close to that of shallow-water waves celerity ( gH).
5
2. Tsunamis
Wind-generated waves usually have periods of five to thirty seconds and a wavelength of 100 to 200 [m]
whereas tsunamis can have periods and wavelengths that are much higher (Helal, 2008).
A major difference between these waves and tsunamis is the mass transport carried by the wave. While
wind-generated waves induce nearly-closed trajectories 2 to the particles a tsunami displaces them several
meters from their initial positions. The volume of water that gets transported and the speeds reached
by tsunami waves allows them to carry enough energy to wipe out entire towns and cities (Helal, 2008).
Wind-generated waves are called oscillatory waves whereas tsunamis are translatory waves (Figure 2.1).
Oscillatory Wave
Translatory Wave
Stokes Drift
Figure 2.1: Oscillatory and translatory waves
Tsunami propagation comprises three main stages: oceanic propagation, near shore shoaling and inland
propagation (Figure 2.2 on page 7). At the first stage, within deep ocean waters, the period of a tsunami
can range from 5 minutes to 2 hours and the wavelength, λt , can be as high as 500 [km].
Waves are governed by shallow-water equations if the ratio between the wavelength and the depth of the
ocean is over 20. Given that the maximum mean depth of the sea-floor, H, is about 5 km, the ratio of the
flow depth to the wavelength is small: H/λt ≈ 0.001, a common value, for instance, in river flows. Hence,
tsunamis will always behave as shallow-water waves no matter where on the globe they are being modeled
√
or observed. The velocity of a shallow-water wave is also equal to gH, where g is the acceleration of
gravity (9.8 [m/s2 ]) and H is the depth of the water column.
Thus, in very deep water, a tsunami will propagate at high speeds: in the Pacific Ocean, for instance,
the average depth is approximately 4000 [km], which means that the velocity of a tsunami in this water
mass is nearly 200 [m/s]. Considering that the distance across this ocean, for instance from Seattle to
Hong Kong, is roughly 8000 [km] this would yield a total of 11 hours for travel time.
2 Trajectories are not exactly closed because the final position of the particles is not coincident with the initial position.
This is called the Stokes Drift (Stokes, 1847). Furhter details on surface wave dynamics can be found in the Appendix 1.
Nonetheless, the mass transport in this drift is negligible when compared to that of a tsunami.
6
2.3. Generation mechanisms
Offshore tsunami amplitudes are likely below 1 [m] and wavelengths are certainly a few hundred kilometers
long. Given that the order of magnitude of the hydraulic gradient can be given by the wave amplitude
to wavelength ratio, this leads to practically null energy dissipation: the hydraulic gradient is 10−5 for a
rather small wavelength of 100 [km]. This yields basically no head loss for a tsunami propagating over
deep waters, reaching distant shores with high energy levels.
As a tsunami approaches shallower waters near shore the depth of the water decreases, meaning that the
velocity of the tsunami will also decrease. However, its energy and mass are maintained, suggesting that
the height of the tsunami must grow enormously. It is at this stage that the tsunami becomes visible,
resembling the atrocious footage spread out on the media in 2004 (Indian Ocean) and 2011 (Japan).
The third stage of a tsunami deals with breaking as it approaches shore and depends greatly on the
bathymetry. Overland tsunami propagation may form bores: abrupt wave forms mathematically represented as discontinuities. These discontinuous and highly turbulent flows are likely to incorporate debris,
whether natural sediment or human built environment remains. The dynamics of this final stage of a
tsunami is somewhat similar to that of flood waves caused by dam breaking, dyke breaking or overtopping.
Very high concentrations of debris are carried on the water column and may change the hydrodynamic
properties and rheological behaviour of the fluid (Dias, 2006).
Translatory wave
low amplitude
high wavelength
mass is conserved
bores may form
amplitude increases
mean sea level
velocity
shore
continental shelf
sea bed
Figure 2.2: Three main stages of tsunami propagation
2.3
Generation mechanisms
Tsunamis can are mainly generated by seismic activity (seismicity), landslides or volcanic activity. The
seismic generation mechanism is responsible for most tsunami occurrences and, as a major source, must be
discussed with further detail. Tsunamis can be generated by seismicity when the sea floor abruptly rises
and vertically displaces the overlying seawater. Tectonic earthquakes are a particular kind of earthquake
that are associated with the continental drift. When these earthquakes occur in the ocean, the water
above the deformed area is displaced from its initial position (Dias, 2006). Namely, a tsunami can be
generated when thrust faults, associated with convergent or destructive plate boundaries, suddenly slip
and displace large volumes of water due to the vertical component of momentum involved (see Figures
2.3a to 2.3d).
7
2. Tsunamis
Mathematical modeling of this mechanism requires a displacement field to be parametrized at the sea
bottom. Volterra (1907) originally solved this problem in the case of an entire elastic space. Hence, and
for simplicity’s sake, that will be the approach herein presented.
Considering the thrust mechanism as a fracture, the discontinuities in the displaced components across
the fractured surface can be described using the theory of elasticity, namely surfaces across which the displacement field is discontinuous (Steketee, 1958). For simplicity’s sake some assumptions are introduced.
The curvature of the earth, its gravity, temperature, magnetism and non-homogeneity are neglected. A
semi-infinite, homogeneous and isotropic medium is considered instead. Further details on how the earth’s
curvature, topography, homogeneity and isotropy influence the displacement fields can be found on the
works published by McGinley (1969), Ben-Mehanem et al. (1970a,b) and Masterlark (2003). Furthermore,
one assumes that the the linear elasticity theories are valid within the fractured region.
Let O be the origin of a Cartesian coordinate system, xi the Cartesian coordinates with i = 1, 2, 3, and ~ei
a unit vector in the positive xi direction. A force F~ = F ~ek at O generates a displacement field uki (P, O)
at point some point P , given by,
uki (P, O) =
F
8πµ
δik
∂2r
∂2r
−α
2
∂xn
∂xi ∂xk
(2.1)
where δik is the Kronecker delta, r is the distance from P to O and α is a coefficient defined as,
α=
1
λ+µ
= (1 − ν)
λ + 2µ
2
(2.2)
where λ and µ are the Lamé’s constans and ν is the Poisson’s ratio
3
. Hooke’s Law easily yields the
stresses due to the displacement fields in equation (2.2)
k
σij
= λδij
∂uk
+µ
∂xk
∂ui
∂uj
+
∂xj
∂xi
(2.3)
−
The dislocation is defined as a surface Γ across which one finds a discontinuity ∆ui = u+
i − ui in the
displacement fields. This yields the following formula (Volterra, 1907) describing the displacement field
for a particular point (x1 , x2 , x3 )
uk (x1 , x2 , x3 ) =
1
F
ZZ
k
∆ui σij
· nj dS
(2.4)
Γ
. In general, the displacement field represented by equation (2.4) is reproduced in the ocean free surface
(Figure 2.3c on page 9), posing the initial conditions for tsunami propagation (Imamura et al., 2006).
Tsunamis often have a small amplitude offshore which is why they generally pass unnoticed at sea, forming
only a slight swell usually less than half a meter above the normal sea surface. A tsunami can occur in
3 Transverse
8
to axial strain ratio when a material is compressed.
2.3. Generation mechanisms
any tidal state and even at low tide can still severely inundate coastal areas. The 1964 Alaska earthquake
(Magnitude, Mw, 9.2), the 2004 Indian Ocean earthquake (Mw 9.2) and the 2011 Tohoku earthquake
(Mw 9.0) are recent examples of powerful thrust earthquakes that generated powerful tsunamis, capable
of crossing entire oceans and still inflict severe devastation in coastal areas. Other potential generation
mechanisms may occur, although with fewer recorded occurrences. These include landslides, explosive
volcanic eruptions, large atmospheric depressions and even meteor impacts.
Tsunamis generated by landslides are called sciorrucks. These phenomena rapidly displace large water
volumes, as energy from falling debris or expansion is transferred to the water at a rate faster than it can
absorb. Their existence was confirmed in 1958, when a giant landslide in Lituya Bay, Alaska, caused the
highest wave ever recorded, which had a height of 524 [m]. The wave didn’t travel far, as it struck land
almost immediately. Two people fishing in the bay were killed, but another boat amazingly managed
to ride the wave (Fine et al., 2006). Scientists named these waves megatsunami and discovered that
extremely large landslides from volcanic island collapses can also generate these kind of waves.
Meteotsunamis are generated by deep depressions that cause tropical cyclones, generating a storm surge
often called a meteotsunami that can raise tides several meters above normal levels. The displacement
comes from low atmospheric pressure within the centre of the depression.
SWL - Still Water Level
SWL - Still Water Level
Overriding Plate
Subducting Plate
Overriding Plate
Subducting Plate
(a) Converging plates before earthquake
(b) Overriding plate bends under strain
Overriding Plate
Subducting Plate
Overriding Plate
Subducting Plate
(c) Accumulated energy is released
(d) Tsunami propagates radially
Figure 2.3: Schematics of a seismic generation mechanism
9
2. Tsunamis
2.4
Tsunami dynamics
A simple mathematical formulation to demonstrate the physics behind a tsunami wave is herein presented Helal (2008). In this approach a tsunami is conceptualized as a non-linear bi-dimensional wave
propagating over a horizontal ground.
The governing law for this solution is the Korteweg-de Vries equation (KdV). It is a nonlinear, dispersive
partial differential equation for a function of two real variables, space, x, and time, t. For some function
φ the Korteweg-de Vries equation is stated as,
∂t φ + ∂x3 φ + 6φ∂x φ = 0
(2.5)
where the constant 6 in front of the last term is merely conventional: multiplying t, x, and φ by constants
can be used to vary the constants preceding each of the featured terms.
For hydrodynamic applications this equation is often presented in the form,
∂η
∂η
3 ∂η
ξ ∂3η
+
+ η
+
=0
∂t
∂x 2 ∂x 6 ∂x3
(2.6)
and featuring the following non-dimensional variables,
η=
δ−h
;
H
=
h
;
H
ξ=
h2
L2
(2.7)
where is the wave amplitude corresponding to the depth h; L is the wavelength; η(x, t) is the relative
movement of the free surface relating to the maximum height and ξ is the squared depth to wavelength
ratio. This form of the KdV equation is an numerical simulation standard when it comes to offshore
propagation and its solutions are provided in the form of solitons (Figure 2.4) or other dispersive waves.
L
H
h
mean sea level
δ
sea bed
Figure 2.4: A soliton or solitary wave
The dispersive term (fourth in equation (2.6)) counterbalances the quasi-linear terms (second an third
terms in equation (2.6)) allowing for the self-similar propagation of the wave-form shown in Figure 2.4.
As the wave approaches shore non-linear terms will become more important and the wave will cease to
be decribed by the KdV equation.
10
2.5. Tsunami consequences
When a tsunami reaches the seashore and propagates overland ( last stage in Figure 2.2) it may develop
into an abrupt wave form, or bore, which is mathematically considered as a discontinuity. Bores are
non-linear wave forms withour relevant dispersive terms, especially if debris are incorporated into the
water column.
t0
t0+dt
t0+2dt
Figure 2.5: A regular wave evolving into a bore
In this case, the dynamics of the tsunami is governed by shallow-water equations. In fact, dispersive
effects are often disregarded and shallow-water equations are frequently used to describe the propagation
of tsunamis in open waters (Imamura & Shuto, 1990).
The shallow-water equations will not be described in this section. They constitute the core of the
mathematical simulation tools employed in this dissertation and will be carefully derived in Chapter 3,
p. 17 to 21.
2.5
2.5.1
Tsunami consequences
The 1755 Lisbon tsunami
Around 9:40 AM on the 1st November of 1755, the All Saints’ Day Catholic holiday, Lisbon was struck
by a violent earthquake with a magnitude ranging from 8.5 Mw to 9 Mw 4 , (Abe, 1989) with an epicenter
in the Atlantic Ocean
5
that became known as the Great 1755 Lisbon Earthquake. This earthquake was
followed by several fires and most importantly, for the purpose of this dissertation, by a major tsunami.
At that time Lisbon was inhabited not only by Portuguese people but also by a considerable foreign
population, mainly French, Spanish, English and German. This contributed to the wide diversity of
publications, in both Portuguese and English, regarding the Great 1755 Lisbon Earthquake and made it
the most well documented historic tsunami event (Baptista et al., 1998b).
The following are amongst the most relevant historic reports (Baptista et al., 1998a, 2006),
"There was another great shock after this that pretty much affected the river, but, I think not
so violent as the preceding, though several people afterwards assured that as they were riding on
horseback, in the great road to Belem, one side of which lays upon the river, the waters rushed
4 Moment
5 Located
Magnitude Scale, Mw
south of the west coast of Algarve and northeast of the Archipelago of Madeira(Baptista et al., 1998b)
11
2. Tsunamis
in so fast that they were forced to gallop as fast as possible to the upper grounds for fear of being
carried away" (British Accounts, 1755)
"A large quay, piled up with goods near the Custom House sunk the first shock, with about 600
persons upon it, who all perished" (British Accounts, 1755)
"(...) a lot of people run to the river bank trying to escape from the ruins. Suddenly the sea came
in through the bar and flooded the river banks..." (Moreira de Mendonça, 1758).
"(...) waters with their ebb and flow flooded the Customs, the Square and the Vedoria".
Figure 2.6: 1755 copper engraving depicting the 1755 Lisbon earthquake and tsunami (City Museum)
To avoid interpretation and subsequent distortions Baptista et al. (1998a) carried a direct analysis of
historical data between 1755 and 1759 where all the documents were original from the 18th century and
were analyzed in their original language. This study pointed out very useful aspects about the earthquake
and tsunami that destroyed Lisbon on that morning of the November 1st .
The earthquake shook the city at about 9h30 a.m. and lasted about 9min. When most people were
heading to the riverfront in a desperate effort to avoid the fires and ruining buildings in the inner city
and also out of curiosity about the unusual low tide, which was in fact a drawback and the only warning
of an incoming tsunami, a huge water wall hit the riverfront. The tsunami traveled for about 90 minutes
in the open sea and should have hit the city at about 11h00 a.m. climbing the shores up to a mean
run-up level of 5 [m] and penetrating the city as deep as 250 [m]. The death toll only due to the tsunami
is estimated in around 900 people whereas the earthquake itself should be held responsible for a number
lives ranging from 10.000 up to 100.000 (Álvaro S. Pereira, 2006).
About eighty-five percent of the city’s built environment collapsed, including famous palaces and libraries
constituting some of the most brilliant examples of the unique Manueline architectural style. Even
buildings that had suffered little damage from the earthquake were destroyed in the subsequent fires.
The Royal Ribeira Palace, which was located near the riverfront in the modern square of Terreiro do
Paço, was destroyed by the earthquake and tsunami. Inside, the 70.000 volume royal library, also filled
12
2.5. Tsunami consequences
with hundreds of works of art by notorious painters and sculptors of that time were forever lost. The
scientific and historian communities were also deprived of priceless intelligence on the early campaigns
led by the famous navigators Vasco da Gama, Fernão de Magalhães and Pedro Alvares Cabral.
All the major churches in Lisbon were ruined by the earthquake, namely the Lisbon Cathedral, the
Basilicas of São Paulo, Santa Catarina, São Vicente de Fora, and the Misericórdia Church. The Royal
Hospital of All Saints in the Rossio Square was consumed by fire and many patients burned to death.
Visitors to Lisbon may still walk the ruins of the Carmo Convent, which were preserved to remind
Lisbonners of that day’s destruction.
2.5.2
Recent tsunamis in the Indian Ocean and Japan
Two recent tsunami events in the last decade have surely awaken people and governments to the reality
and unpredictability of these natural hazards, namely the 2004 Indian Ocean and the 2011 Tohoku
tsunamis. The 2004 Indian Ocean earthquake was an undersea megathrust earthquake
6
that occurred
on Sunday, 26 December 2004, with an epicenter off the west coast of Sumatra, Indonesia. The resulting
tsunami is given various names, including 2004 Indian Ocean tsunami or Boxing Day tsunami.
The earthquake was caused by subduction and triggered a series of devastating tsunamis along the coasts
of most landmasses bordering the Indian Ocean, killing over 230.000 people in fourteen countries, and
inundating coastal communities.
Asia
Hanimaadhoo
Male
Africa
Lamu
0˚
Colombo
Gan
Zanzibar
Pointe
La Rue
Australia
5 hr
40˚
Hillarys
St. Paul
Pointe La Rue
9 hr
Kerguelen
11 hr
Colombo
Cocos
Port Louis
7 hr
Cocos
Port Louis
M w =9.3
Diego Garcia
3 hr
20˚S
Hillarys
200 cm
Salalah
Sea level
20˚N
Hanimaadhoo
Salalah
Male
Lamu
Gan
Indian Ocean
13 hr
60˚
Mawson
Syowa Zhong Shan
0˚
20˚
40˚
60˚
Casey
Davis
Diego Garcia
Zanzibar
Dumont d'Urville
Antarctica
80˚
100˚
120˚
140˚E
26
27
28
26
27
28
Figure 2.7: Tide gauge measurements of the 2004 tsunami at several locations (Rabinovich & Thomson,
2007)
Figure 2.7 displays sea level measurements from several buoys along the Indian Ocean. As expected,
offshore amplitudes were relatively small. For instance, the tide gauge located at the Cocos Islands, 1700
[km] from the epicentre, revealed a small 30 [cm] high wave followed by a long sequence of oscillations,
6 A megathurst earthquake is an earthquake generated by seismicity at convergent plate boundaries, where the subducting
plate shallowly dips under the overriding plate. All recorded major earthquakes, with Mw over 9.0, are of the megathurst
type.
13
2. Tsunamis
Higher readings were provided by the gauges located at Gan and Diego Garcia Islands where amplitudes
were just below 1 [m]. As the tsunami approached continental platforms the wave height grew considerably. Gauge data yields amplitudes of nearly 3 [m] (Titov et al., 2005) at Pointe La Rue, Salalah and
Colombo. At Banda Aceh, Indonesia, water was reported to have risen by 3 [m] in under 1.5 minutes
(Chanson, 2005).
Important morphological changes were observed in Sri Lanka and Indonesia with some bed variations
reaching over 5 [m] (Chanson, 2005; Goto et al., 2011) at relevant harbors, for instance, the Lhoknga
harbour depicted in Figure 2.8, and at most seafront mobile reaches. This highlights the solid transport
capacity of the overland propagation stage, clearly visible on the aftermath pictures in Figure 2.8, and
the importance of effectively accounting for this feature in numerical simulations. Waste and debris are
mostly deposited inland but can also be carried offshore, when draw-down occurs, dispersing along the
shore.
(a) Lhokonga
before
the tsunami impact
(Chanson, 2005)
(b) Lhokonga after the
tsunami
impact
(Chanson, 2005)
(c) Phuket before the
tsunami
impact
(NASA)
(d) Phuket after the
tsunami
impact
(NASA)
Figure 2.8: Morphological impact of the 2004 Indian Ocean at Lhoknga, northwest coast of Sumatra, and
Phuhket, Thailand
This tsunami was one of the deadliest natural disasters in recorded history. Indonesia was the hardesthit country, followed by Sri Lanka, India, and Thailand. With a magnitude of Mw 9.1 – 9.3, it is the
third largest earthquake ever recorded on a seismograph. The plight of the affected people and countries
prompted a worldwide humanitarian response.
The 2011 earthquake off the Pacific coast of Tohoku and also known as the the Great East Japan
Earthquake or the 3/11 Earthquake, was a magnitude 9.0 (Mw) undersea megathrust earthquake off the
coast of Japan that occurred on Friday, 11 March 2011, with the epicenter approximately 70 kilometres
east of the Oshika Peninsula of Tohoku and the hypocenter at an underwater depth of 32 [km]. It was the
most powerful known earthquake ever to have hit Japan, and one of the five most powerful earthquakes
in the world since modern record-keeping began in 1900. The earthquake triggered powerful tsunami
waves that travelled as far as 10 [km] inland with up to 4 [m] high.
14
2.5. Tsunami consequences
Tohoku
Fukushima
Figure 2.9: Pressure gage readings for the 2011 Tohoku tsunami (Maeda et al., 2011)
Pressure gauge readings for this earthquake, such as the ones shown in Figure 2.9, show offshore heights
up to 1.2[m]. As previously referenced these small waves are prone to shoaling and the water height
that hits the shore is expected to be much higher. Several amateur footage was spread out on the media
showing flow depths ranging from 2 to 4 [m] high (Biggs & Sheldrick, 2011).
Significant morphological impacts occurred throughout all Japanese coast. Mobile reaches suffered severe
bed configuration changes as depicted in Figure 2.10. Built areas were completely wiped by the advancing
front of a turbulent, debris laden, flow.
(a) Natori before the tsunami impact
(b) Natori after the tsunami impact
(c) Fukushima before the tsunami impact
(d) Fukushima after the tsunami impact
Figure 2.10: Morphological impact of the 2011 tsunami at Natori and Fukushima, south of Tohoku (ABC
News, 2012)
15
2. Tsunamis
The Japanese National Police Agency confirmed 15,867 deaths. The tsunami also caused a number of
nuclear accidents, primarily the level 7 meltdowns at three reactors in the Fukushima Daiichi Nuclear
Power Plant. Many electrical generators were taken down, and at least three nuclear reactors suffered
explosions due to hydrogen gas that had built up after the cooling system failed.
Early estimates placed insured losses from the earthquake alone at US$14.5 to $34.6 billion.[31] The Bank
of Japan had to offer $15 trillion (US$183 billion) to the banking system on 14 March in an effort to
normalize market conditions. The World Bank’s estimated economic cost was US$235 billion, making it
the most expensive natural disaster in world history.
16
Chapter 3
Conceptual Model
3.1
Basic equations of fluid and sediment dynamics
The purpose of the present section is to introduce the fundamental equations of fluid dynamics within
the continuum hypothesis (Atkin & Crane, 1976). Under this hypothesis individual sediment grains are
ignored and the mass of sediment is uniformly distributed over the volume of the system. Furthermore,
this assumption is only valid at scales larger than the grain diameter and assumes homogeneity within
the sediment-fluid mixture. Along these lines, variables such as density, velocity or suspended solids
concentration are to be described as continuous time-dependent scalar or vector fields on R3 .
The physical phenomena that govern flow dynamics are to be mathematically represented by a set of
three laws, namely conservation of mass, momentum and energy. The derivation of these equations can
follow one of two distinct paths: a space-averaged estimation of flow effects over a finite region, usually
referred to as the integral form, or a more refined point-to-point description of the flow, the differential
form. The conceptual model, described in section 3, is significantly more consistent with the differential
form and so this will be the undertaken approach to derive the referred conservation equations.
3.1.1
Material derivative and transport theorem
The time-derivative of a given property φ(xi , t) of a fluid element within a velocity field ui (xi , t) can
be obtained by tracing this same element along its trajectory. In such case there will be two distinct
contributions for the total variation of φ: the first is due to the time dependence of the φ field, named the
local term, and the other is induced by the motion of the fluid element through spatial gradients of φ, the
advective term. The space coordinates xi are x1 = x, x2 = y and x3 = z, where x, y are the horizontal
coordinates and z is a vertical coordinate aligned with the gravity field. The time coordinate is t.
The material derivative, denoted by Dφ/Dt, is thus given by,
Dφ
∂
= φ(t, xi (t))
Dt
∂t
(3.1)
17
3. Conceptual Model
which in turn can be expanded, by means of the chain rule, to,
Dφ
∂φ
∂φ dxi
∂φ
∂φ
=
+
=
+ ui
Dt
∂t
∂xi dt
∂t
∂xi
(3.2)
where t stands for time and ui are the velocity vector components. This form clearly distinguishes the
two contributory terms for the total variation of φ: the local term, ∂t φ, and the advective term, ui ∂xi φ.
Another rather important concept in fluid mechanics is the Reynold’s Transport Theorem which states
that for any material volume V (t) and differentiable scalar field ϕ it is valid that,
d
dt
Z
Z
ϕ dV =
∂ϕ
∂
+
(ϕui ) dV
∂t
∂xi
(3.3)
V (t)
V (t)
which provides the necessary framework to formulate the conservation laws for mass, momentum and
energy (Wesseling, 2009).
3.1.2
Conservation of mass
This conservation law states that the rate of change of mass must equal the rate of mass production in
an arbitrary material volume V (t). Mathematically this is expressed by,
Z
d
dt
Z
ρ dV =
V (t)
σ dV
(3.4)
V (t)
where ρ is the apparent density of an homogeneous granular-fluid mixture and σ is the rate of mass
production per unit volume. The definition of ρ is as follows,
ρ = lim
δ∀→0
δM
δ∀
(3.5)
where δ an arbitrary difference operator, δM is the mass of the fluid-granular mixture and δ∀ is the
volume occupied by that mass. The mass of the mixture is M = Mw + Ms where M is the total mass
and Ms is the sediment mass and hence equation (3.5) becomes,
ρ = lim
δ∀→0
δMs δ∀s
δMw δ∀w
+
δ∀ δ∀w
δ∀ δ∀s
(3.6)
where ∀w and ∀w are the volumes of water and sediment, respectively. Note that ∀w = ∀ − ∀s and (3.6)
converges to ρ = ρw (1 − C) + ρs C. Defining the specific gravity as s = ρs /ρw on obtains the final form
of the mixture density,
ρ = ρw [1 + (s − 1)C]
18
(3.7)
3.1. Basic equations of fluid and sediment dynamics
By using the transport theorem, equation (3.3), one can expand equation (3.4) to,
Z
Z
∂ρ
∂
+
(ρui ) dV =
∂t
∂xi
V (t)
σ dV
(3.8)
V (t)
where ui is the velocity vector field and t stands for time. The mass production rate per unit volume, σ,
accounts for singular sinks or sources inside V (t). In this case there are no sinks or sources since all inputs
and outputs to the control volume can be expressed in terms of incoming or outgoing fluxes, respectively.
Provided that equation (3.8) holds for every arbitrary V (t), the differential form of the mass conservation
law can be obtained by extracting the left-hand integrand,
∂ρ
∂
+
(ρui ) = 0
∂t
∂xi
(3.9)
which is commonly known as the continuity equation as it requires no further assumptions other than
that both density and velocity are continuum functions (White, 2002).
The conceptual model is aimed at river flows which are greatly influenced by the transported sediment in
what concerns hydrodynamics. Since both water and sediment are incompressible so must be their mixed
flow but this alone does not imply a constant density value. Following the material derivative definition,
equation (3.2), one further obtains,
∂ui
Dρ
+ρ
=0
Dt
∂xi
3.1.3
(3.10)
Conservation of momentum
The momentum conservation equation is derived from Newton’s 2nd Law, the fundamental physical principle governing dynamics. Employing this law to balance the rate of change in linear momentum with all
the applied external forces on a fluid element will yield,
d
dt
Z
Z
ρui dV =
V (t)
Z
ρgi dV +
V (t)
fSij nj dS
(3.11)
∂V (t)
where ∂V (t) is the boundary surface of V (t), gi is the gravitational force per unit mass, nj is the outward
unit vector and fSij are the surface forces per unit area. Expanding the left term according to the
Reynold’s Transport Theorem will further demonstrate the implications of this conservation equation,
Z
V (t)
∂ρui
∂ρui
+ uj
dV =
∂t
∂xj
Z
Z
ρgi dV +
V (t)
fSij nj dS
(3.12)
∂V (t)
where ui stands for the velocity component being conserved.
19
3. Conceptual Model
The source terms on the right hand side are gravity (body force), pressure and viscous stresses (surface
forces). It can be shown that the surface forces are defined by the following 2nd order tensor, (Aris, 1989)
fSij = −p δij + τij
(3.13)
where p is the hydrostatic pressure, δij is the Kronecker delta and τij is the deviatoric stress tensor
(viscous stress tensor, in the case of pure water). By applying the Divergence Theorem to the surface
integral in (3.12) one obtains,
∂fSij
dV
∂xj
Z
Z
fSij · nj dS =
(3.14)
V (t)
∂V (t)
which, in further combination with equation (3.12), must give,
Z
Z
∂ρui
∂
+ ui
(ρuj ) dV =
∂t
∂xj
V (t)
ρgi +
∂
fS dV
∂xj ij
(3.15)
V (t)
which, being valid for any material volume V (t), can be condensed into a differential form by extracting
the integrands. Expanding the surface forces according to equation (3.13) will further yield,
∂ρui uj
∂ρui
∂p
∂τij
+
=−
+ ρgi +
∂t
∂xj
∂xi
∂xj
(3.16)
which completes the momentum conservation law. For the special case of incompressible flows, as the
null divergence hypothesis still holds, (3.16) can be further simplified to,
∂ρui
∂ρui
∂p
∂τij
+ uj
=−
+ ρgi +
∂t
∂xj
∂xi
∂xj
(3.17)
The Navier-Stokes equations
In order to complete the system of equations it is still necessary to relate the deviatoric stress tensor
τij with the velocity field by means of a constitutive relation. For clear water, the simplest of these
constitutive relations establishes that the viscous stresses are proportional to the strain rates and to the
dynamic viscosity, µ, (newtonian fluid) and is then given by, (Batchelor, 1967)
τij = µ
∂uj
∂ui
+
∂xj
∂xi
2
− ∇ · ~u δij
3
(3.18)
which, for incompressible flows, simplifies to,
τij = µ
20
∂ui
∂uj
+
∂xj
∂xi
(3.19)
3.2. Governing equations for stratified flows
Combining the momentum conservation equation for incompressible flows, equation (3.17), with the
previous consititutive relation (3.19) will yield the widely known Navier-Stokes equations.
∂ρui
∂ρui
∂p
∂
+ uj
=−
+ ρgi + µ
∂t
∂xj
∂xj
∂xj
∂ui
∂uj
+
∂xj
∂xi
(3.20)
which describe the flow in case of absence of suspended sediments.
3.2
Governing equations for stratified flows
The conceptual model underlying the present dissertation idealizes river flow as a stratified planar flow
of both water and sediment. Holding the hypothesis that the layers composing this stratified flow are
continua, the conservation equations suitable for describing such media are also applicable to each layer
and ultimately to the whole layer set. Each of these layers is a mixture of sediment and water characterized
by its respective sediment concentration and ensuing density.
As the model is designed in 2D further work is required on the equations presented in Chapter 3.1.
Specifically these equations must be brought from their generic 3D arrangement down to a 2D nearly
equivalent form. A consistent path to performing such transition is depth-averaging those conservation
equations. For a given function, for instance f (x, y, z, t), the depth-averaging is defined as
FL (x, y, t) =
1
hL
Z
y0 +hL
f (x, y, z, t) dz
(3.21)
y0
where FL (x, y, t) is the depth-averaged f (x, y, z, t) function within a layer L with depth hL .
It can be shown that integrating any function f (x, y, z, t) over a finite domain of the dependant variable
z will result in some function of the form Fz (x, y, t) which may not depend on z itself. Hence,
∂FL
= −1
∂z
(3.22)
Neglecting the vertical velocity (2D) in the conservation equations can be acceptable if the studied
phenomena have an horizontal scale that happens to be much greater than the vertical scale, which is
consistent with the hydrodynamics of a tsunami wave. Although the propagation of this type of wave
occurs within the deep waters of the oceans, its wave length is great enough so that the tsunami will
always interact with the seabed, thus generating friction, leading to a shallow-water wave behaviour.
The applicability of this conceptual model is restricted to flows where there is a clear distinction between
suspended sediment, transported within the water column, and bedload, composed by coarser sediment
transported near the riverbed Ferreira (2005). The layers featured in this conceptual model are those
depicted in Figure 3.1. The first layer, from herein onwards Layer [1] or clear water layer, is the upper
water column in which fine suspended sediments are transported. Underneath lies the contact load layer,
21
3. Conceptual Model
Layer [2] or bedload layer, which has a high concentration of coarse and suspended sediment. One must
note that the thickness of Layer [2] is expected to be much smaller than the one of Layer [1], and so the
set-up presented in Figure 3.1 only serves illustrative purposes. The bed, or Layer [b], was modeled as a
passive, non-moving repository of sediment that is allowed to undergo either deposition or scour.
Clearwater Layer
Layer [1]
Bedload Layer
Layer [2]
Bed
Layer [b]
Figure 3.1: Layer set featured in the conceptual model
As previously introduced in section 3.1 the governing equations are a set of partial differential equations
(PDEs) that describe fluid dynamics. The featured conservation equations consist of two mass conservation equations, for either total and solid phases, two momentum conservation equations, one for each
direction of the x − y plane and a set of closure equations. Since two distinct phases are featured, the
bedload layer and the clear water layer, the 2D conservation equations must be manipulated with special
care, so that all fluxes are accounted for and all quantities remain conserved.
Each layer is bounded by two interfaces. These are defined as continuum surface functions of the form,
SL (x, y, z, t) = IL (x, y, t) − z = 0
(3.23)
These surfaces can be defined as material or non-material. Material surfaces are not prone to fluxes which
implies that the fluid velocity vector at every point of these surfaces must equal the velocity vector of
the surface itself. On the contrary, non-material surfaces can be crossed and so no equality is imposed
between the velocity vectors of the fluid and surface (see Figure 3.2).
At some instant, t, the interface is defined by the surface z = St and the fluid particles lying on it are
defined by the surface z = St0 . After an infinitesimal time step, dt, both of these surfaces will change their
0
configuration into St+dt and St+dt
, respectively. The particles initially lying on the interface may now be
completely detached from it, meaning that the surface and these particles will have different velocities,
respectively ~uI and ~u. However these can be related to one another by an equivalent flux velocity, ~uφ ,
according to ~u = ~uI + ~uφ . The particular case where ~uφ = 0 defines a material surface, as it is the case
of the free-surface.
Whether a surface is defined as material or non-material it should always be a continuous surface and by
such definition one can use equation (3.2) and obtain,
SL (x, y, z, t) = 0 =⇒
22
DSL
=0
Dt
(3.24)
3.2. Governing equations for stratified flows
uF
u
uS
uF,t
nS
uF
S‘t+dt
uF,n=F
St+dt
uS
St=S’t
Figure 3.2: A non-material surface
that leads to the corresponding interface continuity equation,
DSL
∂IL
∂IL
∂IL
=
+ uI
+ vI
− wI = 0
Dt
∂t
∂x
∂y
(3.25)
If fluid velocity is the vector sum of the surface velocity and the equivalent flux velocity then u~I = ~u − u~φ
and one can rewrite equation (3.25) as,
∂IL
∂IL
∂IL
+ (u − uφ )
+ (v − vφ )
− (w − wφ ) = 0
∂t
∂x
∂y
(3.26)
∂IL
∂IL
∂IL
∂IL
∂IL
+ u
+v
− w + uφ
+ vφ
− wφ = 0
∂t
∂x
∂y
∂x
∂y
(3.27)
which will simply yield,
The flux across the interface, ΦIL , is defined as the projection of the equivalent flux velocity u~φ on the
interface normal since the tangential component of u~φ represents a displacement along the surface and
not through it, thus carrying no flux, as displayed in Figure 3.2. The surface normal is herein defined as
a unit vector collinear with the surface gradient.
ΦIL = u~φ · ~nIL = ~uφ ·
~ L
∇I
= ~uφ · r
~ Lk
k∇I
∂IL
∂x
∂IL
∂y
−1
2
∂IL 2
∂IL
+
+1
∂x
∂y
+
(3.28)
For the sake of simplicity an additional hypothesis is proposed: the interfaces are surfaces who happen
2
L
to be smooth enough so that ∂I
≈ 0 and so equation (3.28) can be written as,
∂x
~ L = uφ ∂IL + vφ ∂IL − wφ
ΦIL = u~φ · ∇I
∂x
∂y
(3.29)
Merging equations (3.29) and (3.27) and multiplying all the resulting terms by the fluid density at the
23
3. Conceptual Model
interface will yield
∂IL
∂IL
ρ|z=IL +
[ρui ]|z=IL − [ρw]|z=IL = ΦIL ρ|z=IL
∂t
∂xi
3.2.1
(3.30)
Conservation of mass
Subscript notation can be used to write a more compact form of the mass conservation equation for
incompressible flows derived in Chapter 3.1, (equation (3.10))
∂ρ ∂ρui
+
=0
∂t
∂xi
(3.31)
To integrate the terms of this equation over the depth of each layer it is necessary to make systematic
use of the Leibniz’ Integral Rule, which states that,
Z
b(β)
a(β)
∂
∂
f (x, β) dx =
∂β
∂β
Z
b(β)
f (x, β) dx +
a(β)
∂a(β)
∂b(β)
f (a(β), β) −
f (b(β), β)
∂β
∂β
(3.32)
and where all variables are generic. For presentation and readability purposes the depth-averaging will
be carried out separately for each term. The integration sequence will also follow the order in which these
terms appear in the mass conservation equation.
Recalling the depth-averaged function definition and the Leibniz Integral Rule, equations (3.21) and
(3.32), respectively, the depth-averaging of the local term yields,
Z
Lu
Ll
Z
∂ Lu
∂Ll
∂Lu
∂
ρ dz =
ρ dz +
ρ|z=Ll −
ρ|z=Lu
∂t
∂t Ll
∂t
∂t
∂
∂Ll
∂Lu
=
[ρL hL ] +
ρ|z=Ll −
ρ|z=Lu
∂t
∂t
∂t
(3.33)
where t is time, xi are the spatial coordinates (x, y, z) assuming that z is a vertical coordinate, ~ez = ~g /|~g |,
ρL is the depth-averaged density within a layer L enclosed by lower and upper interfaces, denoted by Ll
and Lu respectively. This procedure can be also used to integrate the convective term, thus obtaining,
Lu
Z
Ll
Z Lu
∂
∂Ll
∂Lu
∂
ρui dz =
ρui dz +
[ρui ]|z=Ll −
[ρui ]|z=Lu
∂xi
∂xi Ll
∂t
∂t
∂
∂Ll
∂Lu
=
[ρL UiL hL ] +
[ρui ]|z=Ll −
[ρui ]|z=Lu
∂xi
∂xi
∂xi
(3.34)
which is accurate for both xx and yy axis. The zz axis is a particular case of this integration,
Z
Lu
Ll
∂
∂
ρw dz =
∂z
∂z
Z
Lu
ρw dz +
Ll
∂Ll
∂Lu
[ρw]|z=Ll −
[ρw]|z=Lu
∂z
∂t
(3.35)
because the first term on the right hand side of equation (3.35) is obviously null. Following equation
24
3.2. Governing equations for stratified flows
(3.22) and provided that the interfaces are defined as surfaces of the form FL (x, y, z, t) = IL (x, y, t) − z,
equation (3.35) simplifies to,
Lu
Z
Ll
∂
ρw dz = [ρw]|z=Lu − [ρw]|z=Ll
∂z
(3.36)
concluding the depth-averaging of the continuity equation. Finally, summing equations (3.33), (3.32) and
(3.36) will yield,
∂
∂
∂Ll
∂Ll
[ρL hL ] +
[ρL UiL hL ] +
ρ|z=Ll +
[ρui ]|z=Ll − [ρw]|z=Ll
∂t
∂xi
∂t
∂xi
∂Lu
∂Lu
−
ρ|z=Lu
[ρui ]|z=Lu − [ρw]|z=Lu = 0
∂t
∂xi
(3.37)
which is solely an extended form of the mass conservation equation when integrated over the depth of a
given layer L,
Z
Lu
Ll
∂ρ ∂ρu ∂ρv ∂ρw
+
+
+
dz = 0
∂t
∂x
∂y
∂z
Equation (3.37) must still undergo further manipulation to clearly reveal the vertical fluxes through the
interfaces of overlying layers. Most of the interfaces featured in this model are non-material, being the
free-surface, the topmost surface, the only exception. Finally, by recalling equation (3.37)
∂
∂
∂Ll
∂Ll
[ρL hL ] +
ρ|z=Ll +
[ρL UiL hL ] +
[ρui ]|z=Ll − [ρw]|z=Ll
∂t
∂xi
∂t
∂xi
∂Lu
∂Lu
ρ|z=Lu +
[ρui ]|z=Lu − [ρw]|z=Lu = 0
−
∂t
∂xi
one can notice that the terms inside the brackets are in fact the vertical fluxes, as defined in (3.30).
Performing such substitution will yield the mass conservation equation for the conceptual model in its
final differential form,
∂
∂
∂
[ρL hL ] +
[ρL UL hL ] +
[ρL VL hL ] = ΦLu ρ|z=Lu − ΦLl ρ|z=Ll
∂t
∂x
∂y
(3.38)
where ΦLu ρ|z=Lu and ΦLu ρ|z=Ll are the net mass fluxes between Layer L and the layer above, through
interface Lu , and also with the layer underneath, through interface Ll .
Another particular feature of this formulation is that every layer that happens to have an appreciable
concentration of suspended sediment justifies two distint mass conservation equation: one for the solid
phase - suspended sediment - and another for the liquid phase - clean water. Recalling the definition of
mixture density, introduced in (3.5),
ρL = ρw (1 + CL (s − 1))
(3.39)
25
3. Conceptual Model
and bearing in mind that it evidences two distinct phases,
ρL = ρw (1 − CL ) + ρs CL
(3.40)
substituting (3.40) in (3.38) will result in a single equation that implicitly contains both solid and water
mass conservation equations. For the water mass conservation equation one has,
∂
∂
∂
[(1 − CL )ρw hL ] +
[(1 − CL )ρw UL hL ] +
[(1 − CL )ρw VL hL ] =
∂t
∂x
∂y
= (1 − CL ) ΦLu ρ|z=Lu − (1 − CL ) ΦLl ρ|z=Ll
(3.41)
while the solid mass conservation equation is simply written as the supplementary part,
∂
∂
∂
[CL ρs hL ] +
[CL ρs UL hL ] +
[CL ρs VL hL ] = CL ΦLu ρ|z=Lu − CL ΦLl ρ|z=Ll
∂t
∂x
∂y
(3.42)
Applying equations (3.41) and (3.42) to the layers featured in the conceptual model will complete the
final set of mass conservation equations.
For the suspended sediment layer, Layer [1] in Figure 3.1, one can assume that the sediment concentration
is negligible and therefore the depth-averaged layer density is roughly the same as clean water density
and so only equation (3.41) is relevant. If one considers the water density to be constant in both time
and space then it can be removed from the conservation equation,
∂
∂
∂
h1 +
[U1 h1 ] +
[V1 h1 ] = − Φ1,2
∂t
∂x
∂y
(3.43)
where U and V are the depth-averaged velocities in the ~ex and ~ey directions, respecively, and I1,2 is
the interface between Layer [1] (clear water) and Layer [2] (bedload). High concentration of sediment is
expected in Layer [2] which leads to using both (3.38) and (3.42),
∂
∂
∂
[h2 ] +
[U2 h2 ] +
[V2 h2 ] = Φ1,2 − Φb
∂t
∂x
∂y
(3.44)
∂
∂
∂
[C2 h2 ] +
[C2 U2 h2 ] +
[C2 V2 h2 ] = C2 Φ1,2 − (1 − p) Φb
∂t
∂x
∂y
(3.45)
where b is the bed, p is its respective porosity, Φ1,2 is the mass net flux between layers [1] and [2], Ib is
the bed interface and Φb is the net flux between the bed and Layer [2]. Figure 3.3 can further clarify the
physical meaning of the previous variables.
Since the bed is assumed to be saturated and to have no velocity at all there is only need for a sediment
26
3.2. Governing equations for stratified flows
mass conservation equation in this layer. Hence the only relevant equation is (3.42)
(1 − p)
∂
Zb = (1 − p) Φb = D − E
∂t
(3.46)
where p is the bed porosity, Zb is the bed elevation and D and E are equivalent deposition and erosion
rates, respectively.
Φ1,2
I1,2
Ib
U1
C1 ≈ 0
h1
h2
Clearwater Layer
Layer [1]
I1,2
C2
Zb
U2
Φb
Ib
Bedload Layer
Layer [2]
Bed
Layer [b]
Figure 3.3: Layer set featured in the conceptual model
3.2.2
Conservation of momentum
In this section the momentum conservation equation derived in Chapter 3.1, equation (3.16), will undergo
the same depth-averaging technique as the mass conservation equation, in the previous section. Recalling
the momentum conservation equation,
∂ρui
∂ρui
∂p
∂τij
+ uj
=−
+ ρgi +
∂t
∂xj
∂xi
∂xj
one can use it to justify the hydrostatic pressure field. If the vertical velocity is negligible so must be the
vertical acceleration which leads to that both left-hand side of equation (3.16) and stresses τzj are null.
This will simply resume (3.16) to,
∂p
= ρg ⇔
∂z
Z
0
h
∂p
dz =
∂z
Z
h
ρg dz ⇔ p(h) = ρgh
(3.47)
0
while showing why the pressure field is assumed to be hydrostatic.
Again, to improve the readability and understanding of the following steps, all integrations will be
performed separately for each hand of the momentum conservation equation, and also a term at a time
whenever deemed necessary. One should also note that the horizontal components of gi are null. Applying
the Leiniz’ Rule to the integration of the left-hand member over the depth of a given layer L will yield,
27
3. Conceptual Model
Z
Z
Lu
Ll
Lu
Ll
∂ρui
∂
∂Ll
∂Lu
dz =
[ρL UiL hL ] +
[ρui ] |z=Ll +
[ρui ] |z=Lu
∂t
∂t
∂t
∂t
(3.48)
∂ρui
∂
∂Ll
∂Lu
dz =
[ρL UiL hL ] +
[ρui ] |z=Ll +
[ρui ] |z=Lu
∂xj
∂xj
∂xj
∂xj
(3.49)
where the i index now stands for the 2D horizontal space directions - i = [1, 2]. One should note the
existing similarities between these two equations and equations (3.33) and (3.34). The steps used to
derive the vertical mass fluxes can also be applied to full extent on equations (3.48) and (3.49). After
those manipulations one obtains,
∂Lu
∂Lu
[ρui ] |z=Lu +
[ρui uj ] |z=Lu − [ρui w] |z=Lu = ΦLu [ρui ] |z=Lu
∂t
∂xj
∂Ll
∂Ll
[ρui ] |z=Lu +
[ρui uj ] |z=Ll − [ρui w] |z=Ll = ΦLl [ρui ] |z=Ll
∂t
∂xj
(3.50)
(3.51)
which represents the vertical net flux of momenta, in the ~ei direction, through the interfaces of some layer
L. The depth-averaged left-hand member of (3.16) can be therefore written as
∂
∂ [ρL UiL hL ] +
ρL UiL ULj hL + ΦLl [ρui ] |z=Ll − ΦLu [ρui ] |z=Lu
∂t
∂xj
(3.52)
For the right-hand member of (3.16), the applied external forces, the same logic must be followed in order
to achieve a consistent set of depth-averaged quantities.
The first term to be integrated is the pressure field, which was shown to be hydrostatic,
Lu
Z
Ll
∂p
∂
dz =
∂xi
∂xi
Z
Lu
p(z) dz +
Ll
∂Ll
∂Lu
p|z=Ll −
p|z=Lu
∂xi
∂xi
(3.53)
The depth-averaged pressure is obtained by using the variable substitution ξ = z−Ll while also accounting
for for the weight of all the nL layers overlying layer L,
Lu
Z
Z
p(z) dz =
Ll
hL
p(ξ) dξ =
0
nL
X
1
γL h2L +
γLl hi hL = PL
2
(3.54)
i=L+1
and γL is the same as ρL g. The depth-averaging of the viscous stresses is performed by integrating the
deviatoric stress tensor over the depth of a layer, τij , yielding,
Z
Lu
Ll
Z Lu
∂
∂Ll
∂τij
∂Lu
dz =
τij dz +
τij |z=Ll −
τij |z=Lu
∂xj
∂xj Ll
∂xj
∂xj
∂Ll
∂Lu
∂
=
TLij hL +
τij |z=Ll −
τij |z=Lu
∂xj
∂xj
∂xj
(3.55)
Adding equations (3.52) to (3.55) will yield a 2D depth-averaged equation of the momentum conservation
28
3.3. Closure equations
equation that, even while apparently more complicated, is in fact simpler to solve as it involves fewer
unknowns, due to the removal of the vertical velocity, and thus reduces the governing system to three
conservation laws. So, the 2D version of equation (3.16) is given by,
∂
[ρL UiL hL ] +
∂t
∂Ll
−
p|z=Ll +
∂xi
∂ ∂
ρL UiL ULj hL + ΦLl [ρui ] |z=Ll − ΦLu [ρui ] |z=Lu = −
PL
∂xj
∂xi
∂Lu
∂
∂Ll
∂Lu
p|z=Lu +
TLij hL +
τij |z=Ll −
τij |z=Lu
∂xi
∂xj
∂xj
∂xj
(3.56)
Applying equation (3.56) to Layer [1] will yield the following result,
∂ ∂
[ρ1 U1i h1 ] +
ρ1 U1i U1j h1 + ΦI1,2 [ρui ] |z=I1,2 =
∂t
∂xj
∂
∂I1,2
∂
∂I1,2
=−
P1 −
p|z=I1,2 +
T1 hL +
τij |z=I1,2
∂xi
∂xi
∂xj ij
∂xj
(3.57)
and similarly for Layer [2] one should obtain,
∂ ∂Zb
∂
∂
[ρ2 U2i h2 ] +
ρ2 U2i U2j h2 − Φ1,2 [ρui ] |z=I1,2 = −
P2 −
p|z=Zb
∂t
∂xj
∂xi
∂xi
∂
∂Zb
∂I1,2
∂I1,2
p|z=I1,2 +
T2ij h2 +
τb,j −
τij |z=I1,2
+
∂xi
∂xj
∂xj
∂xj
(3.58)
A total momentum conservation equation, accounting for both layers, is obtained if one sums equations
(3.57) and (3.58). Considering the total values of momentum, momentum flux and applied external forces
for both layers as a whole one obtains,
∂
∂
∂Zb
∂
∂Zb
∂
[ρUi h]T +
[ρUi Uj h]T = −
PT −
ghT +
Tij T hT +
τb,j
∂t
∂xj
∂xi
∂xi
∂xj
∂xj
(3.59)
where all the variables with index T represent the sum of their equivalents over both layers 1 and
P2
2: generically ϕT =
i=1 ϕi , where ϕ can be ρUi h, ρUi Uj h, P , h or Tij hi . Due to the hydrostatic
pressure field it is also equivalent to write ∂xi Zb ghT instead of ∂xi Zb p|z=Zb . Provided that there are no
discontinuities in the velocity field across the interface that bounds both layers 1 and 2, Ferreira et al.
(2009) showed that the vertical fluxes of momentum Φ1,2 [ρui ] |z=I1,2 and the terms ∂xi I1,2 p|z=I1,2 are
cancelled when (3.57) and (3.58) are summed.
The governing system of equations has thus been formally derived and its complete definition depends
only on adequate closures or closure equations, which will be derived in the next section.
3.3
Closure equations
The proposed system of governing equations, namely (3.43), (3.44), (3.45), (3.57) and (3.58), requires
closures for the dynamics of the bedload layer, h2 and U2 , for the evolution of bed morphology, Φb or
D − E, and finally for flow resistance, τb,j and TLij (Canelas, 2010; Ferreira et al., 2009).
29
3. Conceptual Model
The closures for the dynamics of the bedload layer were derived under the premise that greater
kinetic energy fluxes will lead to an increase of the bedload layer thickness, thus allowing for a complete
dissipation of the kinetic energy fluctuating components (Ferreira, 2005).
The proposed closure for the bedload layer thickness, h2 , is given by, (Ferreira, 2005)
h2
= m1 + m2 θ
ds
(3.60)
where ds is the significant particle diameter, m1 and m2 are parameters that are dependent upon the
mechanical properties of the particles and fluid viscosity and lastly θ is the Shield’s parameter, calculated
as θ = Cf k~uk2 / (gds (s − 1)) where Cf is the friction coefficient. The bedload layer depth-averaged velocity,
U2 , is given by a power law derived by Ferreira (2005),
U2 = U
h2
h
1/6
(3.61)
Bed morphology closures relate the existing imbalance between solid discharge and transport capacity
with bed elevation changes. The solid discharge, qs , is given by qs = CρU while several formulas can be
used for the transport capacity, qs∗ , namely Meyer-Peter & Muller (1947), Bagnold (1966) and a specific
formula for stratified flows and debris flows by Ferreira (2005). The closure equation for ∂t Zb is thus
given by,
qs − qs∗
∂Zb
=
(1 − p)−1
∂t
Λ
(3.62)
where p is the bed porosity and Λ is an adaptation length defined as a two horizontal asymptotic curve
Λ(Λmin , Λmax , θref , θ) whose arguments Λmin , Λmax and θref are parameters used to ensure that the
value of Λ will never grow too large nor will it decrease to nearly zero.
Closures for flow resistance feature bottom friction, viscous stresses and turbulent stresses. The
bedload layer is expected to carry a high concentration of sediment leading to the same rheological
behavior of an hyperconcentrated flow.
In such conditions the total shear stress is the sum of three main components: the yield stress, τy , the
viscous stresses, τν , and the turbulent stresses, τt . The cohesive nature of very fine particles is taken into
account through constant values of the yield stress, since most of the presently available formulae are
only supported by very singular experiments.
A useful property of a Newtonian fluid is that the viscous stresses can be expressed in terms of both
dynamic viscosity, µ, and vertical velocity gradients, ∂z ui . As for the turbulent stresses induced by
bottom interaction they can be associated to the second power of the vertical velocity gradient as these
30
3.3. Closure equations
are related to the kinetic energy of the flow (Takahashi, 2007). Total shear stress is thus given by
τi = τy + τν + τt = τy + µ
∂ui
+ ρ d2s cf
∂z
∂ui
∂z
(3.63)
The bed shear stress is given by,
τbi = Cf ρk~ukui
(3.64)
where the friction coeficient Cf can be given by the Manning-Strickler formula,
1
Cf =
gh 3
Ks2
(3.65)
or by a specific formula concerning stratified flows or debris flows (Ferreira et al., 2009),
k~ukds
hωs
Cf =
(3.66)
where ωs is the settling velocity of the sediment particles (Jimenez & Madsen, 2003) and Ks is the
Manning coefficient.
For the turbulent stresses a simple κ − model is applied (Pope, 2000). This model ignores small-scale
eddies in the motion and calculates a large-scale motion with an eddy viscosity that carries and dissipates
the energy of the smaller-scale flow.
Under these assumptions one can define the depth-averaged stress tensor Tij as,
Tij = ρνT
∂uj
∂ui
+
∂xj
∂xi
(3.67)
where νT is the eddy viscosity defined in terms of the depth-averaged kinetic energy, κ, and turbulent
kinetic energy dissipation rate, ,
4.375u2∗
κ2
= 0.09
νT = Cµ
5u3∗ /h
2
where u∗ is the shear velocity or friction velocity, calculated as u∗ =
(3.68)
√
τb /ρ.
The system of equations is now completely defined featuring governing equations (3.43), (3.44), (3.45),
(3.57), (3.58) and closure equations (3.63) to (3.68).
31
3. Conceptual Model
32
Chapter 4
Numerical Model
4.1
Discretization scheme
4.1.1
Introduction
Eulerian mesh-based methods featuring finite difference, finite element or finite volume schemes are
amongst the most widely used tools in computational fluid dynamics (CFD) (Hircsh, 1988). The employed
numerical model (Canelas, 2010; Murillo & García-Navarro, 2010) stands on a fully conservative finite
volume method (FVM), that has now became a standard within the CFD world. Its important property
of providing a proper discrete representation of the conservation laws when discontinuities are present in
the domain (Hircsh, 1988; LeVeque, 1990; Toro, 2009) is a valuable advantage over other methods.
Only one flow layer is currently implemented in the computational model and hence the governing equations derived in Chapter 3 can be simplified. The absence of adjacent flow layers vanishes all the vertical
flux terms obtained when depth-averaging the governing equations with Leibniz’s Integral Rule, as demonstrated for the total momentum conservation equation (3.59) in chapter 3. The only exceptions are the
total and sediment mass fluxes between the flow layer and the bed. The governing system features total
mass conservation, equation (4.1), total momentum conservation for the x and y directions, equations
(4.2) and (4.3), respectively, and sediment mass conservation, (4.4).
∂t h + ∂x (hu) + ∂y (hv) = −∂t Zb
1
∂x hTxx −
∂t (uh) + ∂x u2 h + 12 gh2 + ∂y (uvh) = −gh∂x Zb −
ρm
1
∂t (vh) + ∂x (vuh) + ∂y v 2 h + 12 gh2 = −gh∂y Zb −
∂x hTyx −
ρm
(4.1)
1
τb,x
∂y hTxy −
+ Sx
ρm
ρm
1
τb,y
∂y hTyy −
+ Sy
ρm
ρm
∂t (Cm h) + ∂x (Cm hu) + ∂y (Cm hv) = −(1 − p)∂t Zb
(4.2)
(4.3)
(4.4)
where x, y are the space coordinates, t is time, h is the depth, u and v are the depth-averaged velocities
in the x and y directions, respectively, Zb is the bed elevation, ρm and Cm represent the depth-averaged
33
4. Numerical Model
density and concentration of the mixture, respectively, Tij are the depth-averaged turbulent stresses
and τb represents the friction exerted by the bed on the fluid. Sx and Sy are source terms induced
by 2nd order perturbations due to stratification1 (Ferreira, 2005). The bed variation ∂t Zb in given by
(1−p)∂t Zb = (qs −qs∗ )/Λ where qs is the sediment discharge, qs∗ is its capacity value and Λ is an adaptation
length. Sediment inertia and added pressure are neglected in the momentum equations.
Bed slope and friction source terms are of special importance for well-balanced CFD models. The employed formulation (Murillo & García-Navarro, 2010) represents an effort to correctly introduce these
terms while ensuring compatibility with steady solutions 2 , a feature lacked by several other models.
4.1.2
Finite volume schemes: The case of the 2D shallow-water equations
FVM’s discretize the physical domain into control (finite) volumes where the conservation laws are applied
using an Eulerian perspective of the flow: fluxes and accumulated quantities must be balanced for each
finite volume according to the system of equations derived in Chapter 3. One should note that the
numerical model herein presented uses a different notation from the Conceptual Model presented in
chapter 3: vectors are represented by bold characters, for instance the conservative variables vector, U ,
the i index refers to cells and the k index refers to cell edges.
Considering some quantity U : R × R → Rm , where m is the system dimension, defined in an arbitrary
volume Ω, fixed in both space and time, bounded by a closed set of surfaces Γ then the local value of U
is dependent on incoming and outgoing fluxes, denoted by E ∈ Rm×m , expressing the contributions of
the neighboring volumes. Singular sources also contribute and are herein denoted by H : R × R → Rm .
The fluxes are functions of the quantity in question, which is to say E = E(U ). A conservation law for
the quantity U can be integrated over the the domain Ω yielding,
Z
∂t
I
∇·E =
U dΩ +
Ω
Z
Γ
HdΩ
(4.5)
Ω
Expanding this equation for a fixed control volume, Ωi , a discrete form of the conservation law is obtained,
Ai ∂tU i +
N
X
(E · n)Lk = Ai H i
(4.6)
k=1
where U i is the dependant variable vector, H i is the source term vector, both defined as cell-averages, Ai
is the area of cell Ωi which is bounded by N permeable boundaries. The system composed by equations
(4.1) to (4.4) can be written in a compact form as,
∂t U (V ) + ∇ · E(U ) = H(U ) ⇔ ∂t U (V ) + ∂x F (U ) + ∂y G(U ) = H(U )
(4.7)
where the dependant non-conservative variables, V , dependant conservative variables U and fluxes F
1 For
further details refer to Appendix 3, page 107
as the case of still water in hydrostatic equilibrium over steep bottom topography.
2 Such
34
4.1. Discretization scheme
and G are defined as,

h


h




uh
vh

 2


 



u h + 21 gh2 
 uh 
 u 


vuh








; F =
V =  ; U = 
; G= 2


1
2
uvh

 vh 
 v 
v h + 2 gh 

Cm h
Cm
Cm hv
Cm hu
(4.8)
and the source term vector H is H = R+T+S, where R is the parcel of source terms not susceptible to be
treated as non-conservative fluxes, T are non-physical fluxes and S are stratification induced source terms.
These terms are treated as fluxes in order to obtain a well-balanced numerical scheme (Vázquez-Cedón,
1999) and are defined as,

− (D−E)
1−p



0




b,x
−gh∂x Zb 
 − ρτ(w)





R=
τb,y
 ; T = −gh∂ Z  .

 − ρ(w) 
y b
0
−(D − E))
(4.9)
To obtain the FVM discretization, the system (4.7) is integrated over a cell i and the Gauss theorem is
applied to the resulting divergence term. This procedure leads to,
Z
∂t
I
E (U ) · ndl =
U (V )dS +
Ωi
Z
Γi
H(U )dS
(4.10)
Ωi
where E · n = F nx + Gny with n = (nx , ny )T being the outward normal to cell Ωi . Denoting the cell
area by Ai and given that all quantities are piecewise constant fields, one obtains,
I
∂t Ai hU i i +
E (U ) · ndl = Ai hH i i
(4.11)
Γi
where the h i represents the spatial average in the cell, given by
1
hξ i i =
Ai
Z
ξ (x, y, t) dS.
(4.12)
Ωi
Performing the surface integrals on the N boundaries of cell i, equation (4.11) becomes
∂t Ai hU i i +
N
X
Lk hE · niik = Ai hH i i
(4.13)
k=1
where Lk is the length of the k edge. The local discretization of the system becomes,
n
Ai
i
∆ hU i i X
+
Lk ∆ik h(E − T ) · ni = Ai hRi i
∆t
(4.14)
k=1
35
4. Numerical Model
where the flux variation across the k th] edge of cell i, ∆ik hE − T i, is expressed as a function of the
local variation of the dependent conservative variables, using Roe’s approximate Riemann solvers (Roe,
1981). Note that the flux vector of a system of shallow-flow conservation laws is not homogeneous and,
hence, it is not possible to perform an exact flux vector splitting. Assuming a local linearization of the
conservative flux vector orthogonal to the edge in question, one obtains (Canelas, 2010; Toro, 2001),
∆ik hE · ni = (hE j i − hE i i) · nik =
NW
X
(n) (n) (n)
λ̃ik αik ẽik
(4.15)
n=1
(n)
(n)
where NW is the dimension of the eigenspace, λ̃ik and ẽik are the approximate eigenvalues and eigenvectors, respectively, and are defined as,
(1)
(2)
(3)
(4)
λ̃ik = (ũ · n − c̃)ik ; λ̃ik = (ũ · n)ik ; λ̃ik = (ũ · n + c̃)ik ; λ̃ik = (ũ · n)ik

(1)
ẽik
1


0


1







−c̃ny 
ũ + c̃nx 
ũ − c̃nx 
(2)
(3)
 ; ẽ(4)





; ẽik = 
; ẽik = 
=
ik



 c̃nx 
 ṽ + c̃ny 
 ṽ − c̃ny 
C˜m
0
C˜m
ik
ik
ik
 
0
 
0

=
0
 
1
(4.16)
(4.17)
ik
(n)
where coefficients αik are the wave-strengths for each of the eigenvectors and ũ, ṽ, c̃ and C̃m are
the approximate values for the two velocity components, (u, v), free-surface disturbance celerity, c, and
(n)
sediment concentration, Cm , respectively. The eigenvector wave-strengths αik are given by,
(1)
αik =
∆ik hhi
2
(2)
αik =
(3)
αik =
−
1
c̃ik
∆ik hhi
2
+
1
2c̃ik
(∆ik huhi − ũik ∆ik hhi) · nik
(∆ik huhi − ũik ∆ik hhi) · tik
1
2c̃ik
(4.18)
(∆ik huhi − ũik ∆ik hhi) · nik
(4)
αik = ∆ik hhCm i − C˜m ∆ik hhi
and the approximate values of u, v, c and Cm are defined as, (Toro, 2001)
ũik
r
p
p
p
√
√
√
vi hi + vj hj
Cmi hi + Cmj hj
ui hi + uj hj
hi + hj
p
p
p
; ṽik = √
; c̃ik = g
; C̃mik =
= √
√
2
hi + hj
hi + hj
hi + hj
(4.19)
Murillo & García-Navarro (2010) expressed the variation of the non-conservative flux vector as a linear
expansion in the basis formed by the approximate eigenvalues,
∆ik hT · ni =
NW
X
(n) (n) (n)
λ̃ik βik ẽik
(4.20)
n=1
(n)
where βik are wave-strengths corresponding to the eigenvalues and pb is the pressure exerted by the
bottom boundary on the flow over a bottom step. Expressions for pb (Canelas, 2010) that guarantee a
well-balanced scheme and energy dissipation at jumps have been derived by Murillo & García-Navarro
36
4.1. Discretization scheme
(2010). The wave-strengths associated with the eigenvalues are given by,
(1)
βik = −
1
2c̃ik
pb
ρw
(2)
(3)
(1)
(4)
; βik = 0; βik = −βik ; βik = 0
(4.21)
The remaining source term, T , is discretized in a point-wise formulation. T n+1
is computed in an
i
intermediate step with the updated variables of the homogeneous part of equation (4.13). The shear-rate
in the definition of the depth-averaged turbulent stresses is calculated from directional derivatives. The
discretized directional derivative ∇ab F , of a generic differentiable variable F , in the direction of the vector
r a b that connects the barycentres of cells a and b is,
Dab F =
F (xb , yb ) − F (xa , ya )
|r i,j |
(4.22)
At each cell, three directional derivatives can be defined, for each of the neighboring cells. In a Cartesian
referential the directional derivatives are written as ∇ab F = ∇F erab , with ∇ = ∂x + partialy and erab
being the directional unit vector from the barycentre of cell a to the barycentre of cell b. With the
unknowns as the gradients in the Cartesian referential, a pair of directional derivatives can be used to
calculate ∂x F and ∂y F . In general
Dab
Dac
!
=
cos(ηab )
sin(ηab )
cos(ηac )
sin(ηac )
!
∂x F
!
∂y F
(4.23)
where ηab and ηac are the angles that the segments linking the barycentres of cells a and b and a and c,
respectively, make with the x direction. Inverting the system allows for the calculation of ∂x F and ∂y F .
The eigenspace dimension for this scheme is NW = 4 and since the cells are of triangular shape then
N = 3. Introducing equations (4.15) and (4.20) into (4.14) and following the Godunov’s method (LeVeque,
2002), the volume integral in the cell at time tn+1 is given by
U n+1
= U ni −
i
3
3
4 −
X
∆t X
(n)
Lk
λ̃(n) α(n) − β (n) ẽik + ∆t Tin+1
Ai
ik
n=1
(4.24)
k=1
which comprises the desired flux-based finite volume scheme.
4.1.3
Stability region
When cell-averaging the solution, the time step is chosen small enough to guarantee that there is no
interaction between waves obtained as the solution of the Riemann Problem (RP) at adjacent cells. The
stability region considering the homogeneous part of the system becomes
3 Further
details on source term treatment and flux-vector splitting techniques can be found in (Canelas, 2010)
37
4. Numerical Model
∆t ≤ CFL ∆tλ̃
∆tλ̃ =
min (χi , χj )
max |λ̃(n) |n=1,2,3
(4.25)
where the CFL value is always lower than 1 and χ is an equivalent 1D distance defined as,
χi =
Ai
max (Lk )i,k=1,2,3
(4.26)
Details on time step restrictions due to source terms can be found in Murillo & García-Navarro (2010)
4.1.4
Wetting-drying algorithm
In wet/dry interfaces with discontinuous bed level, the time step can become practically null since negative
fluid depths can be expected enforcing positivity requires a very small time step. In order to ensure
positivity and conservation in the solution for all cases, the fluxes for the update of the conserved variables
are computed as (Murillo & García-Navarro, 2010),
If hnj = 0 and h∗∗∗
<0
j
−
then (∆E − T )−
ik = (∆E − T )k and (∆E − T )jk = 0
If hni = 0 and h∗i < 0
−
then (∆E − T )−
jk = (∆E − T )k and (∆E − T )ik = 0
(4.27)
otherwise, the update is performed by (4.24). These restrictions efficiently prevent the appearance of
negative fluid depths: the flux is zero at cells where an intermediate step predicts negative flow depths.
In areas were h approaches zero, the computation of the primitive variables u = uh=h leads to potentially
large round-off errors. Instead of defining a threshold for h above which the cell is allowed momentum,
a threshold for the flow depth (hmin ) is set to compute velocities using Kurganov & Petrova (2007)
u= p
√
uh 2h
h4 + max(h4 , )
(4.28)
√
where = hmin max(1, 2A), hmin is a mesh dependant parameter and A is the cell area.
4.1.5
Entropy correction
Non-physical solutions generated by Roe’s approximate Jacobian can occur, since there is no guarantee
that the entropy condition 4 . is satisfied. To enforce a physically admissible solution an entropy correction
must be employed. In cases were the flow is subcritical in one side and supercritical on the other side of
an edge, entropy violating stationary shocks are computed by the described FVM (Canelas, 2010).
After the local linearization of the system the RP’s left and right rarefaction waves are identified by
(1)
(1)
(3)
(3)
λik < 0 < λjk and λik < 0 < λjk , respectively, with the wave speeds being the eigenvalues computed
4 Detailed
38
in (LeVeque, 2002)
4.2. Boundary conditions
with the independent variables of the cells sharing edge k: λik = λ(U i ) and λjk = λ(U j ) as opposed to
the approximated λ̃(n) used in the update.
A version of the Harten-Hyman entropy fix (LeVeque, 2002) is implemented, were the objective is to
forcefully promote the propagation of information in any expansion wave by diffusing the shock. Murillo
(1)
(3)
& García-Navarro (2010) decomposed the initial jump associated to λ̃ik and λik in two jumps.
(n)
The same concept can be applied to the bottom source terms coefficients, βik , but only if a conservative
(n)
splitting is verified. A proportional splitting to the new waves associated to λ̃ik results in bad solutions (Murillo & García-Navarro, 2010) due to numerical instabilities, and the time step suffers drastic
reductions 5 .
4.2
4.2.1
Boundary conditions
Introduction
Ensuring the existence of unique solutions requires suitable conditions at the boundary ∂Ω of any domain
Ω. To accurately model physical phenomena the solution should depend continuously on the variable
values inside Ω and those imposed on ∂Ω (Wesseling, 2009). This leads to the concept of well-posed
problems, according to which, mathematical models of physical phenomena should provide unique solutions, whose behavior is smooth when small disturbances are introduced in either initial or boundary
conditions.
The boundary condition scheme firstly implemented in the numerical model relies on a Riemann’s Invariant solution. Whilst this is quite advantageous from the computational performance and mathematical
well-posedness perspectives it also carries a major drawback: if flow direction at an inlet is opposite to
the inner normal of ∂Ω then a unique solution does not exist.
In the interest of overcoming this limitation, new boundary conditions were developed and implemented.
The derived and implemented scheme is based on the Rankine-Hugoniot shock conditions.
4.2.2
Rankine-Hugoniot conditions
When the characteristics of partial differential equations (PDE) intersect, a shock is formed. These
are very thin surfaces propagating rapid transitions in physical variables, such as velocity or pressure,
and may be mathematically represented by discontinuities. The Rankine - Hugoniot relations or jump
conditions (Rankine, 1869) suitably describe the relationship between the states on both sides of the
shock wave. Considering the following hyperbolic conservation law for a quantity ϕ,
∂ϕ
∂
+
f (ϕ) = 0
∂t
∂x
5 Further
(4.29)
details can be found in Murillo & García-Navarro (2010)
39
4. Numerical Model
and should the solution exhibit a discontinuity at x = xs (t), integrating over the domain x ∈ [xR , xL ],
where xR < xs (t) < xL , will yield,
f (ϕ(xL , t)) − f (ϕ(xR , t)) =
d
dt
Z
xs (t)
ϕ(x, t) dx +
xL
d
dt
Z
xR
ϕ(x, t) dx
(4.30)
xs (t)
where the subscripts L and R denote the sides immediately to the left (upstream) and right (downstream)
of the shock, respectively. Using Leibniz’ Rule, equation (3.32), on the integral derivatives of (4.30) will
further obtain, (LeVeque, 1990; Toro, 2009)
f (ϕ(xL , t)) − f (ϕ(xR , t)) =
d xs
d xs
ϕ(xsL (t), t) −
ϕ(xsR (t), t) +
dt
dt
Z
xR
xL
∂ϕ
dx
dt
(4.31)
where ϕ(xsL (t), t) and ϕ(xsR (t), t) are the limits of ϕ(xs (t), t) as x converges to xs (t) from left or right,
respectively.
By forcing the x domain limits, xL and xR , to converge on the shock location xs (t) one obtains the
following result,
f (ϕ(xL , t)) − f (ϕ(xR , t)) = [ϕ(xsL (t), t) − ϕ(xsR (t), t)] S
(4.32)
where S is the speed of the shock, dxs /dt. Equation (4.32) is commonly referred to as the RankineHugoniot shock conditions and is often represented by a more compact form,
∆f = ∆u S
(4.33)
where the jumps are defined as ∆f = f (ϕ(xL , t)) − f (ϕ(xR , t)) and ∆ϕ = ϕ(xL , t) − ϕ(xR , t).
4.2.3
Shocks
Despite being discontinuities, shocks are still required to obey the conservation equations derived in
Chapters 3.1 and 3, namely the mass conservation equation, (3.10), and the momentum conservation
equation, (3.17). The conservation of mass across a shock is thus given by,
uL hL − uR hR = S (hL − hR )
(4.34)
where uk and hk with k = L, R are the velocities and water depths immediately to the left and right of
the shock wave while S is the shock speed.
The bottom elevation is assumed to have a zero-gradient type of boundary condition and thus, at the
boundary edges, the gravitational force is neglectable. Hence, the momentum conservation equation
40
4.2. Boundary conditions
across a shock at ∂Ω is,
1
1
u2L hL + h2L − u2R hR + h2R = S (uL hL − uR hR )
2
2
(4.35)
Expressing the shock speed, S, in terms of unit discharge and water depth jumps will yield,
S=
uL hL − uR hR
hL − hR
(4.36)
and, additionally, replacing (4.36) in (4.35) further obtains,
1 2
uL hL − uR hR
1 2
2
2
(uL hL − uR hR )
uL hL + hL − uR hR + hR =
2
2
hL − hR
(4.37)
which is a significantly more useful expression since all variables are velocities or water depths. Figure
4.1 illustrates a compressive shock at an edge of the boundary, ∂Ω.
UL
hL
qL
hR
UR
Figure 4.1: A shock (compressive) at the boundary, ∂Ω
4.2.4
Expansion waves
While the Rankine-Hugoniot conditions are suitable to describe compressive shocks, the provided solutions
are not physically feasible for expansive shocks (due to entropy condition violation). To handle such
shocks, for instance when the mass flux entering Ω is lower that the mass flux on the inner side of ∂Ω,
an alternative formulation, consisting of expansion waves, is proposed.
The propagation speed of these expansion waves is assumed to have the same value of the characteristic
λR . The expression itself is very close to the one of the Rankine-Hugoniot condition but, since the speed
of propagation is known, only one conservation equation is required. Since the momentum conservation
law features more information (velocity, water depth and pressure) it was the chosen equation to define
the expansion wave conditions,
1 2
1 2
2
2
uL hL + hL − uR hR + hR = λR (uL hL − uR hR )
2
2
(4.38)
41
4. Numerical Model
and, since λR = uR +
√
ghR , one can further obtain,
1
u2L hL + h2L
2
h
i
p
1
− u2R hR + h2R = uR + ghR (uL hL − uR hR )
2
(4.39)
which defines the expansion wave condition and where all terms have the same meaning of those in
equation (4.37).
4.2.5
Analytical solutions
Equations (4.37) and (4.39) must be expressed in terms of the primitive variables inside Ω, uR
6
and
hR , and the hydrographs imposed at the permeable boundaries. These hydrographs define the discharges
(Q) or the water depths (h) as time-dependent functions, Q(t) at the inlets and h(t) at the outlets.
Multiplying (4.37) by hL will define the compressible shocks at the inlets as,
1
1
2
u2L h2L + h3L − u2R hR hL + h2R hL (hR − hL ) = (uL hL − uR hR ) hL
2
2
(4.40)
where L and R denote variables upstream and downstream of ∂Ω, respectively. Expanding all terms will
yield,
1
1
1
1
2
qL
hR + h4L − h3L hR − u2R hR h2L + hR h2L − h3R hL = −2qL uR hR hL
2
2
2
2
(4.41)
which after further rearrangement will evidence a 4th order polynomial in order to hL ,
1
1
1
1
2
2 2
− ghR h4L + gh2R h3L + qR
+ gh3R h2L + gh4R − 2qL uR hR hL + qL
hR = 0
2
2
2
2
(4.42)
As for the expansion waves, equation (4.39), at the inlets, a similar procedure will obtain,
p
1
1 2
3
2
3
− g hL + (uR hR − qL ) uR + ghR − uR hR − ghR hL + qL
=0
2
2
(4.43)
where qL is the value of the discharge hydrograph, Q(t). Equations (4.42) and (4.43) will provide the
required solutions at inlet boundaries, hL and uL
7
, weather using compressive shocks or expansion
waves.
For the outlet boundaries these polynomials must be expressed in terms of uL , as hL is imposed by the
h(t) hydrograph, and so the compressive shocks are given by the following alternative arrangement of
6u
~ n∂Ω where ~
n∂Ω is the inner normal of
R is defined as the normal component of the velocity vector at ∂Ω, uR = U · ~
~ does not influence compressive shocks (nor expansion waves).
∂Ω. The tangential component of U
7 u is simply given by u = q /h
L
L
L
L
42
4.2. Boundary conditions
equation (4.37),
[−hL hR ] u2L + [2uR hR hL ] uL +
1 3
g hR + h3L − h2L hR + h2R hL − u2R hR hL = 0
2
(4.44)
and also for the expansion waves at the outlets one has,
[−hL ] u2L
i
h
p
p
1
2
2
2
+ uR + ghR hL uL + g hR − hL − uR + ghR uR hR + uR hR = 0
2
(4.45)
In summary, the new boundary condition scheme is built upon equations (4.37) and (4.39) which describe
shocks or expansion waves, respectively. These equations must be solved with respect to hL or uL , whether
the solution applies to an inlet or an outlet, respectively, yielding equations (4.42), (4.43), (4.44) and
(4.45).
Equation (4.42) is a 4th order polynomial that can be solved analytically by the Cardano-Ferrari’s method
(Jao, Version 5) providing four mathematical solutions (hL1 , hL2 , hL3 , hL4 ) ∈ C. For (4.43), a 3rd order
polynomial, the applied solving method was the Cardano’s formula (Jao, Version 8) yielding three possible
solutions (hL1 , hL2 , hL3 ) ∈ C. Both (4.44) and (4.45) are 2nd order polynomials that are ordinarily solved
with the quadratic formula, rendering two real solution candidates (hL1 , hL2 ) ∈ R.
4.2.6
Physical feasibility of the mathematical solutions
A wide range of solutions can be as cumbersome as no solution at all. Well-posed boundary conditions
should provide unique and smooth values for the PDE’s being solved inside the domain, which arises the
problem of selecting the only physically correct value within the mathematical array of possible solutions,
a particularly sensitive matter when dealing with 4th or 3rd order polynomials. The purpose of the present
section is to derive adequate conditions that can be employed as validation criteria for analytic solutions.
As time advances, a shock should have characteristics intersecting on its x−t path (Figure 4.2a). If these
characteristics are emanating from the shock (Figure 4.2b) then its propagation is unstable even with
small perturbations, thus endangering the well-posedness of the solution. This condition is called the
entropy condition (LeVeque, 1990; Toro, 2009) and is given by,
λL ≥ S ≥ λR
(4.46)
where λL and λR are the characteristics on the left and right sides of the shock, respectively, and S is the
shock propagation speed. Considering the Rankine-Hugoniot conditions, equation (4.32), this condition
can be satisfied only if (Figure 4.2a), (LeVeque, 1990; Toro, 2009)
ϕ(xsL (t), t) ≥ ϕ(xsR (t), t)
(4.47)
43
4. Numerical Model
where ϕ(xsL (t), t) and ϕ(xsR (t), t) are the states on the left and right sides of the shock. This behavior was
observed during the debug tests of these boundary conditions as the Rankine Hugoniot jump conditions
consistently failed to fulfill the entropy condition when hL ≤ hR (Figure 4.2b).
t
λL > S > λR
t
λR > S > λL
S
λR
(a) Entropy-satisfying shock
physically admissible
λR
λR
x
S
λL
S
λL
λL
t
x
-
(b) Entropy-violating shock
physically non-admissible
x
-
(c) Expansion wave - physically
admissible
Figure 4.2: Physical admissibility with respect to the entropy condition
This led to the conclusion that the Rankine-Hugoniot solutions are only physically possible when hL ≥ hR .
When hL ≤ hR an expansion wave occurs and its corresponding entropy condition is given by,
ϕ(xsL (t), t) − ϕ(xsR (t), t)
<
xsL (t) − xsR (t)
t
(4.48)
where t is time and is some constant > 0. Figure 4.2c depicts and expansion wave with a fan of
characteristics emanating from the discontinuity’s initial position (x, t) = (x0 , t0 ). This is consistent with
the entropy condition for an expansion wave, which is given by λR > λL .
4.2.7
Results
The results herein presented were obtained using an horizontal fixed bed with 10 [m] x 25 [m] and a
constant discharge of 5 [m3 /s] where flow resistance was modeled using Manning-Strickler’s Law with a
Ks coefficient of 33 [m1/3 s−1 ]. An unstructured triangulated mesh was used containing approximately
10.000 elements with an average characteristic length of 0.25 [m]. The initial conditions were similar to
a dam-break scenario.
Comparison with the Riemann’s invariant solution
Overall, the comparison with the previously existent boundary conditions yielded very consistent results.
The Rankine-Hugoniot boundary conditions exhibited a slightly faster wave front (Figure 4.3c) but the
general behavior of the solutions was, in practice, nearly coincidental.
Figure 4.3 shows the evolution of a dam-break scenario with an upstream inflow. The shocks at the inlet,
located at x = 0, held a solution that was very close to that of a Riemann’s invariant in most of the time
steps. The outlet was forcing the flow to a minimum energy level.
44
4.2. Boundary conditions
0.30
0.30
Riemann’s Invariant Solution
0.25
Riemann’s Invariant Solution
0.25
Rankine-Hugoniot Solution
Rankine-Hugoniot Solution
0.20
h (m)
h (m)
0.20
0.15
0.15
0.10
0.10
0.05
0.05
0.00
0
5
10
15
20
0.00
25
0
5
10
x (m)
15
20
25
20
25
x (m)
(a) Time t = 0 s
(b) Time t = 5 s
0.30
0.30
Riemann’s Invariant Solution
0.25
0.25
Rankine-Hugoniot Solution
0.20
h (m)
h (m)
0.20
0.15
0.15
0.10
0.10
0.05
0.05
Riemann’s Invariant Solution
0.00
0
5
10
x (m)
15
20
25
0.00
Rankine-Hugoniot Solution
0
5
10
15
x (m)
(c) Time t = 10 s
(d) Time t = 15 s
Figure 4.3: Comparison results for the implemented boundary conditions
Handling flow reversal
The last test was performed to ensure that the boundary conditions could successfully handle flow reversal,
which is to say, velocity vectors that point outward at inlet cells (or inward at outlet cells). In this test
a steady condition was attained and a wave was then forced at the outlet (Figure 4.4a).
The wave traveled upstream (Figure 4.4c), propagating negative velocities (Figure 4.4d), and eventually
collided with the inlet (Figures 4.4g to 4.4j). At this stage the wave was reflected and was now propagating
downstream (Figures 4.4k and 4.4l) until steady conditions were recovered once again (Figures 4.4m and
4.4n).
0.5
0.3
v (m/s)
h (m)
0.4
0.2
0.1
0.0
0
5
10
15
x (m)
(a) Water Level at t = 8 s
20
25
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0
5
10
15
x (m)
20
25
(b) Velocity at t = 8 s
45
4. Numerical Model
0.5
0.3
v (m/s)
h (m)
0.4
0.2
0.1
0.0
0
5
10
15
x (m)
20
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0
25
(c) Water Level at t = 11 s
v (m/s)
h (m)
0.4
0.3
0.2
0.1
5
10
15
x (m)
20
0
25
0.3
v (m/s)
h (m)
0.4
0.2
0.1
10
15
x (m)
20
0
25
0.3
v (m/s)
h (m)
0.4
0.2
0.1
10
15
x (m)
(i) Water Level at t = 17 s
46
10
15
x (m)
20
25
5
10
15
x (m)
20
25
20
25
(h) Velocity at t = 15 s
0.5
5
5
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
(g) Water Level at t = 15 s
0.0
0
25
(f) Velocity at t = 14 s
0.5
5
20
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
(e) Water Level at t = 14 s
0.0
0
10
15
x (m)
(d) Velocity at t = 11 s
0.5
0.0
0
5
20
25
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0
5
10
15
x (m)
(j) Velocity at t = 17 s
4.2. Boundary conditions
0.5
0.3
v (m/s)
h (m)
0.4
0.2
0.1
0.0
0
5
10
15
x (m)
20
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0
25
(k) Water Level at t = 20 s
0.3
v (m/s)
h (m)
0.4
0.2
0.1
5
10
15
x (m)
10
15
x (m)
20
25
20
25
(l) Velocity at t = 20 s
0.5
0.0
0
5
20
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
-0.6
0
25
(m) Water Level at t = 60 s
5
10
15
x (m)
(n) Velocity at t = 60 s
Figure 4.4: Steady conditions at time t = 100 [s] and t = 2000 [s]
Attaining steady conditions
The implemented boundary conditions had suitably attained a steady flow after a simulation time of
roughly 100 seconds. At the end of the simulation, 2000 seconds, the free-surface profile was exactly the
same, ensuring the well-posedness of this new feature.
0.5
0.4
h (m)
0.3
0.2
Riemann’s Invariant Solution
0.1
0.0
Rankine-Hugoniot Solution
0
5
10
x (m)
15
20
25
Figure 4.5: Steady conditions at time t = 100 s (still the same at t = 2000 s)
47
4. Numerical Model
48
Chapter 5
Integration with Geographic
Information Systems
5.1
Introduction
In the interest of providing a better (graphical) user interaction and also promoting interoperability
between STAV-2D and geographic information systems (GIS) some auxiliary sub-programs (Add-Ons)
were developed, allowing the user to define all the necessary geometric entities for mesh generation and
refinement within a GIS environment.
STAV2D current Add-Ons only feature ArcMap, a geospatial processing program and a main component
of the broadly known ESRI’s ArcGIS suite. Some of the most powerful attributes of ArcMap include
viewing, plotting, editing, creating and analyzing geospatial data in either matrix or vector formats.
These features highlight the advantages of having prompt information exchanges between STAV-2D, GIS
and a mesh generator which was, in this case, the open-source code Gmsh (Geuzaine & Remacle, 2009).
5.2
5.2.1
Pre-Processing utilities
Generating boundaries
Defining complex domain geometries, for instance isoline delimited boundaries, can become rather cumbersome if carried on a point-by-point basis. The purpose of this Add-On was to quickly convert boundary
polygons drawn in ArcMap to a format compatible with Gmsh. This enables the user to directly define
the boundaries over the altimetric dataset, avoiding the manual transfer of (x, y) coordinates into Gmsh.
The usage of this Add-On is quite simple: (the following steps only apply to ArcMap)
1. Draw the boundary polygon in ArcMap using a shapefile (.shp) format.
2. Apply the Feature Vertices to Points tool to the .shp layer
49
5. Integration with Geographic Information Systems
3. Apply the Add XY Coordinates to the .shp layer
4. Export the Attribute Table to ./Input/Boundary.txt
5. Start GIS2Gmsh.exe and follow on-screen instructions
6. Use the generated file (boundary.geo) in ./Output as a geometry file in Gmsh
7. Use Gmsh to mesh the domain inside the boundary delimited in 1.)
(a) Domain boundary polygon drawn in ArcMap
using .shp format (yellow line)
(b) Domain boundary polygon converted to Gmsh’s
.geo format
Figure 5.1: Pre-processing samples: example of ArcMap isoline defined boundary (left) and conversion
to Gmsh (right)
5.2.2
Refining meshes
Mesh detail levels can greatly influence the performance of STAV-2D. In the case of greater detail requirements at singular locations one should be able to locally refine the mesh without having to do so
to the whole domain. The mesh refinement add-on is able to import several refinement polygons from
ArcMap and convert these into a Post-Processing Gmsh file (.POS) that can be applied as a background
refinement mesh. This enables the user to visually define refinement levels, hence promoting productivity,
since there is no need to explicitly parametrize these regions with coding in Gmsh’s geometry files (.geo).
The usage of this add-on can be summarized as:
1. Draw the refinement polygons in ArcMap using shapefiles (.shp).
2. Name the shapefiles sequentially: 1.shp, 2.shp (...)
3. Apply the Feature Vertices to Points tool to all .shp layers.
4. Apply the Add XY Coordinates to all .shp layer
5. Export the Attribute Table of all .shp files to ./Input/X.txt where X is the shapefile number.
6. Modify all X.txt files according to the sample files in "./Input/Templates".
7. Start GIS2Gmsh.exe and follow on-screen instructions
8. Start Gmsh and open the geometry (.geo) file previously generated for the boundary
9. Apply the Post-Processing (Background_mesh.pos) file in ./Output as a background mesh
10. Use Gmsh to mesh the domain using the prescribed cell lengths in 1.)
50
5.3. Post-processing utilities
The following figures only serve illustrative purposes and may not reflect the actual geometries of the simulations presented in Chapter 6.
(a) Refinement polygons (red and white) drawn in
ArcMap
(b) Target element size throughout the domain
(Gmsh, (Geuzaine & Remacle, 2009)))
(c) Mesh example overview
(d) Local refinement detail
(e) Coarser cells at the inlet
(f) Coarser cells at the outlet
Figure 5.2: Example of mesh refinement procedure
5.3
5.3.1
Post-processing utilities
Resampler and raster converter
The VTK resampler is a tool that interpolates flow variables, such as velocity or water height, between
two different meshes. Its most valuable use lies in warming-up roughly meshed problems until steady or
other desired conditions are attained. The obtained solution can then be interpolated to a more refined
mesh and used as an initial condition for further simulations.
51
5. Integration with Geographic Information Systems
The raster converter allows the user to import STAV-2D output into GIS environments so that this data
can be referenced and displayed alongside with national geographic information datasets.
Geometry processing and interpolation algorithms
In these post-processing utilities two distinct interpolation methods were used: a cell-to-cell assignment or
a space-averaged neighborhood assignment. The cell-to-cell assignment method relies on a Ray-Casting
algorithm (Foley et al., 1996) to determine if a given cell of the target mesh is (or not) inside a cell of
the original mesh.
This is solved as a point-in-polygon (PIP) problem (Foley et al., 1996) where the point is the barycenter
of the target mesh cell and the polygon is the cell of the original mesh.
To solve the problem one must draw a ray towards the barycenter of the target mesh cell (a line that
connects the barycenter of that cell to a point at the "infinity") and intersect this ray with the edges of
the original mesh cell. If the number of intersections is odd then the barycenter of the target mesh cell
is inside the original mesh cell, otherwise it is outside.
PIP
∞
In
PIP
∞
Out
Target Mesh Cell
Original Mesh Cell
Original Mesh Cell
Target Mesh Cell
(a) "Inside" PIP condition: odd intersection count
(b) "Outside" PIP condition: even intersection count
Figure 5.3: Point-in-polygon (PIP) problem examples
In the second method, the neighborhood space-average, an Inverse Distance Weighing (IDW) algorithm
is used. Denoting by ~x the target cell barycenter and by ~xi the barycenter of its N nearest neighbors in
the original mesh, the IDW algorithm is expressed as,
A(~x) =
N
X
wi (~x, ~xi )A(~xi )
PN
x, ~xi )
i=1 wi (~
i=1
(5.1)
where A(~x) is the value of the interpolated quantity A assigned to the target mesh cell and A(~x) is the
value of this same quantity inside the nearest original mesh cells. The weights wi are defined as,
wi (~x, ~xi ) =
52
1
D(~x, ~xi )p
(5.2)
5.3. Post-processing utilities
where D(x, xi ) is the distance between the barycenters x and xi and p is a power parameter. If N and p,
which are user inputs, are defined as one and zero, respectively, the IDW algorithm gets reduced to the
also widely used Nearest Neighbor (NN) method.
The usage of this Add-On can be resumed to:
1. Place the VTK files to be converted or resampled in ./Input
2. Name these files sequentially: Input_1.vtk, Input_2.vtk (...)
3. When resampling also include a VTK file with the refined mesh in ./Input and name it Input_New.vtk
4. Start GIS2Gmsh.exe and follow on-screen instructions
5. Raster and resampled .vtk files are saved in ./Output
(a) Original VTK file (Simple testbed with 5.000
cells)
(b) Resampled VTK file (Same testbed upsampled
to 37.000 cells)
(c) Original VTK file (Ponte 25 de Abril close up
on the warm-up simulation using 750.000 cells)
(d) Resampled VTK file (Same location on the
main simulation upsampled to 2.400.000 cells
Figure 5.4: Examples of VTK file resampling results
53
5. Integration with Geographic Information Systems
54
Chapter 6
A 1755 tsunami in today’s Tagus
estuary topography
6.1
6.1.1
Pre-Processing
Lisbon and Tagus estuary topography
The present section focuses on the geographical information datasets employed in this work and their importance in obtaining accurate results since detailed geometries are required to comprehensively describe
wave run-up and inundated areas.
Two altimetric datasets were employed. The most detailed was a 5x5 [m] resolution grid with approximately ninety million cells featuring poorly defined and diffuse geometries but, nonetheless, satisfactory
elevation accuracy for the built environment (see Figures 6.3a and 6.3e for a sample view). A smoother
10x10 [m] resolution grid with twenty-five million cells was also used despite having no apparent built
environment geometries (Luis, 2007) (Figures 6.1, 6.2, 6.3b and 6.3f). The Tagus estuary bathimetry was
also provided as a 10x10 [m] resolution grid (Luis, 2007).
To properly define the built environment geometries these two Digital Elevation Models (DEM) had to
be used alongside with a third information layer. This auxiliary layer was a vector dataset that featured
all the buildings in Lisbon as polygonal entities
1
(see Figure 6.3c and 6.3g). A simple zonal statistic
procedure was then followed to improve the definition of the existing DEMs: each polygon (or "building")
was converted to a single valued raster using the highest elevation of the 5x5 [m] grid inside that same
polygon whereas outside the polyogns (the "streets") the smoother 10x10 [m] grid was used. That is to
say the buildings defined by the combination of vector and 5x5 [m] grid information were overlaid on the
10x10 [m] grid 2 . The obtained results are displayed in Figures 6.3d and 6.3h.
1 This vector dataset was supplied by the Department of Geographic Information and Cadastre (DIGC) of Lisbon Municipality and featured up-to-date building plan geometries.
2 All grids were previously upsampled (using surface interpolation) to 2.0[m] cells to ensure adequate consistency with
very fine meshes (CL ≈ 2.5[m]). The designation of 5x5 or 10x10 [m] resolution only concerns the data measurement grid
and not the actual cell size used for STAV-2D inputs.
55
6. A 1755 tsunami in today’s Tagus estuary topography
The tower foundations of Ponte 25 de Abril were also carefully modeled since these structures may also
influence an incoming tsunami. As for the south bank of the Tagus, the existing harbors were also given
appropriate detail since these areas are expected to have a major impact on the tsunami run-up.
Elevation (m)
0
50
100
150
250
300
Figure 6.1: DEM Raster with 5x5 [m] resolution
Bathymetry (m)
-50
-40
-30
-20
-10
0
Figure 6.2: Bathimetry with 10x10 [m] resolution
56
6.1. Pre-Processing
(a) Lisbon downtown in the 5x5 [m] DEM
(b) Lisbon downtown in the 10x10 [m] DEM
(c) Lisbon downtown in the vector layer
(d) Lisbon downtown in the final DEM
(e) Alcântara in the 5x5 [m] DEM
(f) Alcântara in the 10x10 [m] DEM
(g) Alcântara in the vector layer
(h) Alcântara in the final DEM
Figure 6.3: Example of DEM built environment enhancement
57
6. A 1755 tsunami in today’s Tagus estuary topography
6.1.2
Initial Conditions
In this work two simulations were computed: a warm-up simulation and a main simulation. The warmup simulation purpose is to build a steady flow consistent with the mean annual discharge in the Tagus
river. This setup is aimed at building initial conditions for the main simulation. To do so, nearly constant
velocity and water level readings were required at 4 virtual gauges scattered throughout the estuary (see
Figures 6.4e to 6.4h on page 59). The resulting output (see Figures 6.4a and 6.4b) was then resampled
with the Add-On introduced in Chapter 5 and used as the initial condition for the main simulation.
The initial condition for this warm-up simulation was a still water mass at high tide level (+2 [m]). The
discharge hydrograph at the inlet was not strictly constant (Figure 6.4d on page 59) because the shocks
produced by a sudden discharge of 450 [m3 /s] would likely introduce small waves that could propagate
inland, corrupting the correct wet area and possibly compromising data resampling on waterfront cells.
The outlet condition was a constant high tide (+2 [m]) water level hydrograph.
6.1.3
Boundary Conditions
Upstream
A simple hydrological analysis was performed to obtain the modular flow to be imposed at the inlet (or
upstream) boundary. A mean annual flow duration curve was derived from a 25 year flow time-series
registered at the Ómnias - Sacavém hydrometric station (see Figure 6.5).
5000
Q (m 3/s)
4000
3000
2000
1000
0
0
Q mod = 450 m3/s
2000
4000
6000
8000
10 000
t (d)
Figure 6.5: Mean annual flow duration curve (Ómnias - Sacavém hydrometric station)
The corresponding modular flow is given by,
R tf
Qmod =
t0
Qd (t) dt
tf − t0
≈ 450 [m3 /s]
(6.1)
where Qd (t) is the mean annual flow duration curve, t is time (ranging from t0 to tf ) and Qmod is the
modular flow. The obtained value corresponds to an event with a 1.5 year return period.
58
6.1. Pre-Processing
(a) Water level at the end of warm-up
(b) Velocity magnitude at the end of warm-up
500
Q (m3/s)
400
300
200
100
0
0
t (s)
6000
8000
10 000
0. 5
2.0
0. 4
2.0
0. 4
1.5
0.3
1.5
Velocity
0.3
1.0
0. 2
1.0
Water Level
0. 2
0
2000
0. 1
0.5
0.0
10 000
0.0
0. 1
Velocity
4000
6000
8000
0
2000
4000
6000
8000
t (d)
t (d)
(e) Virtual gauge 1 readings
(f) Virtual gauge 2 readings
0.0
10 000
0. 5
2.5
0. 5
2.0
0. 4
2.0
0. 4
1.5
0.3
1.5
Velocity
0.3
1.0
0. 2
1.0
Water Level
0. 2
Water Level
0.5
0.0
Velocity
0
2000
4000
6000
8000
h (m)
2.5
0. 1
0.5
0.0
10 000
0.0
V (m/s)
0.0
Water Level
V (m/s)
2.5
h (m)
0. 5
V (m/s)
2.5
0.5
h (m)
4000
(d) Warm-up inlet discharge hydrograph
V (m/s)
h (m)
(c) Virtual gauges location
2000
0. 1
0
2000
4000
6000
8000
t (d)
t (d)
(g) Virtual gauge 3 readings
(h) Virtual gauge 4 readings
0.0
10 000
Figure 6.4: Warm-up setup and results
59
6. A 1755 tsunami in today’s Tagus estuary topography
Downstream
The downstream boundary condition was obtained from offshore propagation simulations performed at
the University of Lisbon facilities. The model employed was a modified version of the COMCOT code
by Wang & Liu (2006). The provided data describes the tsunami at the atlantic outlet arc. This data
was supplied pointwise and was implemented in the perspective of a worst case scenario: the most severe
wave series was imposed at all outlet cells. Several combinations were attempted to ensure that the wave
propagating inside the domain represented the worst case scenario within the supplied data.
Figures 6.6a to 6.7b display the supplied data as a function h(t). In order to save computation time
the first 2000 seconds were ignored since there are no relevant waves within this time interval, which is
to say, the hydrograph shown in Figure 6.6a was shifted by 2000 seconds towards the origin t0 = 0[s].
Simulation time was 10400 seconds (3 hours) and allowed for a complete propagation and draw back of
all the waves from the hydrograph depicted in Figure 6.6a. If there was no data available, t ≥ 7200[s], a
still high tide level was imposed (+2[m]).
10
8
h (m)
6
4
2
0
0
2000
4000
6000
8000
10 000
t (s)
10
10
8
8
6
6
h (m)
h (m)
(a) Tsunami level imposed at the atlantic boundary
4
2
0
2000
4
2
2200
2400
t (s)
2600
2800
3000
(b) First and main wave of the tsunami event
0
4700
4800
4900
5000
(c) Second wave
Figure 6.6: Supplied data for a 1755-like tsunami
60
t (s)
5100
5200
10
10
8
8
6
6
h (m)
h (m)
6.2. Brief framework for result presentation and discussion
4
2
4
2
0
5550
5600
5650
t (s)
5700
(a) Third wave
5750
5800
0
6950
7000
7050
t (s)
7100
7150
7200
(b) Fourth and fifth waves
Figure 6.7: Supplied data for a 1755-like tsunami (2)
6.2
Brief framework for result presentation and discussion
The previous chapters introduced the mathematical and theoretical concepts upon which the STAV-2D
model is built and all the information layers that compose the framework of a potential tsunami scenario
in Lisbon. The results herein presented are mainly focused on the effects of the first wave of a 1755-like
tsunami (Figure 6.6b), by far the most destructive and the earliest of the whole event.
Display and discussion of the obtained results is carried at two distinct levels: an overview of the main
morphological impacts and run-up extents at the seafronts of Carcavelos, Oeiras, Belém, Caparica and
Barreiro and a more detailed street-level analysis at Alcântara, Lisbon Downtown and Cacilhas. This
last set of results ads great value to the tsunami modelling efforts carried by the Portuguese scientific
community since no previous studies were known to explicitly discretize the built environment features.
Another important contribution arises from the inclusion of solid transport capabilities in STAV-2D,
allowing debris flow simulations to be performed. Recent events showed that the hydrodynamics of
tsunami inland propagation is much closer to that of a debris flow than to a clear water flow. The last
section in this chapter is devoted to debris scattering patterns and main morphological impacts.
Figure 6.8: Relevant locations for result presentation and discussion
61
6. A 1755 tsunami in today’s Tagus estuary topography
6.3
6.3.1
Overview of a tsunami scenario in Lisbon
General description of tsunami propagation
Figure 6.9 depicts the propagation of the first tsunami wave for 27 minutes. The description of the first
wave at the Tagus mouth is well documented (Baptista et al., 2011) in a magazine dated 19 November
1755 (Urban, 1755): “[...] Castle of Bugio was so far overcome by the water, that the garrison fired several
guns as signals for help [...]. By my best judgment, the water rose in five minutes about 16 feet [...]”.
From this report the wave height at Bugio Tower should be approximately 5 [m]. After hitting Bugio
the tsunami propagates towards the seafront of Oeiras causing massive run-ups, at very high speeds.
The wave then propagates inward the estuary through the Tagus channel, accelerating due to the higher
depths while flooding Lisbon and Cacilhas, a particularly hazardous stage of this event.
(a) Water level as the tsunami hits Bugio Tower
(b) Water level, t = 7 minutes
(c) Water level, t = 14 minutes
(d) Water level, t = 17 minutes
(e) Water level, t = 20 minutes
(f) Water level, t = 27 minutes
Water level (m)
2.0
3.0
4.0
5.0
6.0
Figure 6.9: Tsunami propagation overview at the Tagus estuary
Tsunami impact at Bugio Castle used as relative time reference
62
6.3. Overview of a tsunami scenario in Lisbon
When the tsunami reaches the estuary it propagates radially and the wave height drastically reduces as
a consequence of the widening cross section. Nonetheless, it still carries immense energy with significant
run-ups occurring at Barreiro, on the south bank of the estuary.
6.3.2
Impact at Carcavelos
Carcavelos is one of the busiest beaches in Lisbon. The prominent shape of the seafront at this location
makes it one of the most hazardous areas in the case of a tsunami. The wave does not undergo any kind
of deflection or attenuation prior to impact, which is nearly head-on with the beach. This causes an
extensive and destructively fast inundation. Maximum run-up is as high as 20 [m] above mean-sea level
(MSL) and may extend to 600 [m] away from shore, with velocities just above 5 [m/s] (Figure 6.10).
2.0
(a) Water level, t = -1 minutes
(b) Velocity magnitude, t = -1 minutes
(c) Water level, t = 1 minutes
(d) Velocity magnitude, t = 1 minutes
(e) Water level, t = 3 minutes
(f) Velocity magnitude, t = 3 minutes
Water level (m)
Velocity (ms-1)
3.0
4.0
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.10: Tsunami run-up at Carcavelos
Tsunami impact at Bugio Castle used as relative time reference
63
6. A 1755 tsunami in today’s Tagus estuary topography
6.3.3
Impact at Caxias and Cruz Quebrada
The head-on impact of the tsunami at this location (Figure 6.9 on page 62) inflicts massive wave heights
and velocities when the tsunami propagates overland. Run-up is expected to propagate inland until the
ground rises to 10 [m]. The very mild topography of the reaches at Caxias and Cruz Quebrada implies
that a 10 [m] run-up can lead to an inundated area as far as 1500 [m] from the seafront. A mild terrain
inflicts much less resistance to the incoming wave leading to a maximum velocity of 6 [m/s]. Algés is also
a vulnerable area with a inundation of 300 [m] inland, up to 8 [m] high.
Inland penetration reaching further than 1 [km] is only observed at Caxias. Poor built environment
definition (and low density) can also favor these extensive inundations since macro flow resistance is
practically ignored.
2.0
(a) Water level, t = 6 minutes
(b) Velocity magnitude, t = 6 minute
(c) Water level, t = 7 minutes
(d) Velocity magnitude, t = 7 minutes
(e) Water level, t = 8 minutes
(f) Velocity magnitude, t = 8 minutes
Water level (m)
Velocity (ms-1)
3.0
4.0
5.0
6.0
0.0
1.0
Figure 6.11: Tsunami run-up at Caxias and Cruz Quebrada
Tsunami impact at Bugio Castle used as relative time reference
64
2.0
3.0
6.3. Overview of a tsunami scenario in Lisbon
6.3.4
Impact at Belém and Algés
Belém is recognized for its concentration of national monuments and public spaces, composing a mixture
of historical monuments and modern symbols of the Portuguese culture. This combination of iconic
landmarks makes it one of the most attractive locations in the whole city, bringing the discussion of a
tsunami impact at this location to a very sensitive level. Mosteiro dos Jerónimos, Torre de Belém, Centro
Cultural de Belém and Padrão dos Descobrimentos are all at jeopardy due to their riverfront location,
facing inundation depths of 3.5 [m] and velocities higher than 5 [m/s]. Run-up extent is limited to 550
[m] from the riverfront, up to 5 [m] high (Mosteiro dos Jerónimos, Figure 6.12).
2.0
(a) Water level, t = 10 minutes
(b) Velocity magnitude, t = 10 minute
(c) Water level, t = 11.5 minutes
(d) Velocity magnitude, t = 11.5 minutes
(e) Water level, t = 13 minutes
(f) Velocity magnitude, t = 13 minutes
Water level (m)
Velocity (ms-1)
3.0
4.0
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.12: Tsunami run-up at Belém and Algés
Tsunami impact at Bugio Castle used as relative time reference
65
6. A 1755 tsunami in today’s Tagus estuary topography
6.3.5
Impact at Caparica and Trafaria
Another rather hazardous area is the Atlantic seafront between Costa da Caparica and Trafaria. At south
of Costa da Caparica the wave will propagate inland as further as 800 [m]. Smaller inundation lengths
occur at Costa da Caparica and at Almada’s northern seafront due to dense built environment and forest
areas. Another concerning detail of this seafront is the presence of many camping sites bordering beach
areas. These heavily populated sites, namely during bathing seasons, are worryingly prone to catastrophic
life losses due to inherent evacuation and access difficulties.
Important wave shoaling is expected due to the mild slope of the seabed. Preceding further inland
propagation this seafront will be devastatingly swept by a 10 [m] wave at a speed of roughly 7 [m/s]. A
later discussion on this chapter is reserved to the severe morphological changes that are to occur.
2.0
(a) Water level, t = 1 minutes
(b) Velocity magnitude, t = 1 minutes
(c) Water level, t = 4 minutes
(d) Velocity magnitude, t = 4 minutes
(e) Water level, t = 7 minutes
(f) Velocity magnitude, t = 7 minutes
Water level (m)
Velocity (ms-1)
3.0
4.0
5.0
6.0
0.0
1.0
Figure 6.13: Tsunami run-up at Caparica and Trafaria
Tsunami impact at Bugio Castle used as relative time reference
66
2.0
3.0
6.4. Street-level results for inland propagation
6.4
6.4.1
Street-level results for inland propagation
Results for Alcântara
As a result of the 1755 earthquake and tsunami, the Portuguese Royal Family and the government moved
to Alto da Ajuda. At that time, and following the catastrophe, it was the best place to hold the national
power entities: shallow bedrocks allowed better seismic structural behaviour as foundations could be
embedded, the elevation was enough to surpass any incoming tsunami and airy surroundings protected
the Royal family from any developing plague due to the many fatalities in the lower city. Nobility, workers
and municipal services followed this relocation. The increasingd traffic between Alto da Ajuda and the
old Lisbon city allowed this otherwise suburban area to expand at a fast pace and develop into what is
nowadays one of the most central spots in Lisbon.
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.14: One-way traffic zones (top) and tsunami approach to Alcântara (bottom; t=10’20”)
Tsunami impact at Bugio Castle used as relative time reference
Known to have suffered severe damage from the 1755 event this area was refined with a 2.5 [m] mesh
and the built environment was modeled with actual city planimetry. The urban mesh is quite dense and
organically adapted to the underlying topography. This poses some difficulties to traffic management
and, to avoid typical rush-hour bottlenecks, many loops are implemented with restrictive one-way traffic
streets. These are especially worrisome as they can drive traffic towards the riverfront. Figure 6.14 shows
5 concerning zones in what respects to restrictive circulation.
67
6. A 1755 tsunami in today’s Tagus estuary topography
Roughly 10 minutes after the tsunami hits Bugio Tower, the water level has risen to the point of riverfront
overtopping (Figure 6.14 on 67) and a minor 0.5 [m] high inundation front takes place at Avenida da
Índia. At this stage velocities range from 0.3 to 2.5 [m/s]. Impact at Cordoaria Nacional is imminent.
Two minutes later (Figure 6.15 on page 69), at 12 minutes and 20 seconds, most of the riverfront between
Cordoaria Nacional and Lisbon harbour is already flooded. At this stage the inundation front already
approaches Hospital Egas Moniz as Rua da Junqueira and Congress Center are swept by a devastating
wave front with nearly 1.6 [m] at a speed of almost 4 [m/s]. Vila Galé Ópera Hotel and Alcântara/Santo
Amaro Docks are also severely hit. At the hotel water depths are comparable to a minor flooding, about
0.6 [m], but with velocities over 4 [m/s] heavy debris can be carried inside the hotel halls, causing many
casualties. At Alcântara Docks the marina dampens the wave and slows down the flow at the cost of
higher depths, which can reach 1.5 [m]. The harbour is a concerning location as high flow depths and
great velocities can carry a large amount of containers, and even ships, towards the streets.
Water is already reaching the inner grounds of Hospital Egas Moniz by the time of 14’20” (Figure 6.16
on page 70). Patient evacuation becomes severely constrained since three of the four surrounding streets
are one-way minor arteries and lead towards the inundated areas, namely Calçada da Boa Hora, Travessa
Conde Ribeira (Zone 3) and Rua da Junqueira (Zone 1 - completely flooded at this stage). Only one
street can be used to effectively reach higher grounds: Rua Filipe Vaz towards Travessa Giestal (Zone
1). Alcântara-Mar Train Station is inundated with flow depths reaching more than 1 [m] and velocities
around 3 [m/s]. According to traffic circulation allowances Zones 1, 2, 3 and 4 are now completely
isolated with possible access to Zones 1, 3 and 4 only possible if such rules are neglected. Zone 2 is
completely surrounded by water. Avenida 24 de Julho gets inundated by a fast advancing front as the
flow breaks through the maintenance facilities on the east area of the harbour and heads towards the
National Museum on Ancient Art. A significant draw-down starts shaping at Cordoria Nacional with
flow velocity towards the river reaching 2.5 [m/s] and flow depths of approximately 1.7 [m].
Two minutes later, 16’20” (Figure 6.17 on page 71), the tsunami crest has already left the harbour area and
now heads towards Downtown. Nonetheless, the tsunami hazard is still far from stabilized at Alcântara.
Draw-downs can be as destructive as run-ups and all the area is now subjected to a generalized drawdown. Velocities higher than 3 [m/s] are observed at the Alcântara Docks and Congress Center. Through
all inundated area the flow patterns are clearly affected by the buildings, acting as obstacles, as seen in the
velocity map. Inundation extent is at its maximum and nearly coincides with the 5 [m] elevation isoline
and Avenida 24 de Julho. Water has now reached Santos, approximately half-way between Alcântara
and Downtown, with water depths around 1.7[m] and velocities of 2.4 [m/s]. Unusual flow patterns occur
inside the harbour as two noticeably large eddies begin to form (Figure 6.17b).
After another four minutes the draw-down is still intense at the harbour and riverfront areas with outgoing
flow reaching 3 [m/s] (Figure 6.17, page 71). Depths are generally below 1 [m]. The flow pattern
displayed in the velocity map clearly shows the main drawing-down flow trajectories. The buildings act
as bottlenecks inducing high velocities mainly at the harbour area, including the Alcântara Docks.
68
6.4. Street-level results for inland propagation
Tsunami impact at Alcântara waterfront. 12’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 12 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 12 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.15: Detailed results for Alcântara, t = 12 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
69
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Alcântara waterfront. 14’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 14 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 14 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.16: Detailed results for Alcântara, t = 14 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
70
3.0
6.4. Street-level results for inland propagation
Tsunami impact at Alcântara waterfront. 16’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 16 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 16 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.17: Detailed results for Alcântara, t = 16 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
71
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Alcântara waterfront. 20’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 20 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 20 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.18: Detailed results for Alcântara, t = 20 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
72
3.0
6.4. Street-level results for inland propagation
6.4.2
Results for Downtown
The Downtown area comprises the Pombaline Town, Cais do Sodré and the Alfama district beneath the
Lisbon Castl e, and extends northwards towards Rossio and Figueira squares and Avenida da Liberdade,
a well-known boulevard noted for the shops and cafes. The Pombaline Lower Town is a distinct district
built after the 1755 Lisbon earthquake and tsunami. The tsunami impact is well documented and is
known to have inundated, for instance, Largo de São Paulo, Terreiro do Paço and Cais do Sodré. Further
details on the historical details of the 1755 event can be found in Baptista et al. (2006).
This area is one of the first examples of earthquake-resistant construction, chareacterised by buildings
featuring a technique known as gaiola pombalina Bento et al. (2004). Its modern layout features a narrow
but regular grid where traffic circulation is, presently, arranged in several looping cells that can become
worrisome when traffic is forced to head towards the riverfront. Figure 6.19 shows several streets that
can pose difficulties in safely reaching higher grounds.
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.19: One-way traffic zones (top) and tsunami approach to Downtown (bottom; t=13’20”)
Tsunami impact at Bugio Castle used as relative time reference
73
6. A 1755 tsunami in today’s Tagus estuary topography
Figure 6.19 shows the tsunami approaching Cais do Sodré. At 13’20” water level is beginning to rise
to the point of overtopping Avenida 24 de Julho. The inundation front rapidly advances and less than
two minutes later, 15’ after the tsunami hit Bugio Tower, Cais do Sodré has been almost completely
overwhelmed by the tsunami. Rua do Instituto Industrial is already flooded half-way despite being
almost 250 [m] away from the riverfront. The vulnerable Avenida Ribeira das Naus is also overtopped,
blocking one of the main exit routes from Downtown to Northern riverfront areas. The inundation front
is advancing toward Terreiro do Paço and a nearby fluvial interface. Arsenal da Marinha is also severely
exposed. The traffic in Zone 1 and Zone 2 is particularly critical since some streets force traffic towards
the incoming wave, namely Boqueirão Ferreiros, Boqueirão do Duro, Rua da Moeda and Rua do Arsenal.
Zone 3 is a converging point for the traffic from Zones 1, 2, 4 and 5 which can lead to a complete clogging
as Zone 3 is one of the first areas to become inundated. Traffic would be completely paralyzed in all five
zones.
Fourty seconds later the extent of the inundation area is quite disturbing (Figure 6.20 on page 75).
Terreiro do Paço is midway submerged with water depths reaching 1.5 [m] and velocities as great as 3.0
[m/s]. Arsenal da Marinha is also severely hit. Run up has penetrated inland as deep as 380 [m] at Rua do
Instituto Industrial and Boqueirão Ferreiros. Cais do Sodré is now completely flooded with water depth
reaching almost 2 [m] and velocities can be higher than 3 [m/s]. Zones 1, 2 and 3 are now completely
isolated. The main concern is the traffic that is driven towards the riverfront and Cais do Sodré due
to a few one-way streets, namely Rua do Arsenal, Rua do Comércio, Rua do Ouro and Rua da Prata.
Many administrative facilities are at jeopardy of imminent tsunami impact, including the Ministry of
Agriculture, the Ministry of Justice, the Ministry of Finance, Lisbon’s City Hall and several courts. The
urban grid next to Terreiro do Paço is also a very important financial area with several major national
and international banks holding their headquarters at Rua do Comércio and Rua S. Julião.
By 16’20” Terreiro do Paço has now been completely inundated as the wave front heads north, towards
the Pombaline Lower Town grid (Figure 6.21 on page 76). Most administrative facilities are, at this stage,
completely overwhelmed by the incoming tsunami. Inundation fronts are now taking over east downtown
with valuable churches and museums being consumed by the incoming water. All the streets are two-way
and should promote an effective evacuation. The tsunami crest is now ashore of Cais do Sodré, meaning
that the inundation is still being fed at this location. The water propagates deeper inside the urban mesh
at Cais do Sodré reaching Rua da Boavista. A moderate draw-dawn at Santos is also evident.
The Tsunami crest is shorewards of Alfama by the 17th minute, leaving extensive inundation and strong
draw-downs at Cais do Sodré and Terreiro do Paço (Figure 6.22 on page 77). Largo de São Paulo was
known to have been flooded by the 1755 tsunami and the obtained results are consistent with this report.
Rua do Comércio, Rua São Julião and Rua Augusta are at this stage suffering minor flooding, with water
depths of 0.6 [m] and velocities of approximately 0.7 [m/s]. Impact at Santa Apolónia Train Station is
considerably strong, with flow depths of nearly 1 [m] and a velocity around 1.9 [m/s]. Inundated extent
is, again, nearly coincident with the 5 [m] elevation isoline. At the 20th minute (Figure 6.24 on page
79) the tsunami is heading towards the estuary interior. Generalized draw-down occurs throughout the
74
6.4. Street-level results for inland propagation
riverfront, especially in narrow sections between buildings, where velocities can reach almost 4 [m/s].
Tsunami impact at Downtown waterfront. 15’ after hitting Bugio Tower
(a) Water level at Alcântara, t = 15 minutes
(b) Velocity magnitude at Alcântara, t = 15 minutes
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.20: Detailed results for Downtown, t = 15 minutes
Tsunami impact at Bugio Castle used as relative time reference
75
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Downtown waterfront. 15’40” after hitting Bugio Tower
(a) Water level at Downtown, t = 15 minutes 40 seconds
(b) Velocity magnitude at Downtown, t = 15 minutes 40 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.21: Detailed results for Downtown, t = 15 minutes 40 seconds
Tsunami impact at Bugio Castle used as relative time reference
76
3.0
6.4. Street-level results for inland propagation
Tsunami impact at Downtown waterfront. 16’20 after hitting Bugio Tower
(a) Water level at Downtown, t = 16 minutes 20 seconds
(b) Velocity magnitude at Downtown, t = 16 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.22: Detailed results for Downtown, t = 16 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
77
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Downtown waterfront. 17’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 17 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 17 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.23: Detailed results for Downtown, t = 17 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
78
3.0
6.4. Street-level results for inland propagation
Tsunami impact at Downtown waterfront. 20’20” after hitting Bugio Tower
(a) Water level at Alcântara, t = 20 minutes 20 seconds
(b) Velocity magnitude at Alcântara, t = 20 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.24: Detailed results for Downtown, t = 20 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
79
6. A 1755 tsunami in today’s Tagus estuary topography
6.4.3
Results for Almada
Almada has a strong tie with the Tagus river and fluvial activities, mainly due to several shipyards and
transportation services across the Tagus. The wave starts to approach Cacilhas 14’30” after the impact
at Bugio (Figure 6.25 on page 81). A very fast inundation front will hit Rua do Ginjal, where velocities
will reach 2.85 [m/s] and flow depths of about 1.6 [m]. The inland propagation towards Largo Alfredo
Dinis is very fast as flow velocity can reach 3.68 [m/s] and flow depths will range from 0.7[m] to 1.12
[m]. At this stage the inundation reaches Avenida Aliança Povo MFA with a water column 0.81 [m]
high and a velocity of 2.14 [m/s]. This advancing front is rather concerning as it will prevent access to
northern higher grounds that would otherwise be easily accessible by taking the Avenida Aliança Povo
MFA towards Rua Cândido dos Reis. The only possible route to escape the inundation is Avenida Alinaça
Povo MFA towards Avenida António José Gomes.
Ninety seconds later, 16’10” (Figure 6.26 on page 82) after Bugio is hit, the inundation front has already
clashed with the Fireman Headquarters at Cacilhas, where speeds can be as high as 2.61 [m/s] and
flow depths will be around 1.54 [m], the highest depths so far observed at Cacilhas. This is a critical
loss for the response mechanisms since Fireman will not be able to help this population even after a
complete draw-down. The advancing front rapidly propagates through Avenida Aliança Povo MFA, with
a maximum velocity of 3.61 [m/s] and flow depths as high as 1.14 [m]. At LISNAVE (naval warehouses)
the inundation front has a flow depth of 0.73 [m] and a velocity of 3.33 [m/s]
Figure 6.27 (page 83) shows the inundation at 17’40”. The advancing front is already claiming the
Professional School of Almada with flow depths of 1.08 [m] and velocities reaching 2.81 [m/s] at Avenida
Aliança Povo MFA. The maximum flow depth and velocity at this stage occur at the LISNAVE shipyards,
reaching 1.73 [m] and 4.21 [m/s], respectively. Lisbon Naval is struck by a water column with 0.83 [m]
propagating at 2.71 [m/s].
Avenida Aliança Povo MFA should be completely flooded by 22’50” (Figure 6.28 on page 84). The
available alternatives are restricted to Rua Manuel Febrero e Avenida António José Gomes. where the
0.99 [m] flood is still advancing at a speed of 2.64 [m/s]. Lisbon Naval Base is now completely overwhelmed
the incoming water that are now heading towards Rua Manuel José Gomes.
The wave-front heads towards Largo 5 de Outubro, leaving Rua Tenente Valadim and Rua Francisco
Ferrer completely flooded in a matter of just two minutes (Figure 6.28 on page 84). Flow depths of 0.69,
0.79 and 1.12 [m] are registered at Rua José Alves de Almeida, Rua António José Gomes and Lisbon
Naval Base, respectively. Also, velocities can range from 0.45 to 2.16 [m/s] at these locations. Significant
draw-downs occur at the northern areas, namely the LISNAVE shipyards.
The tsunami crest leaves Cacilhas shoreline when 27’40” have passed since the impact on Bugio Tower.
The progress of the inundation is almost unnoticeable as there is only one new location to add to the
flood extent, Largo 5 de Outubro. Flow depths read 0.51 [m] and velocities range from 0.4 to 1.45 [m/s].
Figure 6.28 (on page 84) shows the maximum run-up extent.
80
6.4. Street-level results for inland propagation
Tsunami impact at Cacilhas waterfront. 14’30” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 14 minutes 30 seconds
(b) Velocity magnitude at Cacilhas, t = 14 minutes 60 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.25: Detailed results for Cacilhas, t = 14 minutes 30 seconds
Tsunami impact at Bugio Castle used as relative time reference
81
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Cacilhas waterfront. 16’10” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 16 minutes 10 seconds
(b) Velocity magnitude at Cacilhas, t = 16 minutes 10 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.26: Detailed results for Cacilhas, t = 16 minutes 10 seconds
Tsunami impact at Bugio Castle used as relative time reference
82
3.0
6.4. Street-level results for inland propagation
Tsunami impact at Cacilhas waterfront. 17’40” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 17 minutes 20 seconds
(b) Velocity magnitude at Cacilhas, t = 17 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.27: Detailed results for Cacilhas, t = 17 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
83
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Cacilhas waterfront. 22’50” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 22 minutes 50 seconds
(b) Velocity magnitude at Cacilhas, t = 22 minutes 50 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.28: Detailed results for Cacilhas, t = 22 minutes 50 seconds
Tsunami impact at Bugio Castle used as relative time reference
84
3.0
6.4. Street-level results for inland propagation
Tsunami impact at Cacilhas waterfront. 24’20” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 24 minutes 20 seconds
(b) Velocity magnitude at Cacilhas, t = 24 minutes 20 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.29: Detailed results for Cacilhas, t = 24 minutes 20 seconds
Tsunami impact at Bugio Castle used as relative time reference
85
6. A 1755 tsunami in today’s Tagus estuary topography
Tsunami impact at Cacilhas waterfront. 27’40” after hitting Bugio Tower
(a) Water level at Cacilhas, t = 27 minutes 24 seconds
(b) Velocity magnitude at Cacilhas, t = 27 minutes 40 seconds
Water level (m)
2.0
3.0
4.0
Velocity (ms-1)
5.0
6.0
0.0
1.0
2.0
Figure 6.30: Detailed results for Cacilhas, t = 27 minutes 40 seconds
Tsunami impact at Bugio Castle used as relative time reference
86
3.0
6.5. Morphological impacts
6.5
6.5.1
Morphological impacts
The relevance of flow-morphology interaction
Tsunami overland propagation is most of the times modelled with clear water equations where bed
interaction is restricted to flow resistance. Solid transport features implemented in STAV-2D attempt to
simulate a noticeable feature of this stage of a tsunami, as seen in most footage of recent events, which
is the inclusion of debris in the flow. Erodible beds are expected to suffer major morphological changes
if inland propagation attains significant speeds.
This section provides insights on the aftermath of a 1755-like tsunami, focusing on the morphological
changes that were left behind. One should note that these morphological changes might not refer only to
bed and riverfront configurations but also to built environment appliances and vehicles, which are often
carried further inland by the incoming flow.
6.5.2
Vale do Zebro Royal complex
This wheat processing facility borders the Coina river, near Barreiro, on the south of the Tagus estuary.
It was an important asset during the Portugese Discovoeris era, being owned and directly managed by
the Portuguese Crown, comprising over 27 ovens for bread and biscuit production. Its location was
strategically chosen: nearby woods and cereal fields would provide fuel and wheat for the production
process and the Coina river would serve as a convenient fluvial runway. This complex is known to have
been destroyed by the combination of earthquake and tsunami that occurred on the November 1st of
1755. In 1961 it was converted into a military academy and museum for the Portuguese Marines.
Figure 6.31: Tsunami approaching the Royal Complex at Vale do Zebro, t=25’
Important morphological changes that occurred due to the tsunami are of particular interest to calibrate
the solid transport parameters required for STAV-2D. The aim of this procedure was to obtain a run-up
consistent with relevant morphological changes in the river bed along with significant momentum in the
flooding of the complex.
87
6. A 1755 tsunami in today’s Tagus estuary topography
Prior to reaching Barreiro, the tsunami had already suffered major diffraction on the Tagus mouth and
channel. Furthermore, the wider cross section of the estuary also contributed to lower energy fluxes,
decreasing the wave amplitude and speed. The propagation of the tsunami along the Coina river is
depicted in Figure 6.32 (image sequence is subsequent to Figure 6.9).
2.0
(a) Water level, t = 29’
(b) Velocity magnitude, t = 29’
(c) Water level, t = 39’
(d) Velocity magnitude, t = 39’
(e) Water level, t = 49’
(f) Velocity magnitude, t = 49’
Water level (m)
Velocity (ms-1)
3.0
4.0
5.0
6.0
0.0
1.0
2.0
3.0
Figure 6.32: Tsunami run-up at Coina river (Vale do Zebro)
Tsunami impact at Bugio Castle used as relative time reference
Figure 6.33 (page 89) displays the time-evolving bed elevation and water level curves at a cross section
that intersects the Royal complex. The building itself is represented by a noticeable prominence on the
left bank (closeup figures on the left column). The advancing front has a severe impact on the ground due
to aggressive scouring (over 0.5 [m]) around the factory (Figure 6.33, t=37’). Water depths are around
0.5 [m] at this stage and velocities nearly reach 2 [m/s]. Roughly ten minutes after, the factory is almost
completely destroyed (Figure 6.33, t=47’). The draw-down is also an important stage of this event as it
might induce great velocities. At the complex the draw-down does significant damage to the remnants
of the factory while inflicting further erosion on the riverbank (Figure 6.33, t=37’).
88
5
4
Z m
10
8
6
4
2
0
-2
-4
-1000
3
2
1
-500
X m
0
500
0
0
1000
100
X m
200
300
400
500
(a) Bed (−) and water (−−) levels, t = 37’
(b) Bed (−) and water (−−) levels (closeup)
10
8
6
4
2
0
-2
-4
-1000
5
Z m
4
3
2
1
-500
X m
0
500
0
0
1000
100
X m
200
300
400
500
(c) Bed (−) and water (−−) levels, t = 47’
(d) Bed (−) and water (−−) levels (closeup)
10
8
6
4
2
0
-2
-4
-1000
5
4
Z m
Z m
Z m
Z m
6.5. Morphological impacts
3
2
1
-500
X m
0
500
1000
(e) Bed (−) and water (−−) levels, t = 57’
0
0
100
X m
200
300
400
500
(f) Bed (−) and water (−−) levels (closeup)
Figure 6.33: Tsunami impact at the Royal complex
Tsunami impact at Bugio Castle used as relative time reference
6.5.3
Urban debris washed at Downtown and Alcântara
A simplified setup was used to simulate debris transport at Alcântara and Downtown. A layer of mobile
material was added to the base DEM, which was held fixed during this simulation. This layer was
delimited by the lanes of the main riverfront avenues, namely Avenida da Índa, Avenida 24 de Julho and
Avenida Ribeira das Naus (Figure 6.34). The mobile material was set with the same height and width
as those of an average vehicle, roughly 1.5 [m] and 1.75 [m], respectively, and a density of 1.5 [kg/m3 ],
corresponding to the density of an average car filled with water (up to 75 percent of its volume).
89
6. A 1755 tsunami in today’s Tagus estuary topography
Real
Real
Built Env.
1.75
Modeled
1.5
1.5
1.75
Modeled
Built Env.
H = 1.50
1.5 (m)
W = 1.75
4.5 (m)
Continuous
W = 1.75
Overview
W = 1.75
Cross Section
H = 1.50
W = 1.75
Bedrock layer
Mobile material
Mobile material
Figure 6.34: Idealized setup for urban debris, mainly composed of vehicle-like solid materials
The container storage area at the harbour was also modeled as mobile material, using the same vehiclelike density. The height of these elements was set to that of two overlying containers, roughly 4.8 [m].
Increasing the height further would be pointless since the maximum depth expected at this location was
about 3 to 4 [m]. Figure 6.35 displays the locations where mobile material was overlaid to the base DEM.
A debris specific formulation (Ferreira et al., 2009) was employed when modeling solid transport. The
main purpose of this simulation was to assess deposition patterns, a relevant element to define preferable
rescue routes. The outcome can be seen in Figures 6.36 and 6.37 (pages 91 and 91, respectively). The blue
coloring stands for erosion areas, which must be consistent with Figure 6.35, and red indicates debris
deposition. Debris scattering patterns are within the expectations for this kind of phenomena: large
quantities are carried inland and are deposited in streets and squares while smaller quantities are also
deposited close to the riverfront as water returns to the river during the draw-down stage.
Figure 6.35: Mobile material distribution along the riverfront
The aftermath of this scenario showed relevant debris deposition at both designated areas. At Alcântara,
significant deposits are expected near the Centro de Congressos. Particularly, Travessa Conde Ribeira,
Travessa Guarda and Travessa Galé are specially prone to deposition. Avenida 24 de Julho may be blocked
90
6.5. Morphological impacts
due to debris scattering near Docas de Alcântara and at Santos area. At the harbour, major scattering
occurs and due to the large size of the containers, significant damage and accessibility difficulties are
expected.
Figure 6.36: Debris scattering patterns in Alcântara at the aftermath of the tsunami
Figure 6.37: Debris scattering patterns in Downtoen at the aftermath of the tsunami
The downtown area is plentiful in what concerns squares and open spaces. These locations are specially
prone to debris accumulation and may prevent rescue teams from deploying field hospitals and other
victim support services. Terreiro do Paço, Largo da Casa dos Bicos, Arsenal da Marinha, Praça Dom
Luís I and Jardim de Santos are some of the squares where significant deposition occurs.
A brief comparison between the clearwater and debris formulations is herein presented to assess if the
solid material carried in the water column has any significant influence on the wave front depths and
speeds. The results for three monitoring gauges of water levels and velocities are displayed in Figure 6.38
on page 92.
91
6. A 1755 tsunami in today’s Tagus estuary topography
6.5
5
6.0
4
V ms
h m
5.5
5.0
4.5
4.0
2
1
3.5
3.0
3
15
T min
20
25
0
30
(a) Water level at Jardim de Santos: (-) Debris flow
(- -) clear water
15
T min
20
25
30
(b) Velocity magnitude at Jardim de Santos: (-)
Debris flow (- -) clear water
6.5
5
6.0
4
V ms
h m
5.5
5.0
4.5
4.0
2
1
3.5
3.0
3
15
T min
20
25
0
30
(c) Water level at Cais do Sodré: (-) Debris flow (-) clear water
15
T min
20
25
30
(d) Velocity magnitude at Cais do Sodré: (-) Debris
flow (- -) clear water
6.5
5
6.0
4
V ms
h m
5.5
5.0
4.5
4.0
2
1
3.5
3.0
3
15
T min
20
25
30
(e) Water level at Terreiro do Paço: (-) Debris flow
(- -) clear water
0
15
T min
20
25
30
(f) Velocity magnitude at Terreiro do Paços: (-)
Debris flow (- -) clear water
Figure 6.38: Comparison between debris and clearwater formulations
At Terreiro Paço significant differences are observable in velocity measurements. Such differences are
due to a barrier effect on the incoming flow since the gauge was relatively close to Avenida Riberida das
Naus, an avenue were mobile material was added. The remaining readings reflect what was expected. In
the case of a debris flow, the wave propagating inland is slower and higher, which leads to later arrival
times for the debris flow advancing front. Nonetheless, its increased mass and momentum are by far more
devastating. The steady water level in the presented charts represents the inundation that remains after
the wavefront passage. In the case of a debris flow, this inundation level is higher since the deposited
debris have a ground rising effect.
92
6.5. Morphological impacts
6.5.4
Morphological impacts on other locations
Early simulations performed with a mobile bed condition throughout the whole domain showed massive morphological changes on the seafronts of Caparica, Caravelos and Oeiras. However these can be
overestimated since the employed mesh was coarser, 25 [m] cells on average, and no built environment
features were at the time included in the DEM. Boundary conditions were the ones presented earlier on
this chapter. The water level that follows the early overland propagation is higher in the debris flow case,
as debris deposition is, virtually, an increasing ground level.
A clear pattern of deposition-erosion-deposition is noticeable on Figures 6.39a and 6.39c. As the run-up
progressed to higher grounds its momentum and transport capacity were decreasing and the material,
firstly dragged from the riverfront, settled much further from there. When the water recedes back to
the ocean its velocity increases and it incorporates debris again, dragging them towards the waterfront.
During the last stages of draw-down, when velocities are lower and transport capacity decreases, new
deposits of debris occur by the waterfront.
(a) Morphological changes at Caparica - Trafaria seafront
(b) Morphological changes at Carcavelos
(c) Morphological changes at Oeiras
Bed variation (m)
-1.0
-0.5
0.0
0.5
1.0
Figure 6.39: Morphological impacts on remaining locations
93
6. A 1755 tsunami in today’s Tagus estuary topography
94
Chapter 7
Conclusions and recommendations
7.1
Main conclusions
The main objective of the dissertation was to generate numerical data on the inundation and run-up of
designated areas in the the Tagus estuary.
To achieve this goal three main tasks were undertaken. In general terms:
1. The conceptual model of STAV2D was reviewed and complemented with appropriate boundary
conditions. In particular, a new boundary condition scheme had to be implemented. The RankineHugoniot algorithm, which was developed as a secondary objective of the present dissertation,
added significant flexibility to the CFD model. The performed functional tests attained very satisfactory results, with particular emphasis to flow reversal handling, as this feature is a significant
enhancement in the applicability range of the STAV2D model.
2. The geometrical description of the estuary waterfront was considerably enhanced. A very considerable effort was employed to achieve a reality-like discretization of the city’s buildings. These
were carefully modelled to ensure that the streets would behave as small canals for the overland
propagating flow. Considerable effort was also placed in articulating STAV2D with geographical
information tools. The Geographical Information Systems toolkit provides useful tools for pre and
post processing operations. It provides a much required compatibility platform, linking simulation
outputs and geographical information with few user steps. The overall effort allowed for numerical
simulations with a spatial discretization of 2.5 [m] in an exact digital model of the built environment.
3. Numerical simulations were performed to generate a database of water depth and velocity at several
locations in the Tagus estuary. Initial and boundary conditions were built from bibliography (ocean
boundary) and hydrologic analysis (Tagus boundary). Special emphasis was placed in the characterization of the inundation and run-up on the Alcântara and Downtown waterfronts, in Lisbon,
Almada waterfront (Cacilhas) and Vale do Zebro (Barreiro). At these locations, the database allows
for detailed estimates of wave arrival, its hydrodynamic properties and debris deposition patterns.
95
7. Conclusions and recommendations
The conceptual model that provides the theoretical basis for this dissertation adds significant value to
tsunami modelling efforts, in the sense that it allows for a more quantitative assessment of the impact of
debris transport by the wave. In fact, from the simulations, it is clear that the impacts of the tsunami
propagation are different depending on whether it is considered a debris flow or a clear-water wave front.
The results obtained and discussed in chapter 6 for maximum run-up levels and extents are consistent
with previous works, namely thos published in Baptista et al. (2006), with an average run-up height
of around 5 [m] at the waterfronts of Alcântara and Downtown. Nonetheless, the governing equations
feature higher detail concerning debris transport, rheological behaviour changes and flow resistance laws.
When compared to clear water results the debris flow had higher destructive capacity, which highlights
the purpose of the present work.
When using meshes with such detail, numerical stability issues often arise. Even so,the numerical model
successfully performed throughout all scenario attempts, with an average performance of 1:30, which is
to say, CPU time was 30 times the simulation time.
The outcome of these simulations is believed to be particularly useful for the Civil Protection services.
Detailed data is provided in what concerns arrival times, run-up depths and velocities and debis scattering
that may lead to unusable areas on the aftermath of the tsunami. This information can be useful in the
design of evacuation plans, victim support and rescue operations and even simulacra training.
7.2
Directions for future work
The present work provides the primary basis to successfully converge academic research with public
service in what concerns to this kind of natural hazard. In fact, it may even extend to other phenomena
within the fluid dynamics field. A solid interaction between public safety services and research teams can
greatly contribute to better response mechanisms and, ultimately, to human life protection.
Emphasizing on tsunami modelling, lagrangian methods, such as Smoothed Parcticle Hydrodynamics,
are especially fitted for describing the overland propagation stage of tsunamis (Dias, 2006). Their ability
to discretely distinguish between solid and liquid phases make this technique specially appealing for this
kind of numerical study.
Lastly, further improvements can be made on the built environment discretization. Recent methods
for ultrasonic scanning of city streets would provide a very accurate topographic basis for this type of
numerical work. Fine scale detail that is otherwise lost could improve result reliability in particularly
complex locations if such methods were available
This work provided a new contribution for the forecast of the impacts of a large tsunami in the Tagus
estuary, by employing ideas and methods not known to have been previously explored. It is believed that
this research line can be further explored, in conjunction with Civil Protection competent authorities, in
order to promote the safety of Lisbon’s citizens if the worst - a tsunami event - is to happen.
96
Bibliography
ABC News (2012). Japan earthquake: before and after.
Abe, K. (1989). Size of Great Earthquakes of 1837 - 1974 inferred from tsunami data. J. Geophysics.
Aris, R. (1989). Vectors, Tensors and the Basic Equations of Fluid Mechanics. Prentice-Hall.
Atkin, R.J. & Crane, R.E. (1976). Continuum theories of mixtures: Basic theory and historical
development. Q. Jl. Mech. appl. Math., XXIX, 209–244.
Baptista, M.A. & Miranda, J.M. (2009). Revision of the Portuguese catalog of tsunamis. Natural
Hazards and Earth System Sciences.
Baptista, M.A., Miranda, J.M., Miranda, P.M. & Victor, L.M. (1998a). The 1775 Lisbon
Tsunami: Evaluation of the Tsunami Parameters. J. Geodynamics.
Baptista, M.A., Miranda, J.M., Miranda, P.M. & Victor, L.M. (1998b). Contrains on the Source
of the 1755 Lisbon Tsunami Inferred from Numerical Modelling of Historical Data on the Source of the
1755 Lisbon Tsunami . J. Geodynamics.
Baptista, M.A., Miranda, J.M. & Luis, J.F. (2006). Tsunami Propagation along Tagus Estuary
(Lisbon, Portugal) Preliminary Results. Science of Tsunami Hazards.
Baptista, M.A., Miranda, J.M., Omira, R. & Antunes, C. (2011). Potential inundation of Lisbon
downtown by a 1755-like tsunami . Natural Hazards and Earth System Sciences.
Batchelor, G.K. (1967). An Introduction to Fluid Dynamics. Cambridge University Press.
Ben-Mehanem, A., Singh, S.J. & Solomon, F. (1970a). Deformation of an homogeneous earth model
finite by dislocations. Reviews of Geophysical and Space Physics.
Ben-Mehanem, A., Singh, S.J. & Solomon, F. (1970b). Static deformation of a spherical earth model
by internal dislocations. Reviews of Geophysical and Space Physics.
Bento, R., Cardoso, R. & Lopes, M. (2004). Earthquake resistant structures of portuguese old
‘pombalino’ buildings. 13th World Conference on Earthquake Engineering, Paper No. 918.
Biggs, S. & Sheldrick, A. (2011). Tsunami Slams Japan After Record Earthquake, Killing Hundreds.
Bloomberg.com.
97
British Accounts (1755). The Lisbon Earthquake of 1755. Anonymous Letter, November 18th 1755 .
The British Society in Portugal.
Canelas, R. (2010). 2D Mathematical Modeling of Discontinuous Shallow Sediment-laden Flows. Master’s thesis, Instituto Superior Técnico.
Chanson, H. (2005). Le Tsunami du 26 décembre 2004 : un phénomène hydraulique d’ampleur internationale. La Houille Blanche.
Dias, F. (2006). Dynamics of tsunami waves. NATO Advanced Research Workshop.
Ferreira, R.M.L. (2005). River Morphodynamics and Sediment Transport - Conceptual Model and
Solutions. Ph.D. thesis, Instituto Superior Técnico.
Ferreira, R.M.L., Franca, M., Leal, J. & Cardoso, A.H. (2009). Mathematical modelling of
shallow flows: Closure models drawn from grain-scale mechanics of sediment transport and flow hydrodynamics. Canadian Journal of Civil Engineering.
Fine, I., Rabinovichb, A., Bornholdd, B., Thomsonb, R. & Kulikovb, E. (2006). Tsunami
Modeling Manual - TUNAMI Model . Tohoku University.
Foley, J.D., van Dam, A., Feiner, S.K. & Hughes, J.F. (1996). Computer Graphics Principles and
Practice. Addison-Wesley.
Geuzaine, C. & Remacle, J.F. (2009). Gmsh: a three-dimensional finite element mesh generator with
built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering,
79, 1309–1331.
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.
Helal, M. (2008). Tsunamis from nature to physics. Elsevier.
Hircsh, C. (1988). Numerical Computation of Internal and External Flows: Fundamentals of Numerical
Discretization. Wiley.
Imamura, F. & Shuto, N. (1990). Tsunami propagation by use of numerical dispersion. International
Symposium on Computational Fluid Dynamics.
Imamura, F., Yalciner, A. & Ozyurt, G. (2006). Tsunami Modeling Manual - TUNAMI Model .
Tohoku University.
Jao, D. (Version 5). Ferrari-Cardano’s derivation of the quartic formula. PlanetMath.org.
Jao, D. (Version 8). Cardano’s derivation of the cubic formula. PlanetMath.org.
Jimenez, J.A. & Madsen, O.S. (2003). A simple formula to estimate settling velocity of natural sediments. J. Waterway, Port, Coast. and Oc. Engrd.
Kundu, P., Cohen, I. & Downling, D. (1990). Fluid Mechanics. Academic Press.
98
Kurganov, A. & Petrova, G. (2007). A second-order wellbalanced positivity preserving central-upwind
scheme for the saint-venant system. Communications in Mathematical Sciences.
LeVeque, R.J. (1990). Numerical Methods for Conservation Laws. Birkhauser.
LeVeque, R.J. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge University Press.
Luis, J.F. (2007). Mirone: A multi-purpose tool for exploring grid data. Computers and Geosciences.
Álvaro S. Pereira (2006). The Opportunity of a Disaster: The Economic Impact of the 1755 Lisbon
Earthquake. Centre for Historical Economics and Related Research at York University.
Maeda, T., Furumura, T., Sakai, S. & Shinohara, M. (2011). Significant tsunami observed at
ocean-bottom pressure gauges during the 2011 off the Pacific coast of Tohoku Earthquake. Earth, Planets
and Space.
Masterlark, T. (2003). Finite element model predictions of static deformation from dislocation sources
in a subduction zone: Sensivities to homogeneous, isotropic, poisson-solid, and half-space assumptions.
Journal of Geophysical Research.
McGinley, J. (1969). A comparison of observed permanent tilts and strains due to earthquakes with those
calculated from displacement dislocations in elastic earth models. California Institute of Technology,
Pasadena.
Mellor, G.L. (1998). Introduction to Physical Oceanography. Springer.
Moreira de Mendonça (1758). Historia Universal dos Terramotos que tem havido no mundo de que
ha noticia desde a sua creção até ao século presente. Biblioteca Nacional de Lisboa, Portugal.
Murillo, J. (2006). Two-Dimentional Finite Volume Numerical Models for Unsteady Free Surface Flows.
University of Zaragoza.
Murillo, J. & García-Navarro, P. (2010). Weak solutions for partial differential equations with
source terms: application to the shallow water equations. Journal of Computational Physics.
Pope, S. (2000). Turbulent Flows. Cambridge University Press.
Rabinovich & Thomson, R.E. (2007). The 26 December 2004 Sumatra Tsunami: Analysis of Tide
Gauge Data from the World Ocean. Pure and Applied Geophysics.
Rankine, W.M. (1869). On the thermodynamic theory of waves of finite longitudinal disturbance. Philosophical Transactions of the Royal Society of London.
Roe, P. (1981). Approximate Riemann solvers, parameter vectors and difference schemes. Journal of
Computational Physics.
Steketee, J.A. (1958). On volterra’s dislocation in a semi-infinite elastic medium. Canadian Journal
of Physics.
Stokes, G. (1847). On the theory of oscillatory waves. Transactions of the Cambridge Philosophical
Society.
99
Takahashi, T. (2007). Debris Flow, Mechanics, Prediction and Countermeasures. Taylor and Francis.
Titov, V., Rabinovich, A.B., Mofjeld, H.O., Thomson, R.E. & González, F.I. (2005). The
Global Reach of the 26 December 2004 Sumatra Tsunami . U.S. Department of Commerce Publications.
Toro, E. (2001). Shock-Capturing Methods for Free Surface Shallow Flows. Wiley.
Toro, E. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer.
Urban, S. (1755). The Gentleman’s Magazine for the year 1755 . R. St John’s Gate - London.
Volterra, V. (1907). Sur l’équilibre des corps élastiques multiplement connexes. Annales Scientifiques
de l’Ecole Normale Supérieure.
Vázquez-Cedón, M. (1999). Improved treatment of source terms in upwind schemes for the shallow
water equations with irregular geometry. Journal of Computational Physics.
Wang, X. & Liu, P.L.F. (2006). Comcot user manual – version 1.6.
Wesseling, P. (2009). Principles of Computational Fluid Dynamics. Springer SCM.
White, F.M. (2002). Fluid Mechanics. McGraw-Hill.
100
Appendices
Appendix 1: Wave dynamics
A brief description of the classic mathematical description of waves is herein presented. The most basic
assumptions are that water is incompressible, and has a constant density value, ρ, and that dissipation
due to bottom friction can be neglected for this type of theoretical analysis as it is often included as a
source term in most models (Dias, 2006).
The irrotational flow of an incompressible fluid is governed by the conservation of mass,
∇2D ~u +
∂w
=0
∂z
(1)
and the conservation of momentum,
ρw
D~u
= −∇2D p~
Dt
;
ρw
Dw
∂p
=−
− ρw g
Dt
∂z
(2)
where ~u is the horizontal velocity (u, v), w is the vertical velocity, and g is the gravitational acceleration.
A more comprehensive version of the momentum equations will be subjected to an in-depth analysis in
Chapters 3.1 and 3. From the definition of From the definition of irrotacionality,
∇ × (u, v, w) = 0
(3)
and thus there exists a potential function φ(x, y, z, t) such that,
∇2D φ = ~u
;
∂φ
=w
∂z
(4)
Use of the potential function is particularly convenient when dealing with a free surface boundary condition (Mellor, 1998). The continuity equation now becomes the Laplace’s equation,
∇22D φ +
∂2φ
=0
∂z 2
(5)
101
and the momentum equation can be written in terms of φ if integrated over z, yielding Bernoulli’s
equation,
∂φ 1
1
+ |∇2D φ|2 +
∂t
2
2
∂φ
∂z
2
+ gz +
p
=0
ρw
(6)
which are directly used in solving flows bounded with solid surfaces. Since the velocity normal to a
surface is null, the value of w = ∂φ/∂z = 0 must be null at the bottom.
The surface boundary condition, which is moving, is far more complicated since two boundary conditions
are required. The free surface is represented by F (x, y, z, t) = η(x, y, t) − z = 0 and the bottom is given
by z = −h(x, y). The necessary boundary conditions can be derived from either a kinematic or dynamic
point of view. The free-surface kinematic boundary condition (FSKBC) is given by,
∂η
∂φ
DF
=
+ ∇2D φ · ∇2D η −
=0
Dt
∂t
∂z
(7)
and the dynamic boundary condition at the free-surface (FSDBC) and at the bottom (BDBC) are given
by Bernoulli’s equation at z = η(x, y, t) and z = −h(x, y), respectively. The FSDBC is given by,
∂φ 1
1
+ |∇2D φ|2 +
∂t
2
2
∂φ
∂z
2
+ gη = 0
(8)
and the BDBC is expressed as,
∇2D φ · ∇2D h +
∂φ
=0
∂z
(9)
The goal is to solve the Laplace’s equation algonside with the FSKBC, FSDBC and BDBC equations. The
conservation of momentum equation is not required in the solution procedure and it is used a posteriori
to find the pressure p once η and φ have been found (Dias, 2006).
102
Appendix 2: Conservation of energy
The energy conservation equation is herein derived with the single purpose of completing the set of conservation laws in fluid dynamics. This equation is not implemented in the numerical model as it becomes
redundant when used alongside with both mass and momentum conservation equations. Furthermore this
derivation applies only to the mechanical energy as no thermodynamic features are presently implemented
in the numerical model.
The balance between local inertia and applied external forces expressed in the momentum conservation
law, equation (3.17), may be converted into an energy conservation law by taking the corresponding work
of those forces. One might do so by working with the rate of work, defined as,
d
dxi
Wk = Fki
= Fki ui
dt
dt
(10)
where Wk is the work of force F~k along the ~x direction. This suggests that multiplying equation (3.17)
by ui , and summing over i, will yield the energy conservation equation (Kundu et al., 1990). Thus, after
some algebraic manipulations, one should obtain the following result,
D
ρ
Dt
1 2
u
2 i
+ u2i
∂p
∂τij
Dρ
= −ui
+ ui ρgi + ui
Dt
∂xj
∂xj
(11)
which is the energy conservation equation, although only suitable to incompressible flows.
103
104
Appendix 3: Implemented boundary condition algorithm
In the interest of future reference this subsection summarizes the implemented subroutines. The complete
Fortran90 code was replaced by an equivalent pseudo-code.
NOTE: During debug tests it was noticed that for very small ∆h, for instance ∆h < 0.1[mm], the
characteristics λR and λL had virtually the same values and therefore a consistent entropy condition
could not be verified for either shocks or expansion waves. In order to overcome this drawback the
following threshold was implemented: if ∆h < 0.1[mm] then uL = uR and hL = hR , that is to say, the
boundary replicates the values of the inner cells from time ti to ti+1 .
Implemented subroutine for INLETS
Acquire ql from q(t) hydrograph
Acquire hr and ur from the inner cells
if ( abs(hr - hl) < 1E-4) then
!Avoid numerical residue
hr = hl
else if (hl > hr) then
!Shock Condition
Solve equation (3.14), 4th order (Ferrari-Cardano’s Method)
Obtain 4 possible solutions for hl
Reject complex solutions (Im > 1E6)
Reject negative real solutions (hl <= 0)
Test the entropy condition for the remaining solution/s
Use the only validated solution
else if (hl < hr) then
!Expansion Wave Condition
Solve equation (3.15), 3rd order (Cardano’s Formula)
Obtain 3 possible solutions for hl
Reject complex solutions (Im > 1E6)
Reject negative real solutions (hl <= 0)
Test the entropy condition for the remaining solution/s
Use the only validated solution
end if
Return hl
Compute ul = ql/hl
Assign the (u,v) velocity vector accordingly with the edge normal
105
Implemented subroutine for OUTLETS
Acquire hl from h(t) hydrograph
Acquire hr and ur from the inner cells
if ( abs(hr - hl) < 1E-4) then
!Avoid numerical residue
hr = hl
else if (hl > hr) then
!Shock Condition
Solve equation (3.16), 2nd order (Quadratic Formula)
Obtain 2 possible solutions for ul
Reject complex solutions (Im > 1E6)
Test the entropy condition for the remaining solution/s
Use the only validated solution
else if (hl < hr) then
!Expansion Wave Condition
Solve equation (3.17), 3rd order (Cardano’s Formula)
Obtain 2 possible solutions for ul
Reject complex solutions (Im > 1E6)
Test the entropy condition for the remaining solution/s
Use the only validated solution
end if
Return ul
Assign the (u,v) velocity vector accordingly with the edge normal
106
Appendix 3: Stratification Source Terms
The 2nd order perturbations in the shallow-water equations induced by stratification are parametrized in
a source term represented by the symbol S(U ), defined as,
T
S(U ) = [ 0 Sx Sy 0 ] ;
(12)
where,
Sx = −
−
−
−
−
Sy = −
−
−
−
−
1
uh∂t ρm + u2 h∂x ρm + uvh∂y ρm
ρm
1
∂t (−ρm uh + ρ2 u2 h2 + ρ1 u1 h1 )
ρm
1
1
1
1
2
2
2
2
2
2
∂x −ρm u h + ρ2 u2 h2 + ρ1 u1 h1 − gρm h + gρc h2 + gρ1 h1 + gρ1 h1 h2
ρm
2
2
2
1
∂y (−ρm uvh + ρ2 u2 v2 h2 + ρ1 u1 v1 h1 )
ρm
1
(−ghρm ∂x Zb + (ρ2 h2 + ρ1 h1 ) g∂x (Zb ))
ρm
(13)
1
vh∂t ρm + uvh∂x ρm + v 2 h∂y ρm
ρm
1
∂t (−ρm vh + ρ2 v2 h2 + ρ1 v1 h1 )
ρm
1
∂x (−ρm uvh + ρ2 u2 v2 h2 + ρ1 u1 v1 h1 )
ρm
1
1
1
1
∂y −ρm v 2 h + ρ2 v22 h2 + ρ1 v12 h1 − gρm h2 + gρc h22 + gρ1 h21 + gρ1 h1 h2
ρm
2
2
2
1
(−ghρm ∂y Zb + (ρc h2 + ρ1 h1 ) g∂y (Zb ))
ρm
(14)
where ρm , u and v are the depth-averaged density and velocity components of the mixture, ρ1 and ρ2
are depth averaged densities of layers 1 and 2, respectively, h1 and h2 are the depths of layers 1 and 2,
respectively, u1 and u2 are the depth-averaged velocities in the x direction of layers 1 and 2, respectively
and v1 and v2 are the depth-averaged velocities in the x direction of layers 1 and 2, respectively
107