Abstract

Most analyses of storm surge and inundation solve equations of continuity and momentum on fixed finite-difference/finite-element meshes. I develop a completely new approach that uses a momentum equation to accelerate bits or balls of water over variable depth topography. The thickness of the water column at any point equals the volume density of balls there. In addition to being more intuitive than traditional methods, the tsunami ball approach has several advantages. (a) By tracking water balls of fixed volume, the continuity equation is satisfied automatically and the advection term in the momentum equation becomes unnecessary. (b) The procedure is meshless in the finite-difference/finite-element sense. (c) Tsunami balls care little if they find themselves in the ocean or inundating land. (d) Tsunami ball calculations of storm surge can be done on a laptop computer. I demonstrate and calibrate the method by simulating storm surge and inundation around New Orleans, Louisiana caused by Hurricane Katrina in 2005 and by comparing model predictions with field observations. To illustrate the flexibility of the tsunami ball technique, I run two “What If” hurricane scenarios—Katrina over Savannah, Georgia and Katrina over Cape Cod, Massachusetts.

1. Introduction

Traditionally, storm surge and inundation have been modeled using finite-difference/finite-element approaches on fixed meshes. A few forefront computer codes of this type (e.g., ADCIRC; see http://adcirc.org/) have been under group development for decades and are true marvels of complexity. Arguably, however, the downside of finite-difference/finite-element approaches may be complexity itself. Researchers interested in trying their hand at storm surge face a steep learning curve to operate ADCIRC-like codes. Moreover, most large-scale finite-difference/finite-element programs are geared specifically for multiprocessor supercomputers. The need for supercomputing to run traditional storm surge codes may form an insurmountable roadblock for scientists who wish to contribute to, or better understand the subject. Might there be a simpler, more intuitive alternative to storm surge than finite difference/finite elements? Could such an approach be scaled down to laptop computer scale?

This article addresses these questions by developing a new approach to storm surge and inundation modeling. The concept springs from Ward and Day [1] who modeled wave runup and inundation using “tsunami balls.” Simply, tsunami balls are bits of water accelerated over a 3D surface. The volume density of balls at any point equals the thickness of the water column there. In storm surge applications, the forces that accelerate tsunami balls derive from surface wind drag, air pressure, gravity, corilois, and bottom friction. As framed earlier, the motivation of this work is not to compete with numerical codes like ADCIRC, but to provide an alternative; an alternative aimed at researchers who wish to learn from 100 km scale storm surge simulations run on their laptop computers.

2. Storm Surge Basics

Let ̂𝐱 and ̂𝐲 directions be east and north in the horizontal plane and ̂𝐳 be up. At vector position 𝐱=(𝑥,𝑦), let still water depth to seafloor be (𝐱) measured positive downward and the perturbation of the surface about the still level be 𝜁(𝐱,𝑡) measured positive upward (Figure 1). Most storm surge calculations solve depth-averaged continuity and momentum equations for variation in water column thickness 𝐻(𝐱,𝑡)=(𝐱)+𝜁(𝐱,𝑡) and mean horizontal water velocity 𝐯(𝐱,𝑡) at fixed mesh nodes: In (2), 𝜌 is the acceleration of gravity, m3(1000 kg/𝜌𝑎) and m3 (1.2 kg/𝑓) are the densities of water and air, and 𝑓=2 is the corilois parameter 7.2921×105s1 (sin(𝜃))𝐶𝑎, with θbeing latitude. Also, 𝐶𝑑 and 𝑃(𝐱,𝑡) are surface and bottom drag coefficients. Equations (1) and (2) are solved given prescribed surface air pressure 𝐕𝑎(𝐱,𝑡) and surface wind velocity (𝐱). Be aware that shallow water equations (1) and (2) assume that surface deformations have wavelengths much longer than water depth (𝐱). If for any reason, the forcing functions on the right-hand side develop short-wavelength components, then this assumption breaks down and solutions to (1) and (2) become unstable.

3. Tsunami Ball Approach

Replacing the established finite-difference/finite-element approach, I intend to employ a momentum equation to accelerate 𝐻(𝐱,𝑡)=𝑁𝑗=1𝑉𝑗𝐴(𝐱,𝐱𝑗(𝑡)).(3) tsunami balls around on a tabletop of varying topography. The tsunami balls are treated as point masses and the thickness of the water column at any location is calculated from the volume density of tsunami balls near that location: In (3), 𝑗 is the fixed water volume of the 𝐱𝑗(𝑡)th ball, 𝑡 is the ball's location at time 𝐴(𝐱,𝐱𝑗(𝑡)), and 𝐴(𝐱,𝐱𝑗(𝑡)) is an averaging area. For instance, 𝑅 might span a circle of radius 𝐱 about 𝐴(𝐱,𝐱𝑗(𝑡))=𝜋𝑅2|𝐱𝐱𝑗(𝑡)|<𝑅=|𝐱𝐱𝑗(𝑡)|>𝑅.(4), then

The primary advantages of a tsunami ball approach to storm surge are four.

(1)Because I track water balls of fixed volume, the continuity equation (1) is satisfied automatically. In locations where the volume density (3) of balls grows, the water column thickness increases. In locations where the volume density drops, the thickness of the water column falls.(2) The tsunami ball procedure is meshless. Meshless applications offer a huge simplification over finite-difference/finite-element methods in that one can download a coastal DEM and start the calculation immediately without having to worry about mesh density, node locations, and so forth.(3)Tsunami balls care little if they find themselves in the ocean (𝐱)<0 or if they have blown onto land 𝐴(𝐱,𝐱𝑗(𝑡)) (Figure 1). Such indifference is made-to-order for inundation applications because it obviates the need for special “dry cell” and “wet cell” behaviors.(4) When programmed efficiently, large-scale (several 100 km) tsunami ball simulations of storm surge can be implemented on a laptop computer.

The primary disadvantage to the tsunami ball approach is granularity. Questions like “How many balls are needed?;” or, “How should averaging area 𝜁(𝐱,𝑡) be selected?” are always close at hand. Smoothness is a particular concern in evaluating the surface gradient term 𝐱 in (2).

3.1. Advection

As was mentioned, basic storm surge calculations solve (1) and (2) at fixed mesh nodes 𝐱. A portion of the change in water velocity at 𝐱 as described by (2) represents real accelerations of the water in the vicinity. Another part of the velocity change accounts for the fact that the water particles at 𝑡 at time 𝑡+𝑑𝑡 are not the same ones that will be there at 𝐱. If the new particles have different velocities than the old ones, 𝐯(𝐱,𝑡)𝐯(𝐱,𝑡) experiences an apparent (or advected) acceleration. Here, because I intend to track specific bits of water, the advection term 𝐱𝐱𝑐𝐑=𝑅 in (2) is not needed. In the tsunami ball approach, quantities like (3), evaluated at fixed locations automatically, account for the fact that different tsunami balls contribute to the count at different times.

4. Application to Hurricanes

This article simulates storm surge due to hurricanes. The primary goal is to introduce and demonstrate the utility of the tsunami ball approach and not to reproduce the details of any specific situation. Accordingly I use “off the shelf” prescriptions of hurricane driving forces and do not dwell deeply on any choice of parameter.

4.1. Surface Pressure and Wind Velocity

If 𝐱 is the vector from the hurricane center to position 𝑃(𝑅)=𝑃𝑐(𝑃𝑛𝑃𝑐)𝑒(𝑅𝑚/𝑅)𝐵𝐕,(5)𝑎𝑅(𝑅)=𝜌𝑎𝜕𝑃(𝑅)+𝑅𝜕𝑅2𝑓241/2𝑅𝑓2(̂𝐳×𝐑)(6), the Holland wind model [2] gives generic surface pressure and zonal (tangential) wind velocity as (Figure 2) with the 𝑅=𝑅𝑚 in (5) taken such that (6) gives a peak wind speed 𝐵2.7183𝑉2𝑚𝜌𝑎𝑃(𝑅𝑚)𝑃𝑐.(7) at radius 𝐱𝑐, To evaluate (5)–(7) one needs to specify as a function of time:

(a)the storm's central position 𝑉𝑚;(b)the storm's central pressure 𝑅𝑚;(c) the storm's peak wind speed 𝑃𝑛;(d) the storm's radius of maximum velocity 𝑅;(e)the background air pressure 𝐕𝑎(𝑅)=|𝐕𝑎̂(𝑅)|[cos𝛽(𝑅)(𝐳×𝐑)sin𝛽(𝑅)𝐑].(8).

I augment the purely zonal wind model (6) in two ways.

(a) As a consequence of friction, the direction of the encircling wind (6) is rotated inward by positive angle β(𝛽(𝑅)=25;𝑅>1.2𝑅𝑚,𝛽(𝑅)=10+75𝑅𝑅𝑚1;𝑅𝑚<𝑅<1.2𝑅𝑚,𝛽(𝑅)=10𝑅𝑅𝑚;𝑅<𝑅𝑚.(9)) Following Martino et al. [3], I take this inward angle to be

(b) Also following Martino et al. [3], I include in the forcing wind an additional component in the direction of hurricane storm velocity 𝐶𝑎 as The inverse barometer effect (see Section 4.4), together with the inward (8) and along track wind components (10), builds and pulls along a pile of water near the hurricane's low-pressure eye. This water eye or “eye-pile” is larger in diameter than the pressure eye and increases in height as it moves into coastal shallows and washes onto land when the storm crosses a shoreline. Water in the eye-pile, together with the water-driven ashore by the zonal winds, comprises the two components of hurricane storm surge.

4.2. Surface Wind Drag 𝐶𝑎(𝑉𝑎)=(0.75+0.067𝑉𝑎)×103,(11)

Many formulations exist for fixing the wind drag coefficient 𝑉𝑎 as a function of wind speed. I employ a standard linear relation from Garratt [4]: where 𝐶𝑑 is the wind speed in m/s, 10 m above the sea surface. In this article, 𝐶𝑑 will be the absolute value of (10).

4.3. Bottom Resistance 𝐶𝑑

For tsunami balls in the ocean, I select a bottom drag coefficient 𝐶𝑎 = 0.001. When tsunami balls blow onto land, presumably they encounter more resistance to flow from vegetation and structures. I increase 𝐶𝑑 for those balls to 0.005.

4.4. “Fly Apart” Surge—Limits on 𝑅 and 𝑣

When viewed from a particle perspective, certain storm surge relationships emerge that might not be apparent in a fixed node perspective. Consider a stationary hurricane with purely zonal winds in which the ocean has reached a steady state with bits of water smoothly orbiting the hurricane pressure eye at distance 𝐴cent=𝑣2𝑅.(12) with constant speed 𝜌𝑎𝐶𝑎𝑉2𝑎𝜌=𝐶𝑑𝑣2or𝑣2=𝜌𝑎𝐶𝑎𝑉2𝑎𝜌𝐶𝑑.(13). From a particle viewpoint, to orbit stablely, water bits must continually experience an inward centripetal acceleration of If inward accelerations fall below (12), the water particles composing the eye-pile will “fly apart.” Steady-state water velocity will be attained when the surface and seafloor drag accelerations in (2) balance, that is, The steady-state inward accelerations are 𝑔𝜁(𝑅), so critical condition (12) becomes The hurricane pressure eye accompanies a water surface high and an air pressure low, so 𝑔𝜁(𝑅)=𝜌1𝑃(𝑅)or(15a)𝜁(𝑅)=(𝜌𝑔)1[𝑃𝑛-𝑃(𝑅)],(15b) is directed outward and [𝑃𝑛𝑃𝑐] is directed inward. Elementary notions ignore the orbital motions of the water and simply strike a balance in these forces to generate the inverted barometer concept: where the height of the eye-pile is proportional to the pressure low. In deep water, for a pressure difference 𝜕𝑃(𝑅)=𝜌𝜕𝑅𝑎𝐶𝑎𝑉2𝑎𝐶𝑑𝑅.(16) = 100 mbar (𝑉𝑎𝑅(𝑅)=𝜌𝑎𝜕𝑃(𝑅)𝜕𝑅1/2whence𝜕𝑃(𝑅)=𝜌𝜕𝑅𝑎𝑉2𝑎(𝑅)𝑅.(17) Pa), the inverted barometer effect (15b) draws up an eye-pile of about 1 m height—a rather minor contributor to surge. Shoaling of the tracking water as the storm runs ashore, however, amplifies the eye-pile height several times.

Returning to stability condition (14) and dropping the 𝐶𝑎𝐶𝑑<1,(18) term for the moment From (6), air speeds are approximately Under these assumptions, stability condition (16) reduces to a remarkably simple form: for a stable state to exist 𝐶𝑑. If 𝐶𝑑, then centrifugal forces will eventually cause the eye-pile to fly apart. If 𝐻(𝐱,𝑡) is 0.001 as I assume, Garratt's [4] formula (11) for 𝐱 violates condition (18) in winds of just 10 m/s.

Of the several assumptions leading to (18), the strongest is that of steady state. Numerical tests using Garratt's [4] formula and 𝐱 = 0.001 (0.004) indicate that in 100 m of water, it takes 11.7 hours (5.8 hours) of sustained 50 m/s winds to accelerate the fluid to 90% of its limited speed of 3.5 m/s (1.7 m/s). Time to steady state increases in proportion to the water depth, so it is unlikely that a travelling storm will remain over a given location long enough to accelerate the water below to steady state speed (13). Still, the concept of water orbiting a hurricane eye with pressure and topographic forces summing to balance centripetal forces is at least as useful as the inverted barometer concept that ignores zonal water velocities all together.

5. Computational Elements

Like every numerical method, tsunami balls utilize certain techniques to stabilize and to speed the calculation. Moreover, because a primary intent is to keep storm surge calculations to a laptop computer scale, additional shortcuts become necessary. Below I list some “tricks of the trade.”

5.1. Grid Not Mesh

I have touted an advantage of tsunami balls as being meshless in the finite-difference/finite-element sense. Still, it is convenient to interpolate certain quantities to a grid. A grid differs from a mesh in that the former is a simple rectangular frame and its selection does not influence the calculation greatly. Interpolation to a grid involves the weighted density of tsunami balls in the area. Equation (3) for instance, evaluates the water column thickness 1𝐯(𝐱,𝑡)=𝐻(𝐱,𝑡)𝑁𝑗=1𝑉𝑗𝐴(𝐱,𝐱𝑗𝐯(𝑡))𝑗(𝑡),(19) at fixed grid point 𝐱𝑗(𝑡). Interpolating flow velocity to grid point 𝐯𝑗(𝑡) might involve where 𝑡 and 𝑁 are the 𝑀th ball's location and velocity at time 𝐻(𝐱,𝑡). Interpolating quantities associated with 𝜁(𝐱,𝑡) balls to 𝜕𝐯𝑗(𝑡)𝜕𝑡=𝑔𝜁(𝐱,𝑡)𝜌1𝑃(𝐱𝑗+𝜌(𝑡),𝑡)𝑎𝐶𝑎|𝐕𝑎(𝐱𝑗(𝑡),𝑡)|𝐕𝑎(𝐱𝑗(𝑡),𝑡)+𝐶𝜌𝐻(𝐱,𝑡)𝑑|𝐯𝑗(𝑡)|𝐯𝑗(𝑡)𝐻(𝐱,𝑡)+𝑓(𝐯𝑗̂(𝑡)×𝐳)(20) grid locations by direct summations like (19) is computationally costly however. In practice, I smooth and interpolate quantities simultaneously (see Section 5.5 in what follows).

Grid-interpolated quantities 𝑃(𝐱𝑗(𝑡),𝑡) and 𝐕𝑎(𝐱𝑗(𝑡),𝑡) appear in the variant of momentum equation used to accelerate the tsunami balls: 𝐻(𝐱,𝑡) and 𝜁(𝐱,𝑡) are evaluated explicitly at 𝐱𝑗(𝑡) from (5) and (10). Quantities 𝑔 and 𝐯𝑗(𝑡) are smoothed interpolated values read from the grid point nearest to 𝜁(𝐱,𝑡).

5.2. One Way Gravity

Equations (1) and (2) arise from shallow water wave theory. Not surprisingly, they generate waves in addition to flows. By and large, oscillatory waves contribute less to storm surge than do longer period flows. Also, compared to flows, waves vary much more rapidly in time and space. Resolving wave actions requires far finer steps in time and space than does capturing flow actions. The presence of waves also contributes to numerical instabilities if the long wave assumption inherent in (1) and (2) becomes violated. In view of this, I invoke “one way” gravity to damp out waves as quickly as possible. One way gravity selects the value of 𝜁(𝐱,𝑡)𝐯𝑗(𝑡)<0;then𝑔=9.8m/s,(21a)𝜁(𝐱,𝑡)𝐯𝑗(𝑡)>0,|𝜁(𝐱,𝑡)|<105;then𝑔=9.8m/s,(21b)𝜁(𝐱,𝑡)𝐯𝑗(𝑡)>0,|𝜁(𝐱,𝑡)|>105;then𝑔=0m/s.(21c) in momentum equation (20) at each location and time based on the current velocity 𝑔𝜁(𝐱,𝑡) of the tsunami ball being accelerated and the gradient ofthe surface slope 𝐯𝑗(𝑡) at the nearest grid location. Specifically, if In words, if acceleration 𝜁(𝐱,𝑡)=𝑃(𝐱,𝑡)𝑔𝜌[104Pa/100km]9.8×103105.(22) opposes ball velocity 𝜌𝜁(𝐱,𝑡)=𝑎𝐶𝑎|𝐕𝑎(𝐱,𝑡)|𝐕𝑎(𝐱,𝑡)𝑔𝜌𝐻(𝐱,𝑡).(23) (21a) or if the surface slope 𝜁(𝐱,𝑡)=(3.04×104)m𝐻(𝐱,𝑡)or(1.25×103)m𝐻(𝐱,𝑡).(24) is small (21b), gravity acts in full force. If acceleration |𝜁(𝐱,𝑡)|105 reinforces velocity 104 and the surface slope is not small (21c), gravity does not act at all. How does one way gravity suppress waves? Think of a playground swing—by fully opposing the upswing (21a) and by not accelerating most of the downswing (21c), one way gravity effectively damps the oscillation. One way gravity curbs wave actions while hardly affecting long-maintained, low-slope flows and surges associated with wind drag and air pressure.

What surge slopes do we expect from wind drag and air pressure? For pressure alone, (15a) and (15b) give From (2), a static force balance predicts surge slopes due to winds alone of Using (11) under winds of 30 m/s or 50 m/s, storm surge slopes would be In 50–100 m of water, both (22) and (24) predict pressure or wind-driven storm surge slopes to be about 𝑁. In shallow water (10 m), wind-driven storm surge slopes might reach 𝑀 or 1 m in 10 km. Water slopes larger than Δ𝑇 likely associate with wave actions and should be extinguished. This line of reasoning motivated the parameter choice in (21a), (21b), and (21c).

5.3. Inner and Outer Loops

One wishes to evaluate (20) and update tsunami ball positions at time intervals such that 𝜁(𝐱,𝑡x) times the ball velocity is less than spacing of base topography, otherwise, the balls might skip over obstructions. On the other hand, updating ocean surface 𝑡x on 𝑡x grid points from 𝑁 balls (see Section 5.5) takes considerable computational effort, so I want to do this less frequently. Accordingly, this calculation incorporates inner and outer time loops. In the inner loop, tsunami ball positions are updated using (20) with small 𝑀, but with the water surface 𝑁 fixed at recent time 𝑀. Only in the outer loop is the surface shape re-evaluated to a new 𝑁×𝑀. The calculation passes 20 or 30 inner loops for each pass of the outer loop.

5.4. Balls at Domain Boundaries

During the calculation, accelerations (20) may push tsunami balls off of the computational grid. I replace out-of-bounds balls at the closest in-bounds location and then zero the ball's velocity. Near domain walls there will be artifacts, so I try to keep the walls as far as possible from the region of interest.

5.5. Rectangular Smoothing

Interpolating quantities to the grid must be done efficiently to hold this calculation to laptop scale. Summing quantities over 𝑁 balls at 𝑀 grid locations directly by (3) or (19) will not be tenable when 𝐱𝑗(𝑡) and 𝑞 each count many 100,000. Instead, I depend on a two-step smoothing/interpolation that reduces the number of numerical operations from being proportional to 𝐺𝑞=Σ𝑞Norm;Σ𝑞=𝑞+𝑊𝑟=𝑞𝑤𝐺𝑟,(25) to being proportional to only 𝐺𝑟. First, cycle through all 𝑊 ball locations 𝑞+1 and assign that ball's value (volume, velocity, etc.) to the nearest grid point. Certain grid points may have zero value while others may have the summed values of several balls. Second, apply a rectangular smoothing to all grid rows and then to all grid columns. The 𝐺𝑞+1=Σ𝑞+1Norm;Σ𝑞+1=Σ𝑞+𝐺𝑞+𝑤+1𝐺𝑞𝑤.(26)th smoothed value in a grid row might be where 𝑀 are raw values, 𝑊 is the half width of smoothing and Norm is some normalization. The 𝑊th smoothed value is You can see that each value (26) derives from the previous value (25) by only one addition, one subtraction, and one division. The number of numerical operations in rectangular smoothing is proportional to the number of grid points 𝑊 only, and independent of both the number of balls 𝐱𝑐 and the smoothing half width 𝑃𝑐. In most cases, I apply rectangular smoothing of half width 𝑉𝑚 twice to get triangle smoothing of half width 2𝑅𝑚. Proper selection of half width 𝑃𝑛 takes trial and error. I desire a smooth solution, but not so smooth as to lose resolution or to drag the water surface artificially far onto land. Because all of the tsunami balls retain a fixed water volume and no balls leave the grid, ocean volume and mean surface height ought to be conserved in smoothing. I select the normalization in (25) and (26) to make this so.

6. Hurricane Katrina at New Orleans

6.1. Hurricane Katrina Parameters

I have previously listed the hurricane parameters that need to be specified:

(a)storm's central position 𝑃𝑛;(b) storm's central pressure 𝑉𝑚(mph);(c) storm's peak wind speed 𝑃𝑐(mbar);(d)storm's radius of maximum velocity 𝑊; (e)the background air pressure Δ𝑇=10s; The Weather Underground [5] website provided Hurricane Katrina storm track information (a)–(c) (see Table 1). I fixed the radius of maximum velocity at 25 km and background air pressure at 104 = 1008 mbar.

6.2. Hurricane Katrina Simulation-Setup

The Katrina at New Orleans simulations ran under a 2300 × 1600 topographic grid 170 m square with minimum and maximum water depths restricted to 3 m and 400 m, respectively. Spaced about 570 m apart, 250 000 tsunami balls were distributed over all wet locations. The water volume assigned to each ball equaled the ball spacing squared multiplied by the water depth at the initial location. (Actually, tsunami balls might be better described as tsunami columns.) To keep the water surface from spreading far onto land as a consequence of smoothing, the smoothing half width dimension , varied from 50 km at locations far from land, to about 2 km at locations close to, or onshore.

The simulation began at 22 GMT 8/28/05 when the Katrina was well south of the grid and ended at 08 GMT 8/30/05, after the storm moved far north of the grid. I updated tsunami ball location and velocity every in the inner loop and recalculated the ocean surface every 5 minutes in the outer loop. A movie frame was constructed at 15-minute intervals. Primary outputs of the simulation include grid-interpolated flow velocity, current flow depth and peak flow depth versus space and time both off and onshore.

6.3. Hurricane Katrina Simulation-Results

Figure 3 shows five frames from the Katrina simulation. (All of the simulations presented in this article are accompanied by Quicktime movie animations. See Table 2. Please view the Quicktime movie version of Figure 3 at http://es.ucsc.edu/~ward/k-at-no.mov) Panel 3(a) sets the stage with geographic information of the New Orleans region. In the Mississippi Delta, levees along the Mississippi River comprise the local high ground. We will see that these linear levees act as dams for surge incoming from either side. Panel 3(a) also presents the surge color scale used throughout this article—colors red, orange, and yellow being positive, and blue and violet being negative.

In Panel 3(b) (06 GMT 8/29), 8 hours after the calculation began, Katrina moves onto the computational grid at the south. The hurricane eye-pile is 2 m high and 60 km across. In all of the simulations, the wind-driven component of surge is most noticeable where winds blow nearly perpendicular to the coastline. Offshore winds of ~35 m/s have begun to draw down the sea surface west of the Mississippi River Delta. On the north-south trending coasts on the east side of the Delta, directly onshore winds of 30 m/s have begun to send surge overland.

In Panel 3(c) (09 GMT 8/29), Katrina has progressed to within about 30 km of first landfall. Water in the eye-pile-driven under the storm, shoals to nearly 5 m. North of the storm, 10 km of overland flooding reaches the Mississippi River and begins to pond against its eastern levee. Northwest of the storm, offshore winds increasing to 50 m/s draw down the sea level by 4 m.

In Panel (d) (12 GMT 8/29), Katrina has moved 20 km inland and crossed to the east side of the river. The eye-pile component of surge now pushing onto land cannot keep pace with the storm and breaks off from the pressure eye itself. Driven by 40 m/s southwest winds, water from the eye-pile runs over 20 km of flat land and dams to 8 m against the River's western levee. Further north, flood water penetrates north and east of the New Orleans town site.

In Panel (e) (15 GMT 8/29), Katrina has crossed the Mississippi State Coast and lies 20 km inland. Near the pressure eye, a quick wind shift after the storm passes causes a redirection of surge. Westerly winds now blow flood waters on the east side of the river back toward the sea from the Delta lands. About this time, surges at New Orleans town site peak at 3-4 m. Squarely onshore southwest and south winds of 50 m/s now force the highest flows of the entire storm event against the Mississippi State Coast. The model registers 11 m of surge 20–40 km east ofKatrina's track.

In Panel (f) (18 GMT 8/29), Katrina has crossed the upper edge of the topographic grid and continues to travel north and weaken. Maintained westerly winds >20 m/s cause the peak of the surge to migrate eastward along the Mississippi State Coast toward Alabama for several hours. Westerly winds persist in ponding water on the west side of the Mississippi River in the southern Delta. Several meters of surge remain there until the end of the simulation (08 GMT 8/30).

Figure 4 contours peak storm surge from 1 to 7 m. (Please view the Quicktime movie of this Figure at http://es.ucsc.edu/~ward/k-at-no-peak.mov) The “dots-on-a-string” map Katrina's position at 1-hour intervals. In the south, note the residual trail from the eye-pile. The eye-pile component of surge increases in height from 2 m well offshore to over 5 m as its tracking water shoals near the Delta. Highest predicted surges for the entire storm (>9 m) occur along the lower Mississippi River and the Mississippi State Coast. At New Orleans, this simulation predicts 3-4 m of flooding. (At the ~200 m scale of the base topography, I do not expect to reproduce a block by block flooding history in New Orleans. Still, if given higher-resolution topography, I do not see why tsunami balls could not be used for that purpose.) Wind-driven surge, banked against the mainland coast, extends 60 km in the shallow waters offshore and as far east as the Florida panhandle. Surge banked up against the mainland covers the numerous islands of Southern Mississippi and Alabama. The slope of the banked up wind blown surge is consistent with that predicted by (24) in shallow water.

Evident in Figure 4 is the strong asymmetry in peak surge relative to the storm track. Counter clockwise rotation of the zonal winds expose coasts on the right-hand side of the storm to more extensive surge than equally distant locations on the left hand side. The zone of maximum surge extends for 100 km to the right side of the track. From Figure 2, 100 km corresponds to the radius of major storm pressure gradient and the distance at which zonal winds fall to 1/2 their peak strength.

Figure 5 contains three large-scale frames from the Katrina simulation. (Please view the Quicktime movie of Figure 5 at http://es.ucsc.edu/~ward/k-at-no-close.mov) At large scale, it is easier to compare wind and water flow directions. As discussed in Section 4.4, the time required for flows to respond to changes in wind direction increases in proportion to water depth. In shallow water, wind and flow directions track fairly closely, and peak flow speeds reach about 2 m/s. At that speed, overland flows take 2-3 hours to attain maximum inundation distance of 20 km. Offshore, this simulation predicts instantaneous surge to be rather lumpy and bumpy over a several km scale (see Figure 5(c) especially). It is difficult to say for certain if this irregularity attributes to a physical process (e.g., vestiges of wave actions) or is an artifact of granularity and smoothing.

6.4. Comparison with Storm Surge Observations

How well does this tsunami ball simulation using off-the-shelf hurricane parameters predict Katrina's actual surge? The best comparison would be made with time histories of sea level measured at many offshore locations. Multiple measurements make it possible to assess variability of surge in time and space (like the lumps and bumps mentioned earlier) and to average noisy information. Also, offshore measurements suffer less from the vulgarities of overland fluid flow than do onshore measurements. A few permanently moored buoys did survive the storm, but these locate at some distance from the action. I do not know of any satellite radar measurements of sea height during the storm, but possibly some exist. Lacking open water surge data, I fall back to information from post-Katrina surveys that measured peak storm surge height onshore. I find a survey by Fritz et al. [6, 7] to be most complete.

Like tsunami runup data, peak storm surge data has several weaknesses. Firstly, peak surge, being a point measure in space and time, can vary dramatically over short distances. Secondly, peak surge, being an extreme measure not a mean measure, incorporates considerable happenstance. Imagine two waves momentarily interfering constructively here, but not over there. Thirdly, peak surge is the sum of peak flow depth plus topographic elevation. If a point at 5 m elevation experiences a flow depth of 3 m, then surge is listed as 8 m. Discrepancies in observed versus model surge heights may be due to differences in flow depth at that point, differences in model versus observed elevation at that point, or both. (Offshore measurements avoid this added difficulty because in water, peak flow depth equals surge height.) At a given latitude and longitude, elevations extracted from the base topography used in these calculations can differ by several meters from those reported by Fritz et al. [6, 7].

In appreciation of the problems inherent in onshore measurements of peak surge, the comparison proceeds as follows.

(1) Given latitude, longitude, and elevation of an observation location, search the base topography file within a 500 m radius and select a reference location that has an elevation as close as possible to the stated one. Elevations at the reference locations mostly matched observed elevations within 50 cm, but there were cases where no high ground could be found and reference location fell as much as 2 m lower than the observation location.(2)Peak surge height (elevation plus peak flow depth) at a comparison location anywhere within 500 m of the reference location was extracted and compared to field data. Selecting the largest surge anywhere near reference location helps account for happenstance.(3) Lastly, recall that the smoothing process outlined in Section 4.4 ignores topographic elevation. Near the coast, even a 2 km smoothing dimension can artificially lap up water high onto otherwise dry land. I restricted the comparison location to be no more than 2 m higher than the reference location.

The top numbers in Figure 6 lists observed peak surge heights in decimeters at 40 locations selected from the Fritz et al. [6, 7] survey. The highest observed surges (>6 m) are found within a 100 km wide band east of Katrina's track. Like the predicted pattern in Figure 4, surge waters backed against the Mississippi State Coast and extended at least 20 km offshore to engulf the barrier islands in several meters of water. The lower numbers in Figure 6 list computed peak surge heights at the comparison locations, selected following the procedure above. Given all of the weaknesses and caveats regarding peak surge, a first look suggests that the trends of the two datasets, if not the values themselves, agree fairly well.

To quantify “eye ball” agreements of Figure 6, Figure 7 plots observed and modeled peak surge heights, peak flow depths, and elevations for the 40 selected data of Figure 6. Observed and predicted peak surge values (Figure 7(a)) show good correlation—85% of the peak surge values fall within a factor of two of the observed values (red-dashed lines). with 40% of predictions being too large and 60% being too small. Considerable scatter exists, but net bias is low: mean of observations = 5.01 m; mean of predictions = 4.96 m. It is true to say there are noticeable outliers. Some outliers seem to be isolated, others show systematic differences. For instance, the model predicts higher surges for the three lower Delta sites (number 28, number 32, number 15) than observed. Southernmost location number 15 is one of the worst outliers, discrepant by over 6 m. These three sites lie close to the storm track where surge is dominated by the eye-pile, not so much the zonal wind. The large misfit there might suggest that the modeled eye-pile was exaggerated. Possibly too, the selected value for onland bottom friction ( = 0.005) was low. Higher friction might keep more of the surge from reaching these far inland sites.

Investigating the origin for systematic differences in observed versus computed Katrina surge is a topic for follow-on research. Seeing that the hurricane pressure and wind fields were generic, the water surface drag formulation came off-the-shelf, and no tidal corrections were included, there are plenty of avenues for refinement. In fact, the only adjustments that I made were in the bottom friction coefficient and smoothing half-width. Although better fits to data could be made, the thrust of this article is to demonstrate the approach and not necessarily to reproduce any specific situation. From the Katrina at New Orleans example, I am convinced that a tsunami ball approach can capture most of the important physical aspects associated with storm surge and that it can provide a viable alternative to traditional methods.

7. Hurricane Katrina What Ifs

Much of the worth in developing a geophysical simulator lay in transferability. Once a simulator can credibly reproduce an actual event, other similar, but hypothetical, situations can be investigated by rolling over previously established model parameters. What if Katrina had stalled out over the Mississippi Coast for four hours? What if Katrina stuck North Carolina instead of New Orleans? Investigating “What If” scenarios can be especially fruitful on laptop-scale simulators because one can set up the scenario and find answers in just hours.

For the What If cases investigated in this article, I employ the identical storm history and physical parameters that were used for Katrina at New Orleans. I will move the starting point of hurricane track and rotate its azimuth such that hypothetical Katrina storms strike at interesting locations and angles.

7.1. Hurricane Katrina over Savannah

The first target for hypothetical Hurricane Katrina will be the Georgia State Coast near Savannah. Unlike New Orleans, the trend of the Georgia/South Carolina Coast is generally straight with paralleling bathymetric contours. Many low lying drainages intersect this stretch of coast and penetrate far inland. Savannah itself locates at 5–15 m elevation 25 km up the Savannah River. How far will surge be forced up these channels before the storm passes and the winds shift? Would Savannah be flooded in a Katrina event? Although a What If storm track running obliquely onto land might expose more coastline to damage, I run this Katrina ashore nearly at right angle to the coast to make a control case.

Figure 8 shows two frames from the Katrina at Savannah simulation 3 hours apart. (Please view the Quicktime movie of this figure at http://es.ucsc.edu/~ward/k-at-sav.mov) In Panel (a), the storm has been approaching 16 hours and it is now 20 km offshore. Interestingly, compared with Katrina at New Orleans (Figure 3(c), the pattern of surge is still compact, nearly circular, and just 1-2 m high. For most the time of approach, zonal winds have been oblique to the Georgia/South Carolina Coast, so little surge has built up. In Figure 8(b), the storm has run inland about 40 km, and the eye-pile has collided with coast and squashed. The zonal winds now face on the northeast coast and surge has quickly grown to several meters.

Figure 9 contours peak storm surge from 1 to 7 m for Katrina at Savannah. (Please view the Quicktime movie of Figure 9 at http://es.ucsc.edu/~ward/k-at-sav-peak.mov) Figure 9 highlights the right-hand/left-hand asymmetry in surge distribution relative to the storm track. Surges greater than 1 m affect 250 km of coastline northeast of the track yet, just 20 km southwest, positive surge is nonexistent. Like the New Orleans case, peak surges reach 9 meters and the zone of maximum surge extends for 100 km to the right of the track. These two features seem to be dictated by the storm parameters itself and not strongly dependent on the details of coastline. Any surge feature that can be ascribed to storm input, rather than local bathymetry/topography input, simplifies inundation prediction because it carries over case to case.

7.2. Hurricane Katrina at Cape Cod

For the second What If, I run Hurricane Katrina obliquely over Cape Cod, Massachusetts, and into Massachusetts Bay, South and East of Boston. Cape Cod, with coasts trending every direction of the compass, should present more complex flow and surge patterns than the Savannah case. I expect Cape Cod to block the eye-pile and zonal wind surges coming from the open ocean. Will interior coasts in Massachusetts Bay be shielded from surge? Will the eye-pile regenerate quickly enough as the storm passes over the Bay to attack mainland coasts for a second time?

Figure 10 shows two frames from the Katrina at Cape Cod simulation 5 hours apart. (Please view the Quicktime movie of this Figure at http://es.ucsc.edu/~ward/k-at-cod.mov) In Panel (a), the storm has been approaching 9 hours and it has just reached Nantucket Island. In contrast to the Savannah case (Figure 8(a)), the surge pattern is complex and spatially extensive. Because much of the Massachusetts coastline trends at right angles to the zonal winds, 2 m surges have built against coasts as far as 250 km to the Northwest. Surge grows both on the eastern coast of the outer cape and on the western shores inside of Massachusetts Bay. Nantucket Island took a glancing blow, the eye-pile pushed water to 9 m on its eastern side. In Figure 10(b), winds shift rapidly south of the storm. Most locations that suffered onshore zonal winds now experience offshore winds and visa versa. The elbow of Cape Cod gets clobbered, as do interior coasts near Provincetown at the Cape tip. The eye-pile reforms to 5 m height in the shallow waters of Massachusetts Bay. In this example, Cape Cod offers only limited storm surge protection for the interior shores.

Figure 11 contours peak storm surge from 1 to 7 m for Katrina at Cape Cod. (Please view the Quicktime movie of Figure 11 at http://es.ucsc.edu/~ward/k-at-cod-peak.mov) Indeed, an oblique storm approach to a variably angled coast makes a complex pattern of surge. The left-hand/right hand asymmetry relative to storm track, seen so clearly at Savannah, gets disrupted and surge heights very in unexpected ways. For instance, Nantucket receives 9 m on its east side but only about 2 on the west side just 20 km away. The complicated surge patterns in Figure 9 highlight the importance that hurricane path variability has on certain coasts. In these situations, surge would be difficult to forecast beforehand, so underscoring the value in developing rapid simulation schemes of the type that tsunami balls are designed to perform.

8. Conclusions

This article introduces a new approach to storm surge and inundation. It uses a momentum equation to accelerate balls of water over variable depth topography and it computes surge height from the volume density of balls. Compared with traditional finite-element/finite-difference approaches, tsunami balls are more intuitive plus they have the advantages that the continuity equation is satisfied automatically; the procedure is meshless; inundation is simply tsunami balls run onto land; and several hundred km size calculations can be done on a laptop computer.

I demonstrate and validate the tsunami ball method by simulating storm surge and inundation in the New Orleans region from Hurricane Katrina in 2005. Despite the fact that the storm input parameters were off-the-shelf and that the pressure, wind, and drag formulations were generic, predicted, and observed storm surge heights along a wide swath of coast compared well with little mean bias. Storm surge has two components—an element related to the eye-pile and a component related to the zonal winds. The eye-pile component of surge is generally symmetric relative to storm track and confined to about 50 km distance on either side. The zonal wind component can present a much farther flung and asymmetric distribution depending upon the orientation of the coastline relative to the wind.

In addition to Katrina at New Orleans, I consider two hypothetical scenarios—Katrina at Savannah and Katrina at Cape Cod. Although the storm parameters were identical, substantially different surge height outcomes are possible depending on the shape of coastline and the angle of attack of the storm. Such variability is difficult-to-predict beforehand and it highlights the value of rapid surge calculations that tsunami balls provide. I am convinced that tsunami balls capture the important physical aspects associated with storm surge and that the method provides a viable alternative to traditional finite-difference/finite-element approaches.