Modeling, Simulating,
and Analyzing a Tsunami

We will use Mathematica to study the propagation of a slightly simplified tsunami wave, including:

• Building the Model

• Simulating the Model

• Analyzing the Simulation Results

In this application having more memory means you can simulate a more realistic model, in other words use a model that can resolve fine features, and get more accurate analysis results.

Start the Mathematica kernel. This will also ask whether you want to run, initialization cells. Answer yes to this question.

$Version

If you are not familiar with Mathematica, this section provides a short introduction to help get you started. Mathematica provides a unique interface allowing you to combine presentation and calculations all in one document, called a notebook.

The document you are reading is a Mathematica notebook and is formatted as a slideshow-style presentation, complete with navigation buttons at the top of each page. Some sections can be expanded by clicking on the small triangle preceding them, such as this section.  Calculations are performed in cells with a yellow background. You will see more of these in the following slides.

Most of the cells here are "initialization cells" that will automatically evaluate themselves as soon as any calculation is started.  This guarantees that all of the necessary definitions are loaded.  You will be asked to evaluate the necessary cells when the time comes.  This is done by placing the cursor in the appropriate cell and then typing ⇧+⌤.

This section defines a few utilities that will make the demo self-contained in the presentation notebook. This section is not essential reading for following the demo, but contains initialization code necessary to make the demo work.

This provides utilities to compute distances on the earth, given latitude and longitude information.

Needs["Utilities`Notation`"]

Off[General::spell];  Off[General::spell1];

BeginPackage["Vincenty`"];

VincentyDistance::usage="VincentyDistance[{lat1, lon1},{lat2, lon2}] computes the distance between two points with the given latitudes and longitudes.";

Begin["`Private`"];

SetAttributes[ArcTanOrAngle, Listable];
ArcTanOrAngle[tan_, ϕ_] := 
    If[Abs[ϕ] < π/2, ArcTan[tan], ϕ];

VincentyDistanceS[{ϕ1_, θ1_},{ϕ1_, θ1_}]  := 
    {0., {0.,0.}};

VincentyDistanceS[{ϕ1_, θ1_}, {ϕ2_, θ2_}] := 
  Module[{a0, flat, r, b0, tanu1, tanu2, dtmp, cosu1, sinu1, cosu2, 
      sinu2, ω, cs, ss, cosalpha2, c2sm, u2, a, b, dsigma, s, alpha12, 
      alpha21},
    a0 = 6378137.00;
    flat = 1./298.257223563;
    r = 1. - flat;
    b0 = a0*r;
    
    tanu1 = r*Tan[ϕ1];
    tanu2 = r*Tan[ϕ2];
    
    dtmp = ArcTanOrAngle[tanu1, ϕ1];
    cosu1 = Cos[dtmp];
    sinu1 = Sin[dtmp];
    
    dtmp = ArcTanOrAngle[tanu2, ϕ2];
    cosu2 = Cos[dtmp];
    sinu2 = Sin[dtmp];
    
    ω = θ2 - θ1;
    
    itfun[λ_] := Module[{ss1, ss2, tansigma, sinalpha, c},
        ss1 = cosu2*Sin[λ];
        ss2 = cosu1*sinu2 - sinu1*cosu2*Cos[λ];
        ss = Sqrt[ss1*ss1 + ss2*ss2];
        cs = sinu1*sinu2 + cosu1*cosu2*Cos[λ];
        tansigma = ss/cs;
        sinalpha = cosu1*cosu2*Sin[λ]/ss;
        cosalpha2 = Cos[ArcSin[sinalpha]]^2;
        c2sm = cs - 2.*sinu1*sinu2/cosalpha2;
        c = flat/16.*cosalpha2*(4. + flat*(4. - 3.*cosalpha2));
        ω + (1. - c)*flat*
            sinalpha*(ArcSin[ss] + c*ss*(c2sm + c*cs*(-1. + 2.*c2sm*c2sm)))];
    
    λ = FixedPoint[itfun, ω];

    u2 = cosalpha2*(a0*a0 - b0*b0)/(b0*b0); 
    a = 1. + (u2/16384.)*(4096. + u2*(-768. + u2*(320. - 175.*u2))); 
    b = (u2/1024.)*(256. + u2*(-128. + u2*(74. - 47.*u2)));
    
    dsigma = 
      b*ss*(c2sm + (b/
                  4.)*(cs*(-1. + 2.*c2sm*c2sm) - (b/6.)*
                    c2sm*(-3. + 4.*ss*ss)*(-3. + 4.*c2sm*c2sm)));
    
    s = b0*a*(ArcSin[ss] - dsigma);
    
    alpha12 = 
      ArcTan[cosu1*sinu2 - sinu1*cosu2*Cos[λ], cosu2*Sin[λ]]; 
    alpha21 = 
      ArcTan[-sinu1*cosu2 + cosu1*sinu2*Cos[λ], 
        cosu1*Sin[λ]];
    {s, {alpha12, alpha21}}
];

ToRadians[deg_] := π deg/180;

ToDegrees[rad_] := 180 rad/π;

VincentyDistance[{Lat1_, Long1_},{Lat2_, Long2_}] := 
Module[{s, az},
    {s, az} = 
      VincentyDistanceS[ToRadians[{Lat1, Long1}], ToRadians[{Lat2, Long2}]];
    {s, Mod[ToDegrees[az] + 180 {2,1},360]}];

End[];

EndPackage[];

GridSpacing will allow us to explain the actual spacing on the earth.

\(GridSpacing[nx_, \ ny_] := \ 
      Module[{Δλ\  = \ \((λ\_max - λ\_min)\)/n, λm = \((λ\_max + λ\_min)\)/
              2. , 
Δϕ\  = \ \((ϕ\_max - ϕ\_min)\)/n, ϕm = \((ϕ\_max + ϕ\_min)\)/
              2. }, 
{"LongitudeSpacing" -> 
            First[VincentyDistance[{λm, ϕm}, {λm + Δλ, ϕm}]]\ "meters", 
"LatitudeSpacing" -> 
            First[VincentyDistance[{λm, ϕm}, {λm, ϕm + Δϕ}]]\ "meters"}]\)

This is the grid spacing we have chosen for the low resolution simulation.

n  =  5*16;
nx  =  n  +  1; ny  = n  +  1;  
    nxy  =  (nx  -  1)*(ny - 1);

These are options for the NDSolve function used later.

\(\(methodOptions\  := \ 
        Sequence[\ \ DependentVariables -> {h, u, v}, \ 
          Method -> {"MethodOfLines", "SpatialDiscretization" -> 
{"TensorProductGrid", \ "MinPoints" -> {nx, ny}, \ "MaxPoints" -> {nx, 
                      ny}, \ "DifferenceOrder" -> "Pseudospectral"}, 
"DiscretizedMonitorVariables" -> True, 
"Method" -> {"EventLocator", \ 
                  Method -> "Adams", \ "Event" -> 
                    Sin[π\ t/Δt], \ "EventAction" :> action}}, \ Jacobian\  -> \ None\ ];\)\)

Modeling, Simulating, and Analyzing a Tsunami

Building a Model – Partial Differential Equation

These are the shallow water equations (SWE) used to model the wave propagation. Shallow means that the ocean is wider than it is deep. In Mathematica you simply type in the equations.

swe = {(∂_t) h[λ,  ϕ,  
                  t] + (1/(R Cos[ϕ])) ((∂_λ)(u[λ,  ϕ,  
                            t]  (h[λ,  ϕ,  t]    -  
                              b[λ,  ϕ] )) + (∂_ϕ)(v[λ,  ϕ,  
                            t]  (h[λ,  ϕ,  t]    -  
                              b[λ,  ϕ] ) 
                          Cos[ϕ])) == 
            0, (((∂_ϕ) u[λ,  ϕ,  t]  v[λ,  ϕ,  t] )/R) + ((g (∂_λ) h[λ,  ϕ,  t] )/(R Cos[ϕ])) + (∂_t) u[λ,  ϕ,  
                  t] + ((u[λ,  ϕ,  t]  (∂_λ) u[λ,  ϕ,  t] )/(R Cos[ϕ])) == f v[λ,  ϕ,  
                t] , (((∂_ϕ) v[λ,  ϕ,  t]  v[λ,  ϕ,  t] )/R) + ((g (∂_ϕ) h[λ,  ϕ,  t] )/R) + (∂_t) 
                  v[λ,  ϕ,  
                  t]  + ((u[λ,  ϕ,  t]  (∂_λ) v[λ,  ϕ,  t] )/(R Cos[ϕ])) == -f u[λ,  ϕ,  t] };

These are the parameters we will use in our model. All units are in MKS (meter-kilogram-second).

R  =  6378000;   

    ω  =  2 π/(24 3600);   

    f =  2  ω Sin[ϕ];   

    g  =  9.8;

These parameters define the range over which we will simulate the model. Basically we want to simulate for two hours are reasonable, but finite, section of the ocean.

T  =  
        2 3600 ; (* Time -- two hours in seconds *) 
(Subscript[λ, min])  = -15 °  ; (Subscript[λ, max]) =  15 ° ;  (* 
      min and max longitude *) 
(Subscript[ϕ, min]) =  -35 ° ;(Subscript[ϕ, max]) =  
        5 ° ;  (* 
      min and max latitude *)

The SWE model makes the following assumptions:

The last assumption is certainly violated when a tsunami nears land.  The U.S. National Oceanic and Atmospheric Administration (NOAA) uses a model they refer to as MOST (Method of Splitting Tsunami) that uses different computations and models for three phases; generation, propagation, and inundation.  Note: This is a simplified simulation of the propagation phase.  The NOAA tsunami modeling code may take several hours to run a realistic model of a tsunami.

Note: This demonstration is based on a general partial differential equation solver that is good for prototyping models but can never compete directly with a specialized code such as the one that NOAA uses.

Modeling, Simulating, and Analyzing a Tsunami

Building the Model – Initial Value

This function computes the shortest distance to a fault line.

ldist[p1_,  p2_][
        p_?(VectorQ[#,  NumberQ] &)]  := 
If[
        p1  ==  p2,  Norm[p  -  p1], 
 
        Module[{pp1  =  p  -  p1, pp2  =  p  -  p2,  
            p12  =  (p1  -  p2)/
                Norm[p1  -  p2]}, 
If[
            Sign[p12 . pp1]  !=  Sign[p12 . pp2],  
            Norm[pp1  -  pp1 . p12 p12],  
            Min[Norm[pp1],  Norm[pp2]]]]
]

This defines the initial displacement, centered at 0° longitude and going from 10° to 20° south latitude.

(Subscript[λ, e])  =  
      0; (Subscript[ϕ, 1]) =  -10 °; (Subscript[ϕ, 2]) =  -20 °;

    h0[λ_,  ϕ_]  =  
        Exp[-10000 ldist[{(Subscript[λ, e]),  (Subscript[ϕ, 1])},  {(Subscript[λ, e]),  (Subscript[ϕ, 2])}][{λ,  ϕ}]^2];

Here is a plot of the initial displacement.  (The longitude and latitude axis ticks are labeled in radians.)

\(\(Plot3D[
        h0[λ, \ ϕ], \ \ {λ, \(-5\)\ °, \ 
          5\ °}, {ϕ, \(-8\)\ °, \ \(-22\)\ °\ }, \ PlotRange -> All, \ PlotPoints -> 100, \ Mesh -> False, \ 
        AxesLabel -> {\ "λ", \ "ϕ", "h (m)"}]\ ;\)\)

Note: The vertical scale is 1 meter, whereas, if converted to distance, the horizontal scale is thousands of kilometers.  One fascinating thing about a tsunami is that in the open ocean, the amplitude of the wave may be quite small (1 meter is not a big wave). Its destructive power is because its energy is distributed all the way from the top to the bottom of the ocean, unlike wind driven waves that have nearly all their energy very near the water surface.

Modeling, Simulating, and Analyzing a Tsunami

Building the Model – Ocean Floor

This defines a function to simulate an undersea mountain.

Seamount[s_,  w_,  {λ_,  ϕ_}][l_,  p_]  :=  
        s Exp[-w ((l - λ)^2  +  (p - ϕ)^2)];

This defines the bottom depth to be 4000 meters (about the average depth of the Indian Ocean) with some seamounts scattered within the range of the simulation parameters.

b[λ_,  ϕ_]  :=  -4000  +  Seamount[3500,  
              1000,  {-2.5 °,  -25 °}][λ,  ϕ] +  Seamount[3000,  
              1000,  {-5 °,  -23 °}][λ,  ϕ]  +   Seamount[3000,  
              1000,  {-7.5 °,  -27 °}][λ,  ϕ] +   Seamount[3000,  
              1000,  {5  °,  -24 °}][λ,  ϕ] + Seamount[3000,  
              1000,  {5  °,  -21 °}][λ,  ϕ] + Seamount[3000,  
              1000,  {5  °,  -18 °}][λ,  ϕ] + Seamount[3000,  
              1000,  {5  °,  -15 °}][λ,  ϕ] + Seamount[3000,  
              1000,  {5  °,  -12 °}][λ,  ϕ];

Here is a plot of the bottom structure.

Plot3D[
        b[λ,  ϕ],   {λ,  (Subscript[λ, min]),  (Subscript[λ, max])}, {ϕ,  (Subscript[ϕ, min]),  (Subscript[ϕ, max])},  PlotRange -> All,  
        PlotPoints -> 200,  Mesh -> False] ;

One of the important things to simulate for tsunami propagation is how the waves are scattered around features in the ocean.  The codes that NOAA uses take into account the actual sea bottom structure, but for this demonstration, a simulated bottom structure will be considered.

Modeling, Simulating, and Analyzing a Tsunami

Simulating the Model – Solution Monitoring

We want to be able to monitor the solution as it progresses ... so we add a little extra code and use GUIKit to build a new monitor widget in Mathematica. This provides visual feedback during the potentially long computing time.

This is all the code needed to build a new widget in Mathematica ...

Load the package.

Needs["GUIKit`"]

Set up the computation monitor.

monitorwidget = Widget["Frame", {
  "title" -> "Solution Monitor",
  WidgetSpace[5],
  { 
    WidgetSpace[5],
      Widget["TextField", {"text" -> "", "editable"->False}, Name->"time"],
     WidgetFill[],
     WidgetSpace[5], Widget["Button", {"text" -> "Stop",
      BindEvent["action",
        Script[
          stopRequested = True;
          ]
       ]}],
       WidgetSpace[5]},

  Widget["ImageLabel", {
    "preferredSize" -> Widget["Dimension", {"width"->288, "height"->288}]
    }, WidgetLayout -> {"Stretching" -> {Maximize, Maximize}}, Name -> "imagePanel"],
    
  Script[
    stopRequested = False;
    ]
    
  }, Name -> "monitorFrame"];

Define the action to be the display for each frame.

Δt  =  60;   
showstep  =  1;   
\(action\  := \ \((
timem\  = \ 
          StringJoin["t = ", \ 
            ToString[Round[t/60]], \ "min  RAM in use = ", 
            ToString[ .1 
                 Round[MemoryInUse[]/
                    10^5]], \ " MB"]; 
Sow[{t, \ 
            timem}]; 
monitor\ @\ 
          SetPropertyValue[{"time", \ "text"}, 
            timem]; 
monitor\ @\ 
          InvokeMethod[{"imagePanel", "setData"}, 
ExportString[\ 
              ListPlot3D[
                Take[h[λ, \ ϕ, \ t], {1, \(-1\), \ 
                    showstep}, {1, \(-1\), \ showstep}], \ 
                Mesh -> False, PlotRange -> {\(- .25\), 1}, 
                DisplayFunction -> \ 
                  Identity], \ "GIF"]]; 
If[\ 
          TrueQ[\ monitor\ @\ Script[stopRequested]], \ 
          monitor\ @\ Script[stopRequested\  = \ False]; 
          Throw[Null, \ "StopIntegration"]];
)\)\)

This starts the monitor widget.

monitor  =  GUIRun[monitorwidget];

Modeling, Simulating, and Analyzing a Tsunami

Simulating the Model – Low Resolution Simulation

For this faster simulation we will choose a grid spacing.

GridSpacing[nx, ny]

Start the monitoring widget. [Go ahead and evaluate]

monitor  =  GUIRun[monitorwidget];

To simulate the model you simply use NDSolve, which solves the partial differential equations for the given initial conditions. [Go ahead and evaluate]

sol  = 
      First@NDSolve[{swe,  
            h[λ,  ϕ,  0]  ==  
              h0[λ,  ϕ],  
            u[λ,  ϕ,  0]  ==  0,  
            v[λ,  ϕ,  0]  ==  0,  
            h[λ,  (Subscript[ϕ, min]),  
                t]  ==   h[λ,  (Subscript[ϕ, max]),  t],  
            h[(Subscript[λ, min]),  ϕ,  
                t]  ==   h[ (Subscript[λ, max]),  ϕ,  t]}, 
          h,  {λ,  (Subscript[λ, min]),   (Subscript[λ, max])},  {ϕ, (Subscript[ϕ, min]),  (Subscript[ϕ, max])},  {t, 0, T},  methodOptions ]

In this case we can see that there is a lot of noise in the solution.We need more data points... and hence more memory.

Close the monitoring widget. [Go ahead and evaluate]

CloseGUIObject[monitor];

The solution of the system of partial differential equations from a specified equation and initial/boundary values is a multistep process that makes extensive use of Mathematica's symbolic and numerical capabilities along with automatic algorithm selection.  Here are some of the key steps in the process:

Modeling, Simulating, and Analyzing a Tsunami

Simulating the Model – 32-bit vs 64-bit

This summarizes some of the memory requirements and what it means in terms of the resolution of detail in the model.

The largest 32-bit simulation on a 2GB machine uses 176 grid points in each direction. This gives the following spatial resolution.

n  =  11*16;
nx  =  n  +  1; ny  = n  +  1;  
    nxy  =  (nx  -  1)*(ny - 1);
GridSpacing[nx, ny]

The largest 64-bit simulation is really only limited by physical RAM. On an 8GB machine we ran a model using 360 grid points in each direction. With the spatial resolution below.

n  = 5*8*9;
nx  =  n  +  1; ny  = n  +  1;  
    nxy  =  (nx  -  1)*(ny - 1);
GridSpacing[nx, ny]

The algorithm used is based on a pseudospectral discretization in space followed by an integration in time.  The pseudospectral discretization provides very good accuracy for relatvely few grid points and is based on computing fast Fourier transforms (FFTs). FFT's are by no means limited to power of two lengths and Mathematica uses an efficient prime factor algorithm for which the FFTs can be computed more efficiently for lengths that have small primes as factors.

Modeling, Simulating, and Analyzing a Tsunami

Analyzing the Model – Tracking the Wave Crest

We would like to get more than a movie from our simulation. In particular we want to analyze some of the properties such as wave height, wave speed, etc. Having the solution in memory makes this very easy.

Set up the height as a function of latitude and time for the mean longitude, and compute its gradient with respect to longitude

hϕ[ϕ_?NumberQ,  t_]  =  
        h[0, ϕ,  t]  /.  sol;
hϕg[ϕ_?NumberQ,  
          t_]  =  {D[
            h[0, ϕ,  t]  /.  
              sol,  ϕ]};
ϕgridp  =  
        Cases[(h  /.  sol)[[3, 2]],  
          x_  /;  (x  >  -15 °)];

The crest of the wave is the maximum of the function.

Crest[t_]  :=  
      Module[{gvals  =  h[0,  ϕgridp,  t]  /.  sol, max},  
        max  =  
          FindMaximum[
            hϕ[ϕ,  t],  {ϕ,  
              Last[Extract[ϕgridp, Position[gvals, Max[gvals]]]] },  
            Gradient -> hϕg[ϕ,  t]]; {t,   ϕ  /.  
            max[[2]],  max[[1]]}]

The following analysis was pre-generated using the high-resolution data.

crests  =  Table[Crest[t], {t,  0, T,  T/120}];

This shows the heights of the crest, in meters, as the wave travels.

ListPlot[crests[[All, {1, 3}]]];

This shows the position of the crests, in radiance latitude, as a function of time.

ListPlot[crests[[All, {1, 2}]]];

Use the Vincenty package to compute distances, in meters, between points on the ocean.

dist  =  
        Map[VincentyDistance[{ Mean[{(Subscript[ϕ, 1]),  (Subscript[ϕ, 2])}], #[[1]]}/
                Degree, { Mean[{(Subscript[ϕ, 1]),  (Subscript[ϕ, 2])}], #[[2]]}/
                Degree] &, 
Partition[crests[[All, 2]], 
            2]];

ListPlot[dist[[All, 1]]];

The time between points is one minute, so we get the wave travel speed of 720 km/hour (440 mi/hour).

12000 (60/1000) km/hour

Modeling, Simulating, and Analyzing a Tsunami

Conclusion

• We have modeled, simulated, and analyzed a tsunami wave scenario. This is something that takes months or years to do with traditional programming.

• It is clear that more addressable (and physical) memory gives the ability to:

□ Simulate the model at higher resolution.

□ Actually analyze the results that we get from the model.

• Mathematica offers vast numbers of other areas of technical computing that will immediately benefit from larger addressable (and physical) memory ... it allows us to tackle greater challenges.