# Title: The 2004 Indian Ocean tsunami # # Description: # # The 2004 Indian Ocean tsunami was caused by a large-scale fault # rupture ($> 1000$ km) at the Indian--Australian and # Eurasian--Andaman plate boundaries. This example uses the fault # model of Grilli et al. \cite{grilli2007} as initial conditions to a # Saint-Venant solution of the subsequent tsunami. The evolution of # the wave elevation is illustrated in the animation of Figure # \ref{h}.a. Adaptivity is used to track the various wave fronts # (Figure and animation \ref{h}.b) and the terrain is reconstructed # dynamically based on the # \htmladdnormallinkfoot{ETOPO1}{http://www.ngdc.noaa.gov/mgg/global/global.html} # dataset. # # \begin{figure}[htbp] # \caption{\label{h}(a) Animation of wave elevation. Dark red is # larger than 2 m, dark red smaller than -2 m. (b) Animation of # adaptive level of refinement. Dark red is 0.8 nautical miles, dark # blue is 101 nautical miles.} # \begin{center} # \begin{tabular}{cc} # \htmladdnormallinkfoot{\includegraphics[width=0.5\hsize]{h.eps}}{h.mpg} & # \htmladdnormallinkfoot{\includegraphics[width=0.5\hsize]{level.eps}}{level.mpg} \\ # (a) & (b) # \end{tabular} # \end{center} # \end{figure} # # The maximum wave elevation reached over 10 hours after fault rupture # is illustrated on Figure \ref{elevation}. # # \begin{figure}[htbp] # \caption{\label{elevation}Maximum wave elevation reached over 10 # hours, 1 metre interval contours. (a) Bay of Bengal, dark blue zero, # dark red $>5$ m. (b) Detail near Northern Sumatra and Thailand, dark # blue zero, dark red $> 8$ m.} # \begin{center} # \begin{tabular}{cc} # \includegraphics[width=0.5\hsize]{hmax.eps} & # \includegraphics[width=0.5\hsize]{hmax-detail.eps} \\ # (a) & (b) # \end{tabular} # \end{center} # \end{figure} # # Finally Figure \ref{data} gives a comparison of the observed (using # tide gauges) and modelled timeseries at specific locations in the # Indian Ocean. Figure \ref{jason} gives a comparison of the modelled wave # profile with the profile observed by the Jason-1 satellite altimeter. # # \begin{figure}[htbp] # \caption{\label{data}Observed and modelled time-series of wave # elevation (metres) at different tide gauge locations. Horizontal # axis is time in hours after the initial fault rupture.} # \begin{center} # \begin{tabular}{cc} # \includegraphics[width=0.5\hsize]{hani.eps} & # \includegraphics[width=0.5\hsize]{dieg.eps} \\ # Hannimaadhoo, Maldives & Diego Garcia \\ # \includegraphics[width=0.5\hsize]{male.eps} & # \includegraphics[width=0.5\hsize]{colo.eps} \\ # Male, Maldives & Columbo, Sri Lanka \\ # \includegraphics[width=0.5\hsize]{gana.eps} \\ # Gan, Maldives # \end{tabular} # \end{center} # \end{figure} # # \begin{figure}[htbp] # \caption{\label{jason}Observed (Jason-1 satellite altimeter) and modelled wave profiles.} # \begin{center} # \includegraphics[width=0.8\hsize]{jason.eps} # \end{center} # \end{figure} # # More details on this simulation and the method used is given in # \cite{popinet2011}. # # Author: St\'ephane Popinet # Command: gerris2D -m tsunami.gfs | gfsview2D tsunami.gfv # Required files: output.gfs tsunami.gfv h.gfv level.gfv hmax.gfv hmax-detail.gfv hanires.txt diegres.txt maleres.txt colores.txt ganares.txt jason.xy jasonres.txt # Version: 110107 # Running time: # Generated files: h.eps level.eps h.mpg level.mpg hmax.eps hmax-detail.eps hani.eps dieg.eps male.eps colo.eps gana.eps jason.eps # the radius of the Earth (in meters) Define RADIUS 6371220. # the size of the domain Define LENGTH 6000e3 # the domain is centered on 94E, 8N Define LONGITUDE 94. Define LATITUDE 8. # anywhere with a fluid depth < 1 cm is considered dry Define DRY 0.01 # 12 levels of refinement max i.e. a maximum resolution of 6000e3/2^12 = 1464 m Define MAXLEVEL 12 # keep a coarse band, 0.04 wide, on all boundaries of the domain # these act as "sponges" before the waves exit the domain Define LEVEL (fabs (rx) < 0.46 && fabs (ry) < 0.46 ? MAXLEVEL : 5) 1 0 GfsRiver GfsBox GfsGEdge {} { PhysicalParams { # length of the domain (m) L = LENGTH # gravity is 9.81 m/s^2 g = 9.81 # from now on, units have been chosen to be metres and seconds } # use the longitude/latitude metric. x and y coordinates are now # longitude and latitude in degrees MetricLonLat M RADIUS # shift the origin to where we want it MapTransform { tx = LONGITUDE ty = LATITUDE } # 10 hours Time { end = 36000 } # use the terrain module for topography GModule terrain # use the okada module for initial deformations GModule okada Refine 5 # maintain Zb as the terrain elevation defined by the samples in # ETOPO1 # the 'etopo1' database needs to have been created beforehand and # be accessible within GFS_TERRAIN_PATH, see the 'terrain' module # documentation above VariableTerrain Zb { basename = etopo1 } { # preserve lake-at-rest balance (but not volume) when # reconstructing bathymetry reconstruct = 1 } # the default minmod limiter is too dissipative for tsunamis AdvectionParams { gradient = gfs_center_sweby_gradient } # deformation from Grilli et al, 2007, Table 1 # 5 segments triggered at different times # Segment 1 Init {} { D = 0 } InitOkada D { x = 94.57 y = 3.83 depth = 11.4857e3 strike = 323 dip = 12 rake = 90 length = 220e3 width = 130e3 U = 18 } # Initial water level is at z = D Init { start = 0 } { P = MAX (0., D - Zb) } # Segment 2 EventList { start = 212 step = 6 end = 272 } { Init {} { D = 0 } InitOkada D { x = 93.90 y = 5.22 depth = 11.4857e3 strike = 348 dip = 12 rake = 90 length = 150e3 width = 130e3 U = 23 } } # make sure the deformation is well resolved AdaptGradient { start = 212 istep = 1 end = 272 } { cmax = 0.05 cfactor = 2 minlevel = 5 maxlevel = LEVEL } D # add it to the current elevation field (only if wet) Init { start = 272 } { P = (P < DRY ? P : MAX (0., P + D)) } # Segment 3 EventList { start = 528 step = 6 end = 588 } { Init {} { D = 0 } InitOkada D { x = 93.21 y = 7.41 depth = 12.525e3 strike = 338 dip = 12 rake = 90 length = 390e3 width = 120e3 U = 12 } } # make sure the deformation is well resolved AdaptGradient { start = 528 istep = 1 end = 588 } { cmax = 0.05 cfactor = 2 minlevel = 5 maxlevel = LEVEL } D # add it to the current elevation field (only if wet) Init { start = 588 } { P = (P < DRY ? P : MAX (0., P + D)) } # Segment 4 EventList { start = 853 step = 6 end = 913 } { Init {} { D = 0 } InitOkada D { x = 92.60 y = 9.70 depth = 15.12419e3 strike = 356 dip = 12 rake = 90 length = 150e3 width = 95e3 U = 12 } } # make sure the deformation is well resolved AdaptGradient { start = 853 istep = 1 end = 913 } { cmax = 0.05 cfactor = 2 minlevel = 5 maxlevel = LEVEL } D # add it to the current elevation field (only if wet) Init { start = 913 } { P = (P < DRY ? P : MAX (0., P + D)) } # Segment 5 EventList { start = 1213 step = 6 end = 1273 } { Init {} { D = 0 } InitOkada D { x = 92.87 y = 11.70 depth = 15.12419e3 strike = 10 dip = 12 rake = 90 length = 350e3 width = 95e3 U = 12 } } # make sure the deformation is well resolved AdaptGradient { start = 1213 istep = 1 end = 1273 } { cmax = 0.05 cfactor = 2 minlevel = 5 maxlevel = LEVEL } D # add it to the current elevation field (only if wet) Init { start = 1273 } { P = (P < DRY ? P : MAX (0., P + D)) } # we don't want the arrival time to be interpolated from coarse # to fine, so we make it a "boolean" variable VariableBoolean Atime Init {} { P = MAX (0., D - Zb) Pmax = 0. Atime = -1. } Init { istep = 1 } { # elevation of the wet surface Hwet = (P > DRY ? H : NODATA) # maximum amplitude Pmax = (P > DRY && H > Pmax ? H : Pmax) # arrival time for an amplitude of +- 5 cm Atime = (Atime >= 0. ? Atime : P > DRY && fabs (H) > 5e-2 ? t : -1.) # implicit scheme for quadratic bottom friction with coefficient 1e-3 U = P > DRY ? U/(1. + dt*Velocity*1e-3/P) : 0 V = P > DRY ? V/(1. + dt*Velocity*1e-3/P) : 0 } # adapt on gradient of wet surface elevation # with a tolerance of 5 cm AdaptGradient { istep = 1 } { cmax = 0.05 cfactor = 2 minlevel = 5 maxlevel = LEVEL } (P < DRY ? 0. : P + Zb) # also keep enough resolution to resolve the maximum elevation # field with a 5 cm discretisation error max AdaptError { istep = 1 } { cmax = 0.05 minlevel = 5 maxlevel = LEVEL } Pmax Global { #define RESOLUTION (LENGTH/pow (2, MAXLEVEL)) #define close_enough(x,y) (distance(x,y,0.) < RESOLUTION) } # include a list of locations for which to output timeseries Include output.gfs # Jason profiles OutputLocation { start = 6900 } jason.txt jason.xy OutputTime { istep = 1 } stderr OutputBalance { istep = 1 } stderr OutputTiming { istep = 100 } stderr # save a snapshot every 100 iterations EventScript { istep = 100 } { rm -f snapshot-*.gfs } OutputSimulation { istep = 100 } snapshot-%ld.gfs # output solution on standard output every 20 timesteps # for on-the-fly visualisation with GfsView OutputSimulation { istep = 20 } stdout # ouput solution every 900 seconds (15 minutes) OutputSimulation { step = 900 } sim-%g.gfs EventScript { step = 900} { gzip -f sim-*.gfs } OutputScalarStats { istep = 1 } p { v = P } OutputScalarNorm { istep = 1 } u { v = Velocity } OutputScalarNorm { istep = 1 } U { v = U } OutputScalarNorm { istep = 1 } V { v = V } OutputScalarNorm { istep = 1 } hwet { v = Hwet } # animation of the wave elevation GModule gfsview OutputView { step = 60 } { ppm2mpeg > h.mpg } { width = 800 height = 700 } h.gfv OutputView { start = 7200 } { convert ppm:- eps2:h.eps } { width = 800 height = 700 } h.gfv # animation of the level of refinement OutputView { step = 60 } { ppm2mpeg > level.mpg } { width = 800 height = 700 } level.gfv OutputView { start = 7200 } { convert ppm:- eps2:level.eps } { width = 800 height = 700 } level.gfv # animation of the velocity OutputPPM { step = 60 } { ppm2mpeg > velocity.mpg } { maxlevel = 10 v = Velocity min = 0 max = 0.25 } # animation of the topography OutputPPM { step = 60 } { ppm2mpeg > zb.mpg } { maxlevel = 10 v = Zb } # graphics OutputView { start = end } { convert ppm:- eps2:hmax.eps } { width = 1600 height = 1600 } hmax.gfv OutputView { start = end } { convert ppm:- eps2:hmax-detail.eps } { width = 1600 height = 1600 } hmax-detail.gfv EventScript { start = end } { # timeseries at specific locations gnuplot < 19.2 ? 0 : \$8) w l t 'modelled' EOF } } { dry = DRY } GfsBox { # inflow/outflow on all boundaries with a reference water level at # z = 0 top = Boundary { BcSubcritical V -Zb } bottom = Boundary { BcSubcritical V -Zb } right = Boundary { BcSubcritical U -Zb } left = Boundary { BcSubcritical U -Zb } }