# Leveraging Google’s Tensor Processing Units for tsunami-risk mitigation
planning in the Pacific Northwest and beyond

Ian Madden  Institute for Computational and Mathematical Engineering, Stanford
University  Simone Marras  Department of Mechanical Engineering & Center for
Applied Mathematics and Statistics, New Jersey Institute of Technology  Jenny
Suckale  Institute for Computational and Mathematical Engineering, Stanford
University  Department of Geophysics, Doerr School of Sustainability, Stanford
University

###### Abstract

Tsunami-risk and flood-risk mitigation planning has particular importance for
communities like those of the Pacific Northwest, where coastlines are
extremely dynamic and a seismically-active subduction zone looms large. The
challenge does not stop here for risk managers: mitigation options have
multiplied since communities have realized the viability and benefits of
nature-based solutions. To identify suitable mitigation options for their
community, risk managers need the ability to rapidly evaluate several
different options through fast and accessible tsunami models, but may lack
high-performance computing infrastructure. The goal of this work is to
leverage the newly developed Google’s Tensor Processing Unit (TPU), a high-
performance hardware accessible via the Google Cloud framework, to enable the
rapid evaluation of different tsunami-risk mitigation strategies available to
all communities. We establish a starting point through a numerical solver of
the nonlinear shallow-water equations that uses a fifth-order Weighted
Essentially Non-Oscillatory method with the Lax-Friedrichs flux splitting, and
a Total Variation Diminishing third-order Runge-Kutta method for time
discretization. We verify numerical solutions through several analytical
solutions and benchmarks, reproduce several findings about one particular
tsunami-risk mitigation strategy, and model tsunami runup at Crescent City,
California whose topography comes from a high-resolution Digital Elevation
Model. The direct measurements of the simulations performance, energy usage,
and ease of execution show that our code could be a first step towards a
community-based, user-friendly virtual laboratory that can be run by a
minimally trained user on the cloud thanks to the ease of use of the Google
Cloud Platform.

##  1 Introduction

The coast of the Pacific Northwest, from Cape Mendocino in California to
Northern Vancouver Island in Canada as depicted in Fig. 1, is located on the
seismically active Cascadia subduction zone [Heaton and Hartzell, 1987,
Petersen et al., 2002]. Along the 1200-km-long Cascadia subduction zone, there
have been no large, shallow subduction earthquakes over the approximately 200
years of modern-data monitoring, but large historic earthquakes have left an
unambiguous imprint on the coastal stratigraphy [Clague, 1997]. Sudden land
level change in tidal marshes and low-lying forests provide testimony of 12
earthquakes over the last 6700 years [Witter et al., 2003], including one
megathrust event that ruptured the entirety of the current Cascadia subduction
zone in 1700 BC [Nelson et al., 1995, Wang et al., 2013]. The event created a
massive tsunami that swept across the entire Pacific Ocean devastating
communities as far away as Japan [Satake et al., 1996, Atwater et al., 2011].
Current seismic-hazard models estimate that the probability of another
magnitude 9+ earthquake happening within the next 50 years is about 14%
[Petersen et al., 2002].

Figure 1: Map of the Cascadia Subduction Zone in the Pacific Northwest of the
United States. Relative location of Crescent City with respect to the
Megathrust Fault Line, with a more detailed picture of the Crescent City
coastline. Esri provided access to the satellite imagery. Crescent City map at
high resolution provided by Maxar. Pacific Northwest Map provided by Earthstar
Geographics.

A magnitude 9+ Cascadia earthquake and tsunami occurring during modern times
would devastate many low-lying communities along the Pacific Northwest. A
recent assessment suggests that deaths and injuries could exceed tens of
thousands and entails economic damages in the order of several billions of
dollars for Washington and Oregon State [see, e.g., Knudson and Bettinardi,
2013, Gordon, 2012], with potentially severe repercussions for the entire
Pacific coast and country as a whole. The tsunami itself would put tens of
thousands at risk of inundation, and threaten the low-lying coastal
communities specifically in the Pacific Northwest with very little warning
time for evacuation Gordon [2012]. But how to confront this risk?
Traditionally, the most common approach to reducing tsunami risk is the
construction of sea walls, but this hardening of the shoreline comes at a
staggering price in terms of the economic construction costs (e.g., 245 miles
of sea walls in Japan cost $12.74 billion) and in terms of long-term negative
impact on coastal ecosystems [Peterson and Lowe, 2009, Dugan and Hubbard,
2010, Bulleri and Chapman, 2010] and shoreline stability [Dean and Dalrymple,
2002, Komar, 1998].

A potentially appealing alternative to sea walls are so-called hybrid
approaches. Hybrid risk mitigation combines nature-based elements and
traditional engineering elements to reduce risk while also providing co-
benefits to communities and ecosystems. An example of a hybrid approach to
tsunami-risk mitigation is a coastal mitigation park: A landscape unit on the
shoreline built specifically to protect communities or critical infrastructure
and provide vertical evacuation space, in the styles of Fig. 2. Communities
across the Pacific Northwest are increasingly considering these nature-based
or hybrid options [Freitag et al., 2011], but many important science questions
regarding protective benefits and optimal design remain open [Lunghino et al.,
2020, Mukherjee et al., 2023]. This gap is particularly concerning given that
existing models show that a careful design is necessary to avoid potential
adverse effects [Lunghino et al., 2020]. The design of current mitigation
parks, such as the one being built in Constitución, Chile, is not yet
underpinned by an in-depth quantification of how different design choices
affect risk-reduction benefits, partly because numerical simulations of
tsunami impacts are computationally expensive.

Figure 2: Map view (top) and side view (bottom) of a proposed tsunami-
mitigation berm as designed by Project Safe Havens. The berm provides vertical
evacuation space for the adjacent community and could also lower the onshore
energy flux that drives the damage created by tsunami impact. We show this
design as one example of a hybrid approach to tsunami-risk mitigation as it
combines an engineered hill and ramp with natural vegetation. Sketches adapted
from Freitag et al. [2011].

The goal of this paper is to leverage Google’s Tensor Process Units (TPUs) for
enabling a fast evaluation of different mitigation park design and ultimately
advancing evidence-based tsunami-risk mitigation planning rooted in
quantitative assessments. TPUs are a new class of hardware accelerators
developed by Google with the primary objective of accelerating machine
learning computations. They are accessible via the Google Cloud Platform
[Jouppi et al., 2017, Google.com, ] and have recently been used for many
different applications in computational physics and numerical analysis [Hu et
al., 2022, Lu et al., 2020b, a, Belletti et al., 2020]. We build upon and
extend an existing implementation [Hu et al., 2022] to simulate the impact of
idealized tsunamis on the coastline. Our implementation of the shallow water
equations includes a non-linear advective term, not considered in [Hu et al.,
2022], within the software framework based on Google’s TensorFlow necessary to
execute code on the TPU. We discretize the new equations using the Weighted
Essentially Non-Ocillatory (WENO) method [Liu et al., 1994] and a third-order
Runge-Kutta in time.

Numerical simulations of tsunamis have contributed to our understanding of and
ability to mitigate wave impacts for many decades now, starting from the early
work by Isozaki and Unoki [1964] for Tokyo Bay, and Ueno [1960] for the
Chilean coast. The ability to capture the rupture mechanism that generates the
initial condition for tsunami propagation then enabled the reproduction of
many historical tsunamis [Aida, 1969, 1974]. Since then, many numerical models
have been developed to simulate tsunami generation [Borrero et al., 2004,
Pelties et al., 2012, López-Venegas et al., 2014, Galvez et al., 2014, Ulrich
et al., 2019], propagation [Titov et al., 2005, LeVeque et al., 2011, Chen et
al., 2014, Allgeyer and Cummins, 2014, Abdolali and Kirby, 2017, Bonev et al.,
2018, Abdolali et al., 2019], and inundation [Lynett, 2007, Park et al., 2013,
Leschka and Oumeraci, 2014, Chen et al., 2014, Marsooli and Wu, 2014, Maza et
al., 2015, Oishi et al., 2015, Prasetyo et al., 2019, Lunghino et al., 2020]
by solving different variations of the shallow water and Navier-Stokes
equations.

The list of existing numerical models is long and was recently reviewed in
Marras and Mandli [2021] and Horrillo et al. [2015]. Some commonly used ones
are FUNWAVE [Kennedy et al., 2000, Shi et al., 2012], pCOULWAVE [Lynett et
al., 2002, Kim and Lynett, 2011], Delft3D [Roelvink and Van Banning, 1995],
GeoCLAW [Berger et al., 2011], NHWAVE [Ma et al., 2012], Tsunami-HySEA [Macías
et al., 2017, 2020b, 2020a], FVCOM [Chen et al., 2003, 2014]. Our work here
relies on well-known numerical techniques to solve idealized tsunami problems.
Its novelty lies in demonstrating the capability and efficiency of TPUs to
solve the non-linear shallow water equations to model tsunamis.

We intentionally use a hardware infrastructure that is relatively easy to use
without specific training in high-performance computing (See the Google Cloud
TPU page at Google.com ), and may become a standard hardware on which physics-
based machine-learning algorithms will be built [Rasp et al., 2018, Mao et
al., 2020, Wessels et al., 2020, Fauzi and Mizutani, 2020, Liu et al., 2021,
Kamiya et al., 2022]. We propose that our implementation is one step towards a
community-based, user-friendly virtual laboratory that can be run by a
minimally trained user on the cloud thanks to the ease of use of the Google
Cloud Platform. The tool, which is freely available on Github at
[tsunamiTPUlab, 2023] under an Apache License, Version 2.0 for collaborative
open source software development, can be modified to include machine learning
capabilities and, eventually, extended to coupled models for earth-quake
generation, inundation, and human interaction.

## Methods

### Numerical approximation

We model tsunami propagation and runup with the 2D non-linear shallow water
equations in the conservative formulation with a source term in a Cartesian
coordinate system. Letting $\mathbf{x}=(x,y)$ denote position, we define
$u(\mathbf{x},t)$ and $v(\mathbf{x},t)$ as the flow velocities in the $x$ and
$y$ directions, respectively. In our implementation, we solve for $h$, $hu$,
and $hv$. We define $h(\mathbf{x},t)$ as the dynamic water height and
$b(\mathbf{x})$ as the imposed bathymetry, meaning that the quantity $h+b$
represents the water surface level. This leads to the following system of
equations, a set very similar to that suggested by Xing and Shu [2005]

$\displaystyle\frac{\partial}{\partial t}h+\frac{\partial}{\partial
x}(hu)+\frac{\partial}{\partial y}(hv)=0$ (1)
$\displaystyle\frac{\partial}{\partial t}(hu)+\frac{\partial}{\partial
x}\left(\frac{(hu)^{2}}{h}+\frac{1}{2}g(h^{2}-b^{2})\right)+\frac{\partial}{\partial
y}(huv)=-g(h+b)\frac{\partial b}{\partial
x}-\frac{gn^{2}\sqrt{(hu)^{2}+(hv)^{2}}}{h^{7/3}}(hu)$ (2)
$\displaystyle\frac{\partial}{\partial t}(hv)+\frac{\partial}{\partial
x}(huv)+\frac{\partial}{\partial
y}\left(\frac{(hv)^{2}}{h}+\frac{1}{2}g(h^{2}-b^{2})\right)=-g(h+b)\frac{\partial
b}{\partial x}-\frac{gn^{2}\sqrt{(hu)^{2}+(hv)^{2}}}{h^{7/3}}(hv)\,,$ (3)

where $g=9.81$ ms-2 is the acceleration of gravity, and $n$ is the Manning
friction coefficient.

For ease of future notation, we let
$\mathbf{u}=\begin{bmatrix}h&hu&hv\end{bmatrix}^{T}$, and we rewrite the above
equations in a vector form, namely:

$\frac{\partial\mathbf{u}}{\partial t}+\frac{\partial\mathbf{F}}{\partial
x}+\frac{\partial\mathbf{G}}{\partial y}=\mathbf{S}$ (4)

where $\mathbf{F}$ and $\mathbf{G}$ are the fluxes in the $x$ and the $y$
directions for the vector $\mathbf{u}$, and $\mathbf{S}$ is a source term
arising from variations in topography and Manning coefficient.

We implement these shallow-water equations using the finite volume method
whereby the half-step flux and height values are determined through a 5th-
order WENO scheme [Liu et al., 1994, Jiang and Shu, 1996]. We approximate
solutions to cell-wise Riemann problems by formulating fluxes using the Lax-
Friedrichs method as in [LeVeque, 2011]. We formulate the bed source term as
suggested by Xing and Shu [2005], and formulate the friction term explicitly
rather than using the implicit process suggested by Xia and Liang [2018]. We
use a 3rd-order, Total Variation Diminishing Runge-Kutta scheme to step the
numerical solution forward in time.

We begin the discretization of the equation in continuous variables $t$, $x$,
and $y$, using respective step sizes $\Delta t$, $\Delta x$, and $\Delta y$,
which indicate the distance between consecutive integral steps in the discrete
indices $n$, $i$, and $j$, respectively. We use the 5th-order WENO scheme in
the $x$ and $y$ directions, where two values of each quantity $h$, $hu$, and
$hv$ are determined at each half-step of $x$ and $y$. These two values
correspond to a positive and negative characteristic, due to the nature of the
footprint that is chosen at a given point. In other words, given the
conservative form with relevant variable $\mathbf{u}$, $\mathbf{u}_{i,j}$
centered on a finite volume cell, we label these outputs of WENO:

$\mathbf{u}_{i+\frac{1}{2},j}^{+},\mathbf{u}_{i+\frac{1}{2},j}^{-}\quad\text{for
WENO in $x$, \quad
or}\quad\mathbf{u}_{i,j+\frac{1}{2}}^{+},\mathbf{u}_{i,j+\frac{1}{2}}^{-}\text{for
WENO in $y$.}$ (5)

From here, we use the Lax-Friedrichs method to approximate flux values that
serve as solutions to the Riemann problem; i.e. we approximate

$\mathbf{F}_{i+\frac{1}{2},j}=\frac{1}{2}\left[\mathbf{F}(\mathbf{u}_{i+\frac{1}{2},j}^{+})+\mathbf{F}(\mathbf{u}_{i+\frac{1}{2},j}^{-})-\alpha_{\mathbf{u}}\left(\mathbf{F}(\mathbf{u}_{i+\frac{1}{2},j}^{+})-\mathbf{F}(\mathbf{u}_{i+\frac{1}{2},j}^{-}\right)\right]$
(6)

where $\alpha_{\mathbf{u}}$ is the associated Lax-Frierichs global maximum
characteristic speed. Now, we discretize Eq. 4 explicitly as:

$\frac{\mathbf{u}^{n+1}_{i,j}-\mathbf{u}^{n}_{i,j}}{\Delta
t}+\frac{\mathbf{F}^{n}_{i+\frac{1}{2},j}-\mathbf{F}^{n}_{i-\frac{1}{2},j}}{\Delta
x}+\frac{\mathbf{G}^{n}_{i,j+\frac{1}{2}}-\mathbf{G}^{n}_{i,j-\frac{1}{2}}}{\Delta
y}=\mathbf{S}(\mathbf{u}^{n}_{i,j})$ (7)

Note that in our case, we also choose to formulate the source term
$\mathbf{S}(\mathbf{u}^{n}_{i,j})$ explicitly and centered at the grid point.
Since we use an entirely explicit formulation, we can rewrite Eq. 7 as a time
stepping operator for $\mathbf{u}^{n+1}$, namely
$\mathbf{u}^{n+1}=\mathbf{T}(\mathbf{u}^{n})$. Because Runge-Kutta uses
multiple stages within each time step, we reassign the output of the
$\mathbf{T}$ operator to be $\mathbf{u}^{(n+1)}=\mathbf{T}(\mathbf{u}^{n})$,
where$(n+1)$ indicates an intermediate full time step forward. This means a
full Runge-Kutta step progresses as follows:

$\mathbf{u}^{(n+2)}=\mathbf{T}(\mathbf{T}(\mathbf{u}^{n}))\quad\longrightarrow\quad\mathbf{u}^{(n+\frac{3}{2})}=\mathbf{T}(0.25\mathbf{u}^{(n+2)}+0.75\mathbf{u}^{n})\quad\longrightarrow\quad\mathbf{u}^{n+1}=\frac{2}{3}\mathbf{u}^{(n+\frac{3}{2})}+\frac{1}{3}\mathbf{u}^{n}$
(8)

The process outlined by Eq. 8 outputs a final $\mathbf{u}^{n+1}$ representing
a full-step forward in simulation time.

### TPU implementation

To leverage the TPU’s several cores, we divide the domain into multiple
subdomains and independently compute the numerical solution to the governing
equations on each core. While a lot of the computation can take place
independently, each subdomain remains dependent on the others via their
boundaries and the Lax-Friedrichs global maximum in characteristic speed. We
determine global maximum characteristic speed by sharing and reducing the Lax-
Friedrich maximum characteristic speed calculated on each core. We transfer
subdomain boundary information with further care by using a halo exchange. The
data transfer behavior and computation structure is summarized in Fig. 3.

Figure 3: Left: Initialization of implementation takes advantage of CPU to
allocate initial conditions and topography. Center: Regular computation period
occurring on each subdomain, run independently on TPU cores with some data
sharing coordinated by CPU. Right: CPU Gather to write results to output
files.

Our implementation is inspired by Hu et al. [2022], who chose halo exchange as
an instrument for the TPU to communicate information across subdomain
boundaries in their formulation of the shallow water equations. In the halo
exchange process, we transfer slices of the domain from one core to the others
immediately adjacent. While Hu et al. [2022]’s methodology only involved the
exchange of a single slice from one core to the other, we transfer several
slices in order to take full advantage of the high accuracy and larger
footprint of the WENO scheme. These halo exchanges are then performed in every
stage of the Runge-Kutta scheme, meaning that they occur multiple times in a
single time step.

The initial conditions and results are communicated from the remote program,
which resides on the CPU, to the TPU workers by means of tpu.replicate which
sends TensorFlow code to each TPU. We refer to Hu et al. [2022] for further
details on the TPU implementation.

##  2 Model verification and validation

We differentiate between model verification and validation in the manner
suggested by Carson . Specifically, we check for model and implementation
error by quantifying the extent to which numerical solutions compare to
correct analytical solutions [Carson, ]: wet dam break (Section 2.1),
oscillations in a parabolic bowl (Section 2.2), and steady state flow down a
slope with friction (Section 2.3). Following this, we validate by checking how
well numerical solutions reflect the real system and apply to the context
[Carson, ]. To do this, we compare against an existing numerical benchmark
from the Inundation Science and Engineering Cooperative [ISEC, 2004] and
results from an investigation of nature-based solutions [Lunghino et al.,
2020] in Section 2.4, and consider the propagation of a computed tsunami over
the observed topography of Crescent City in Section 2.5. To quantify the
accuracy of the solutions, we test our numerical solver against some classical
analytical solutions to the shallow water equations. We assess the model’s
ability to capture key physical processes relevant to inundation, including
steep wave propagation, friction, and topography dependence. We use relative
errors in the $L_{\infty}$ and $L_{2}$ sense as the metric to determine model
accuracy. These are approximated in this paper in the following manner:

$\text{$L_{\infty}$}=\frac{\max_{\Omega}|h_{c}-h_{a}|}{\max_{\Omega}|h_{a}|},\quad\text{$L_{2}$}=\sqrt{\frac{\sum_{\Omega}(h_{c}-h_{a})^{2}\Delta\Omega}{\sum_{\Omega}(h_{a})^{2}\Delta\Omega}}\,,$
(9)

where $h_{c}$ is the computed solution at the discretized cells, $h_{a}$ is
the analytical solution at the corresponding cells, and $\Omega$ denotes the
computational domain.

###  2.1 Wet dam break

Figure 4: On the left, several instances in time of the computed (purple)
water heights to wet dam break compared with the analytical (orange, dashed)
water heights. The rightmost figure plots the $L_{2}$ and $L_{\infty}$
relative norms of the error between the analytical and computed solutions.

The classical one-dimensional Wet Dam Break [Stoker, 1957] provides us an
opportunity to test the ability of our code to capture shock propagation and
advection. In this case, there is no friction ($n=0$) and the topography is
flat ($b(x)=0$). The boundaries are set at a constant height with zero flux.
We impose the following initial condition:

$(hu)=0,(hv)=0,h(x)=\begin{cases}h_{l}&x\leq x_{0}\\\
h_{r}&x>x_{0}\end{cases}\,,$ (10)

where $h_{l}$ and $h_{r}$ are the constant water heights on either side of a
shock front $x_{0}$. We compare our numerical solution for water height
against the dynamic analytical solution from Delestre et al. [2013]:

$\displaystyle h(x,t)$ $\displaystyle=\begin{cases}h_{l}&x\leq x_{1}\\\
\frac{4}{9g}\left(\sqrt{gh_{l}}-\frac{x-x_{0}}{2t}\right)^{2}&x_{1}(t)<x\leq
x_{2}(t)\\\ \frac{c_{m}^{2}}{g}&x_{2}(t)<x\leq x_{3}(t)\\\
h_{r}&x>x_{3}(t)\end{cases}\,,$ (11) $\displaystyle x_{1}(t)$
$\displaystyle=x_{0}-t\sqrt{gh_{l}}\,,$ (12) $\displaystyle x_{2}(t)$
$\displaystyle=x_{0}+t(2\sqrt{gh_{l}}-3c_{m})\,,$ (13) $\displaystyle
x_{3}(t)$
$\displaystyle=x_{0}+t\frac{2c_{m}^{2}(\sqrt{gh_{l}}-c_{m})}{c_{m}^{2}-gh_{r}}\,,\quad\text{and}$
(14) $\displaystyle c_{m}$ $\displaystyle\text{is the solution
to}-8gh_{r}c_{m}^{2}(\sqrt{gh_{l}}-c_{m})^{2}+(c_{m}^{2}-gh_{r})^{2}(c_{m}^{2}+gh_{r})=0\,.$
(15)

A qualitative comparison of the computed and analytical solutions for times
$t=$0, 2.5, and 9 seconds is shown in the left plots of Fig. 4. The relative
error between the analytical and computed solutions in the infinity and
2-norms at a small distance away from the shock front are plotted on the
right. We interpret the converging relative error norms to a low magnitude as
verification of our implementation to sufficiently capture shock propagation
and advection.

###  2.2 Planar parabolic bowl

Figure 5: On the left, several instances in time of the computed (purple)
water heights to the one-dimensional parabolic bowl compared with the
analytical (orange, dashed) water heights. The rightmost figure plots the
$L_{2}$ and $L_{\infty}$ relative norms of the error between the analytical
and computed solutions.

The classical one-dimensional planar parabolic bowl originally suggested by
[Thacker, 1981], is an oscillating solution allowing us to test the source
term for topography without friction ($n=0$). We enforce homogeneous Dirichlet
conditions in both flux and water height, at a resolution of 1 m. Once again,
we take the test directly from Delestre et al. [2013], resulting in the
following description of the base topography:

$b(x)=h_{0}\left(\frac{1}{a^{2}}\left(x-\frac{L}{2}\right)^{2}-1\right)\,,$
(16)

corresponding with the following initial condition:

$(hu)=0,(hv)=0,h(x)=\begin{cases}-h_{0}\left(\left(\frac{2x-L+1}{2a}\right)^{2}-1\right)&\frac{1-2a+L}{2}<x<\frac{1+2a+L}{2}\\\
0&\text{otherwise}\end{cases}\,.$ (17)

This leads to the following dynamic analytical solution for the water height:

$h(x,t)=\begin{cases}-h_{0}\left(\left(\frac{2x-L}{2a}+\frac{1}{2a}\cos\left(\frac{\sqrt{2gh_{0}}t}{a}\right)\right)^{2}-1\right)&x_{1}(t)<x<x_{1}(t)+2a\\\
0&\text{otherwise}\end{cases}\,,$ (18)

where
$x_{1}(t)=\frac{1}{2}\cos\left(\frac{\sqrt{2gh_{0}}t}{a}\right)-a+\frac{L}{2}$.
A qualitative comparison of the parabolic bowl solution at the time instances
$t=$3.5 sand $t=$12.5 s can be seen on the left of Fig. 5. The analytical and
computed solutions appear to correspond to one another well. For a more
quantitative analysis,the relative error-norms of the solutions are depicted
on the right of Fig. 5. We interpret the converging relative error norms to a
low magnitude as verification of our implementation to sufficiently capture
the source term of the shallow water equations induced by topography.

###  2.3 Steady flow down a slope with friction

Figure 6: On the left, several instances in time of the computed (purple)
fluxes given a water level and a slope, compared to the analytical (orange,
dashed) flux. The rightmost figure plots the $L_{2}$ and $L_{\infty}$ relative
norms of the error between the analytical and computed solutions.

We do a short test in order to assess the correctness the discretized friction
source term, focusing on a relatively simple flow down a slope with finite
friction ($n=0.033$) as tested by Xia and Liang [2018]. The steady state flow
down a slope then becomes

$(hu)=\frac{\sqrt{b_{x}}}{n}h^{\frac{5}{3}}$ (19)

where $b_{x}$ is the slope. In this test, we initialize the problem close to
the steady state solution for a wave height of 0.5 and a slope of
$\frac{1}{20}$. This specific example and its convergence toward steady state
is shown in Fig. 6. The left plots show the flux at $t=100~{}$s and
$t=200~{}$s. The results are almost identical, indicating that a steady state
has been reached. On the right plot, the error norm of the steady state flux
takes some time to reach steady state, but reaches a very small level. Because
we approach the appropriate steady state solution and achieve a very small
error norm, our implementation is verified in capturing a manning friction
law.

###  2.4 Validation for tsunami simulations

To assess the ability of the code to capture tsunami propagation, we start
with a popular numerical benchmark from the Inundation and Science Engineering
Cooperative (ISEC) [ISEC, 2004] that represents tsunami runup over an
idealized planar beach that provides solutions for tsunami runup at times
$t=$, 180 s, 195 s, 220 s. We formulate the initial condition for water height
using Lunghino et al. [2020]. The solutions from the benchmark (dashed,
orange) are qualitatively compared with the numerical solution produced by our
code (solid, purple) in Fig. 7. We take the qualitative agreement as
validation of the model’s ability to model the runup of a Carrier N-Wave
[Carrier et al., 2003].

Figure 7: Qualitative comparison of the computed solution with resolution 1 m
compared to the ISEC benchmark at three different time instances.

Since we are interested in leveraging TPUs for tsunami-risk mitigation
planning, we take a look at the ability of our shallow water equation code to
reproduce a few particular results by Lunghino et al. [2020] who investigated
the effects of hills on a tsunami running up on a planar beach. The tsunami is
initialized as Carrier’s N-wave [Carrier et al., 2003]:

$\eta=2(a_{1}\exp\\{-\hat{k}_{1}(x-\hat{x}_{1})^{2}\\}-a_{2}\exp\\{\hat{k}_{2}(x-\hat{x}_{2})^{2}\\}),$
(20)

where $\eta=h+z$, $\hat{x}_{1}=1000+0.5151125\lambda$,
$\hat{x}_{2}=1000+0.2048\lambda$, $\hat{k}_{1}=28.416/\lambda^{2}$,
$\hat{k}_{2}=256/\lambda^{2}$, $a_{1}=A$, and $a_{2}=A/3$. While this is the
analytically correct form, the flow origin in the code is not the shoreline,
so there are some effective shifts $\hat{x}_{1}$ and $\hat{x}_{2}$ that we
need to do. An example of the Carrier wave initial condition and offshore
propagation behavior for $A=15$ and $\lambda=2000$ is shown in Figure 8. We
apply free slip, no-penetration boundary conditions to the four domain
boundaries, which means that that the component of the boundary-normal
component of the velocity vector is zero whereas its tangential component is
unaltered.

Figure 8: Several snapshots in time of the tsunami’s propagation over a
modeled ellipsoidal hill on a slope. From left to right, the formulation of
the initial Carrier N-wave at $t=0$, followed by the propagation of a wave
front toward the hill at $t=50$, collision of the wave front of the hill at
$t=95$, and the formation of a reflected wave at $t=135$.

The shallow water equation model presented in this study is able to reproduce
the wave reflection provided by a hill, consistent with results from Lunghino
et al. [2020]. Because this simulation is possible by the implementation,
other further analysis can be conducted to understand the mitigative benefit
of other nature-based solutions.

###  2.5 Real-world scenario

Past tsunamis impacting the West Coast of the United States have caused more
damage around the harbor of Crescent City in California than elsewhere along
the Pacific Coast [Arcas and Uslu, 2010]. For this reason, we chose an area of
approximately 105 ${\rm km^{2}}$ around Crescent City to demonstrate the
code’s ability to capture the impact of an idealized tsunami event for a real
location at high resolution. To approximate the actual bathymetry and
topography, we use a Digital Elevation Model for this area with uniform grid
spacing of 4 m provided by NOAA [NOAA National Geophysical Data Center, 2010,
Grothe et al., 2011]. For the sake of providing a proof of concept, we
initialize the tsunami with the idealized waveform described above, with
slightly adjusted parameters: $A=10$, $\lambda=2000$,
$\hat{x}_{1}=6000+0.5151125\lambda$, and $\hat{x}_{2}=6000+0.2048\lambda$. The
chosen parameters lead to maximum inundation patterns similar to that seen in
one modeled extreme scenario from Arcas and Uslu [2010]. In Fig. 9, we start
with an absence of any nearshore wave (including at $t=50s$) and then a
development of a tsunami front that is visible to the shoreline by $t=140s$.
That front penetrates the harbor by $t=220s$, and is soon followed by the
inundation of the coastline as well as reflection of wave energy back to the
ocean. We also observe that the mountain range on the upper part of the figure
clearly provides a significant protective benefit to the land beyond it.

Figure 9: Several snapshots of modeled tsunami propagation over terrain and
geological features of Crescent City, CA, where any level of blue indicates
water cover and green depicts a stylized map of the topography above surface
level. From left to right, then top to bottom, we have steady state near shore
at $t=50$; followed by the propagation of a wave front at $t=100$ and $140$;
contact with Crescent City harbor at $t=180$; inundation of the harbor and
some of the coastline at $t=220$ and $260$; and tsunami reflection and
inundation at $t=300$ and $340$.

##  3 Performance analysis on TPU

###  3.1 Number of TPU cores

In communities where users may not have access to high performance computing
facilities, the Cloud TPU Platform provides a unique ability to perform large-
scale computations, and perform them rapidly. To address this potential
benefit, we first measure the wall-clock time of a simulation of a tsunami
reaching Crescent City using a varying number of TPU cores on one TPU device.
As shown in Table 1, the problem size posed by the realistic scenario is
sufficient to see rapid improvements in runtime based on the number of cores.

Number of Cores | 1 | 2 | 4 | 8
---|---|---|---|---
TPU Runtime [s] | 6894 | 3876 | 2000 | 1036
Table 1: Approximate TPU Runtimes (in seconds) with varying numbers of TPU
cores for Crescent City Configuration using time step of $\Delta t=5\cdot
10^{-3}$, Total Array Size of approximately 1802 by 3984 elements (4 meter
resolution), cores all aligned in the y-direction as suggested in [Hu et al.,
2022]. This runtime excludes transfer times between the CPU and TPU.

###  3.2 Geophysical problem resolution

The typical realistic scenario that decision-makers will face will involve
large problem sizes due to both the extent of their spatial domain, but also
the level of resolution necessary to model tsunami propagation and inundation
over complex topography. Therefore, we check convergence and runtime under
varying degrees of resolution for the current realistic scenario as well,
depicted graphically in 10 and in table 2.

Figure 10: Graphical depiction of the TPU runtimes (purple) and relative error norms (blue) under varying resolutions for computing the tsunami propagation over the Crescent City DEM. Resolution [m] | 2 | 4 | 6 | 8 | 10 | 12
---|---|---|---|---|---|---
Runtime [s] | 3805.167 | 1035.719 | 506.0652 | 337.9506 | 241.0804 | 187.3016
Number of Elements | 28713068 | 7179168 | 3192512 | 1794792 | 1149274 | 798128
Efficiency | 0.000133 | 0.000144 | 0.000159 | 0.000188 | 0.00021 | 0.000235
$L_{2}$ Error Norm | * | 0.062475 | 0.092618 | 0.115671 | 0.093901 | 0.132735
$L_{\infty}$ Error Norm | * | 0.212536 | 0.340507 | 0.291389 | 0.265669 | 0.259101
Table 2: Approximate TPU Runtimes (in seconds) with varying resolutions for
Crescent City Configuration using time step of $\Delta t=5\cdot 10^{-3}$,
using the 2 m resolution as a benchmark for correctness. Ran on a single TPU
with 8 cores.

We finally perform the same analysis under varying degrees of resolution using
the benchmark from the Inundation Science and Engineering Cooperative [ISEC,
2004] that we previously validated in Fig. 7. A qualitative comparison of the
tsunami propagation under different resolutions are graphically depicted in
Fig. 11 in the top two and bottom left figure. In the bottom right plot of
Fig. 11, we see the expected fall in runtime based on coarser resolution
(purple), and a rise in relative error (with the exception of the highest
resolution 1 m). The corresponding values are documented in 3.

Figure 11: Graphical depiction of the TPU runtimes (purple) and relative error norms (blue) under varying resolutions for computing the tsunami propagation in the ISEC Benchmark. Resolution [m] | 1 | 2 | 4 | 8
---|---|---|---|---
Runtime [s] | 150.0193 | 81.74382 | 56.74928 | 47.3934
Number of Elements | 1060521 | 277761 | 75756 | 18939
Efficiency | 0.000141 | 0.000294 | 0.000749 | 0.002502
$L_{2}$ Error | 0.000206 | 7.19E-05 | 0.00031 | 0.000282
$L_{\infty}$ Error Norm | 0.000305 | 0.000155 | 0.000398 | 0.000326
Table 3: Approximate TPU Runtimes (in seconds) with varying resolutions for
the ISEC Tsunami Benchmark using time step of $\Delta t=5\cdot 10^{-3}$. Ran
on a single TPU with 8 cores.

###  3.3 Comparison with GeoClaw

Figure 12: TPU solution (top row) at several time instances compared to the
GeoClaw solution (bottom row). The arrival of the tsunami front ($t=100,180$),
the inundation of the harbor ($t=260$), and coastal inundation and reflection
is depicted, and relatively comparable.

We used GeoClaw [Clawpack Development Team, 2020, Mandli et al., 2016, Berger
et al., 2011] ran on a single thread of a CPU (Intel i7-8650 with a base
frequency of 1.9 GHz) to compare numerical solutions to a tsunami propagation
in order to assess performance enhancements provided by the TPU. While our
TPU-based code completes the 400 second simulation in approximately 17 minutes
of wall-clock time, the GeoClaw implementation on the CPU takes approximately
630 minutes. Comparisons of the two numerical solutions can be seen in Fig.
12, where the top row includes several instances in time of the TPU numerical
solution, and the bottom row depicts the GeoClaw numerical solution at the
same instances in time. Although some differences can be seen in the geometry
of inundation by $t=380s$ in the rightmost plots, the solutions do generally
appear similar over time, lending credibility to the validity of the numerical
solution.

###  3.4 Energy utilization

Estimates of energy efficiency of computing operations are becoming
increasingly popular, especially in response to progressing climate change
[Fuhrer et al., 2018, Fourestey et al., 2014]. To address this potential
benefit provided by TPUs, we conduct a heuristic analysis of the energy
savings of running these operations on a TPU rather than a CPU. We base an
order-of-magnitude estimate on the claims of Google in the maximum power
efficiency of the TPU, 2 trillion operations per second (TOPS) per Watt
[Google.com, ]. We use this heuristic for power efficiency, consider the
tsunami propagation problem under a 8 m resolution for Crescent City, and
manually approximate the number of floating point operations that our
implementation performs in each quadrature step at this resolution. We
estimate that the TPU performs approximately 11.9 million floating-point
operations for each simulated time step, and we see that each time-step takes
an average of about 4.2 milliseconds, corresponding to about 2.83 TOPS.
Assuming constant power efficiency regardless of capacity, this translates to
1.41 Watts and an energy usage of 1.18 J for each simulated second given our
current time-step configuration. Combining with the total runtime for a 400
modeled-second simulation of Crescent City, we see a total cost of
approximately $4.76\times 10^{2}$ J = $1.32\times 10^{-1}$ Wh of energy. At a
price of 21 cents/kWh in the U.S. at the time of writing this article, this
simulation has a monetary cost of $2.7\times 10^{-3}$ cents.

We use the aforementioned GeoClaw run as our CPU comparison on energy
utilization in order to get an order of magnitude estimate of the power
savings. The particular processor for this CPU comparison, the Intel i7-8650
with a base frequency of 1.9 GHz, has a Thermal Design Power of 15 W for 8
total threads. Each simulated second took approximately 96.4 seconds of CPU
runtime, which translates to approximately 181 Joules for each simulated
second of a grid with 8 m resolution. Given the order of magnitude estimates,
we note that the 400 modeled-second simulation would see a total cost of
approximately $7.24\times 10^{4}$ J = $2.01\times 10^{1}$ Wh of energy, or a
monetary cost of 0.04 cents, a cost multiplicative factor of 20. We linearly
extrapolate to estimate the number of CPU threads needed to match the runtime
speed of our model using the TPU, and find that approximately 37 CPU threads
would be needed. With this in mind, we find that the cost multiplicative
factor for a CPU simulation of performance equivalent to that of a TPU would
be closer to an order of 700. TPU energy savings for high performance are
clearly substantial, and not ignorable.

##  4 Discussion

Sustainable tsunami-risk mitigation in the Pacific Northwest is a challenging
task. Some challenges come from beneath, because previous large subduction
zone earthquakes at Cascadia led to $0.5-1$ m of co-seismic subsidence, the
sudden sinking of land during an earthquake [Wang et al., 2013]. Strong
shaking can also lead to liquefaction [Atwater, 1992, Takada and Atwater,
2004]. Other challenges come from the ocean, where sea-level rise [Church and
White, 2006, Bindoff et al., 2007] and intensifying winter storms [Graham and
Diaz, 2001] have increased wave heights [Ruggiero et al., 2010, Ruggiero,
2013] and accelerated coastal erosion [Ruggiero, 2008]. A recent USGS report
documented rapid shoreline changes at an average rate of almost 1 m/yr across
9,087 individual transects [Ruggiero et al., 2013], suggesting the possibility
that the shoreline might change significantly during the century-long return-
period of large earthquakes in Cascadia [Witter et al., 2003].

The picture that emerges is that of a highly dynamic coastline – maybe too
dynamic for an entirely static approach. Nature is not only continuing to
shape the coastline, but is also a fundamental component of the region’s
cultural heritage, identity and local economy. So, it is maybe not surprising
that the Pacific Northwest is a thought-leader when it comes to designing
hybrid approaches to sustainable climate adaptation through the Green Shores
program [Dalton et al., 2013] and to vertical tsunami evacuation through
Project Safe Haven [Freitag et al., 2011].

Project Safe Havens is a grass-roots approach to reducing tsunami risk mostly
by providing accessible vertical-evacuation options for communities. Many
proposed designs entail reinforced hillscapes like the one shown in figure 2,
intended to dissipate wave energy and provide vertical evacuation space during
tsunami inundation. To build confidence in such a solution and its mitigation
effects, risk managers must be able to quickly and precisely forecast a
tsunami inundation, preferably via a publicly available, centralized modeling
infrastructure.

This paper is meant to be a first step towards a major community based
infrastructure that will allow local authorities around the world to readily
execute tsunami simulation once a tsunami in their proximity has been
detected. We aim to provide a proof-of-concept rather than a complete
implementation. As such, we used a very similar base framework used by Hu et
al. [2022] of halo exchange in combination with a WENO and Runge-Kutta scheme,
which may not be optimally taking advantage of the TPU’s computing structure
and capabilities. We originally chose these schemes to maintain higher order
accuracy and ease of implementation but, eventually, a convolution-based
implementation of the quadrature of the shallow water equations should be
tested for maximum performance utilization of the TPUs.

Because our code is specifically an implementation of the shallow water
equations for the TPU, it is currently unable to model tsunami initiation, or
any fluid structure interactions that may be desired to accompany analysis of
nature-based solutions. Instead, it requires an initial condition for wave
heights and fluxes, meaning a full tsunami simulation would require coupling
the results of a tsunami initiation model as an input. While our
implementation is relatively limited in scope, the model is certainly able to
provide a starting point for a more complete software package for communities
as they evaluate nature-based options for tsunami mitigation.

Though just a starting point for a remotely-available package, we successfully
replicated some of the results found about coastal mitigation parks posed by
Lunghino et al. [2020], and we were able to model tsunami runup over with the
real bathymetry around Crescent City, in California, included in the code by
means of a DEM file from the USGS digital elevation database; for comparison
purposes, we ran the same test using the popular open-source solver GeoClaw
[Clawpack Development Team, 2020, Berger et al., 2011]. The results of the two
models are in good agreement, as shown in Fig. 12. Our contribution lies in
demonstrating an enhanced ability to run high quality simulations using the
TPU available remotely using Google Cloud Platform. As a result, high quality
tsunami simulations are available to remote communities for rapidly evaluating
mitigation or evacuation options when faced with coastal flooding.

We argue that TPUs are preferable to large, heavily parallel simulations on
CPU or GPUs, because the TPU-based simulations we show here do not require
access to large computing clusters. These are usually made available to
scientists and engineers by supercomputing centers around the world by means
of competitive grants for computing time or by use of the cloud offered by
private companies. However, an expert user knowledge of these systems from a
scientific computing perspective is necessary to design, run, and interpret
model results, and the compute infrastructure itself may not be available to
early warning centers in many parts of the world. In contrast, our code is
available on Github and fully implemented in Python, can be ran through a web
browser, and visualized through a simple notebook file using Google Colab.
While performance can be enhanced with some knowledge about TPU architectures,
community risk managers do not need this knowledge to run high quality tsunami
simulations rapidly for real, physical domains with associated DEMs.

Finally, in the face of rising energy costs in both a monetary and climate
sense, the TPU infrastructure allows the support of more energy-efficient
simulations over those of CPU-based clusters. This means that those coastal
flood risk managers in remote communities immediately are able to support
design decisions with model results in a climate-friendly manner.

Although not our focus here, we note our approach may also contribute to early
tsunami warning. Once triggered, tsunamis move fast; this fact makes it
necessary to model and assess their potential for damage ahead of time once
they have been detected offshore. For a sufficiently fast early warning and
prompt evacuation, the tsunami modeling infrastructure has an important time
constraint [Giles et al., 2021] to be considered, and Faster Than Real Time
(FTRT) simulations are necessary [Behrens et al., 2021, Løvholt et al., 2019].
To make FTRT simulations a reality, tsunami models are being rewritten or
adapted to run on Graphical Processing Units (GPUs) [Løvholt et al., 2019,
Behrens and Dias, 2015, Satria et al., 2012]. A TPU-based implementation as
proposed here might be another meaningful step into that direction.

## Author contributions

Ian Madden: Software, Analysis, Writing. Simone Marras, PI: Conceptualization,
Methodology, Writing, Supervision. Jenny Suckale: Conceptualization,
Methodology, Writing, Supervision.

## Competing interests

The authors declare that they have no conflict of interest.

## Data availability statement

Our work is available as a GitHub release at
https://github.com/smarras79/tsunamiTPUlab/releases/tag/v1.0.0 or on archive
at 10.5281/zenodo.7574655.

## Acknowledgements

This work was supported by the National Science Foundation’s Graduate Research
Fellowships Program (GRFP) awarded to Ian Madden.

## References

  * Abdolali and Kirby [2017] A. Abdolali and J. T. Kirby.  Role of compressibility on tsunami propagation.  _J. Geophys. Res.: Oceans_ , 122:9780–9794, 2017.
  * Abdolali et al. [2019] A. Abdolali, U. Kadri, and J. T. Kirby.  Effect of water compressibility, sea-floor elasticity, and field gravitational potential on tsunami phase speed.  _Scientific Reports_ , 9:1–8, 2019.
  * Aida [1969] I. Aida.  Numercal experiments for the tsunami propagation of the 1964 Niigata tsunami and 1968 Tokachi-Oki tsunami.  _Bull. Earthquake Res. Int. Univ. Tokyo_ , 47:673–700, 1969\.
  * Aida [1974] I. Aida.  Numerical computational of a tsunami based on a fault origin model of an earthquake.  _J. Seismol. Soc. Jpn._ , 27:141–154, 1974.
  * Allgeyer and Cummins [2014] S. Allgeyer and P. R. Cummins.  Numerical tsunami simulation including elastic loading and seawater density stratification.  _Geophys. Res. Letters_ , 41:2368–2375, 2014.
  * Arcas and Uslu [2010] D. Arcas and B. Uslu.  PMEL Tsunami Forecast Series: Vol. 2. A Tsunami Forecast Model for Crescent City, California.  Technical Report 3341, NOAA, 2010.
  * Atwater [1992] B. F. Atwater.  Geologic evidence for earthquakes during the past 2000 years along the Copalis River, southern coastal Washington.  _Journal of Geophysical Research: Solid Earth_ , 97(B2):1901–1919, 1992.
  * Atwater et al. [2011] B. F. Atwater, S. Musumi-Rokkaku, K. Satake, Y. Tsuji, and D. K. Yamaguchi.  _The orphan tsunami of 1700: Japanese clues to a parent earthquake in North America_.  University of Washington Press, 2011.
  * Behrens and Dias [2015] J. Behrens and F. Dias.  New computational methods in tsunami science.  _Phil. Trans. R. Soc. A_ , 373:20140382, 2015.
  * Behrens et al. [2021] J. Behrens, F. Løvholt, F. Jalayer, S. Lorito, M. A. Salgado-Gálvez, M. Sørensen, S. Abadie, I. Aguirre-Ayerbe, I. Aniel-Quiroga, A. Babeyko, et al.  Probabilistic tsunami hazard and risk analysis: a review of research gaps.  _Frontiers in Earth Science_ , 9:628772, 2021.
  * Belletti et al. [2020] F. Belletti, K. King, D. Yang, R. Nelet, Y. Shafi, and J. Shen, Y.-F.and Anderson.  Tensor processing units for financial monte carlo.  In _Proceedings of the 2020 SIAM Conference on Parallel Processing for Scientific Computing. SIAM_ , pages 12–23, 2020.
  * Berger et al. [2011] M. Berger, D. George, R. LeVeque, and K. Mandli.  The geoclaw software for depth-averaged flows with adaptive refinement.  _Adv. Water Res._ , 34(9):1195–1206, 2011.
  * Bindoff et al. [2007] N. L. Bindoff, J. Willebrand, V. Artale, A. Cazenave, J. M. Gregory, S. Gulev, K. Hanawa, C. Le Quere, S. Levitus, Y. Nojiri, et al.  _Climate change 2007-the physical science basis: Working group I contribution to the fourth assessment report of the IPCC_ , chapter Observations: oceanic climate change and sea level.  Cambridge University Press, 2007.
  * Bonev et al. [2018] B. Bonev, J. S. Hesthaven, F. X. Giraldo, and M. A. Kopera.  Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction.  _J. Comput. Phys._ , 362:425–448, 2018.
  * Borrero et al. [2004] J. C. Borrero, M. R. Legg, and C. E. Synolakis.  Tsunami sources in the southern california bight.  _Geophysical Research Letters_ , 31(13), 2004.
  * Bulleri and Chapman [2010] F. Bulleri and M. Chapman.  The introduction of coastal infrastructure as a driver of change in marine environments.  _J. Appl. Ecol._ , 47, 2010.
  * Carrier et al. [2003] G. F. Carrier, T. T. Wu, and H. Yeh.  Tsunami run-up and draw-down on a plane beach.  _Journal of Fluid Mechanics_ , 475:79–99, 2003.
  * [18] J. Carson.  Model verification and validation.  In _Proceedings of the Winter Simulation Conference_. IEEE.  doi: 10.1109/wsc.2002.1172868.  URL https://doi.org/10.1109/wsc.2002.1172868.
  * Chen et al. [2003] C. Chen, H. Liu, and R. C. Beardsley.  An unstructured grid, finite-volume, three dimensional, primitive equations ocean model: application to coastal ocean and estuaries.  _J. Atmos. Ocean Technol._ , 20:159–186, 2003.
  * Chen et al. [2014] C. Chen, Z. Lai, R. C. Beardsley, J. Sasaki, J. Lin, H. Lin, R. Ji, and Y. Sun.  The March 11, 2011 Tohoku M9.0 earthquake-induced tsunami and coastal inundation along the Japanese coast: A model assessment.  _Progress in Oceanography_ , 123:84–104, 2014.
  * Church and White [2006] J. A. Church and N. J. White.  A 20th century acceleration in global sea-level rise.  _Geophysical research letters_ , 33(1), 2006.
  * Clague [1997] J. J. Clague.  Evidence for large earthquakes at the cascadia subduction zone.  _Reviews of Geophysics_ , 35(4):439–460, 1997\.
  * Clawpack Development Team [2020] Clawpack Development Team.  Clawpack software, 2020.  URL http://www.clawpack.org.  Version 5.7.1.
  * Dalton et al. [2013] M. M. Dalton, P. W. Mote, and A. K. Snover.  _Climate Change in the Northwest_.  Springer, 2013.
  * Dean and Dalrymple [2002] R. G. Dean and R. A. Dalrymple.  _Coastal processes with engineering applications_.  Cambridge University Press, 2002.
  * Delestre et al. [2013] O. Delestre, C. Lucas, P.-A. Ksinant, F. Darboux, C. Laguerre, T.-N. Vo, F. James, S. Cordier, et al.  Swashes: a compilation of shallow water analytic solutions for hydraulic and environmental studies.  _International Journal for Numerical Methods in Fluids_ , 72(3):269–300, 2013.
  * Dugan and Hubbard [2010] J. Dugan and D. Hubbard.  Ecological effects of coastal armoring: A summary of recent results for exposed sandy beaches in southern california.  In H. Shipman, M. Dethier, G. Gelfenbaum, K. Fresh, and R. Dinicola, editors, _Puget Sound Shorelines and the Impacts of Armoring, U.S. Geol. Surv. Sci. Invest. Rep._ , 2010.
  * Fauzi and Mizutani [2020] A. Fauzi and N. Mizutani.  Machine learning algorithms for real-time tsunami inundation forecasting: A case study in Nankai Region.  _Pure and Applied Geophysics_ , 177:pages1437–1450, 2020\.
  * Fourestey et al. [2014] G. Fourestey, B. Cumming, L. Gilly, and T. C. Schulthess.  First experiences with validating and using the Cray Power Management Database tool.  _arXiv_ , 2014.  URL http://arxiv.org/abs/1408.2657.
  * Freitag et al. [2011] B. Freitag, J. Wiser, A. Engstfeld, K. Killebrew, C. Scott, R. Kasprisin, T. DeMarco, J. Vitulli, O. El-Anwar, K. Hochstatter, J. Schelling, D. Nelson, J. Mooney, B. Walker, T. Walsh, T. Biasco, N. Wood, F. González, T. Wilde, S. Fritts, D. Rowlett, C. Nelson, L. Shipman, and G. Miles.  Project safe haven: Tsunami vertical evacuation on the washington coast, pacific county.  Technical report, University of Washington and Washington Emergency Management Division, 2011.
  * Fuhrer et al. [2018] O. Fuhrer, T. Chadha, T. Hoefler, G. Kwasniewski, X. Lapillonne, D. Leutwyler, D. Lüthi, C. Osuna, C. Schär, T. C. Schulthess, and H. Vogt.  Near-global climate simulation at 1 km resolution: establishing a performance baseline on 4888 GPUs with COSMO 5.0.  _Geosci. Model Dev._ , 11:1665–1681, 2018.
  * Galvez et al. [2014] P. Galvez, J.-P. Ampuero, L. Dalguer, S. S.N., and T. Nissen-Meyer.  Dynamic earthquake rupture modelled with an unstructured 3-D spectral element method applied to the 2011 M9 Tohoku earthquake.  _Geophys. J. Int._ , 198:1222–1240, 2014.
  * Giles et al. [2021] D. Giles, D. Gopinathan, S. Guillas, and F. Dias.  Faster than real time tsunami warning with associated hazard uncertainties.  _Frontiers in Earth Science_ , 8, 2021.  doi: 10.3389/feart.2020.597865.
  * [34] Google.com.  Google tpu.  URL https://cloud.google.com/tpu/docs/tpus.
  * Gordon [2012] C. Gordon.  Analytical baseline study for the cascadia earthquake and tsunami.  Technical report, National Infrastructure Simulation and Analysis Center, Department of Homeland Security, 2012.  URL https://mil.wa.gov/asset/609ebfd1b819f.
  * Graham and Diaz [2001] N. E. Graham and H. F. Diaz.  Evidence for intensification of North Pacific winter cyclones since 1948\.  _Bulletin of the American Meteorological Society_ , 82(9):1869–1894, 2001.
  * Grothe et al. [2011] P. R. Grothe, L. A. Taylor, B. W. Eakins, K. S. Carignan, R. J. Caldwell, E. Lim, and D. Z. Friday.  Digital elevation models of crescent city, california : procedures, data sources, and analysis.  Technical report, National Geophysical Data Center, Marine Geology and Geophysics Division, 2011.  URL https://repository.library.noaa.gov/view/noaa/1188.
  * Heaton and Hartzell [1987] T. H. Heaton and S. H. Hartzell.  Earthquake hazards on the Cascadia subduction zone.  _Science_ , 236:162–169, 1987.
  * Horrillo et al. [2015] J. Horrillo, S. Grilli, D. Nicolsky, V. Roeber, and J. Zhang.  Performance benchmarking tsunami models for NTHMP’s inundation mapping activities.  _Pure Appl. Geophys._ , 172:869–884, 2015.
  * Hu et al. [2022] L. R. Hu, D. Pierce, Y. Shafi, A. Boral, V. Anisimov, S. Nevo, and Y.-F. Chen.  Accelerating physics simulations with tensor processing units: An inundation modeling example.  _The International Journal of High Performance Computing Applications_ , 36(4):510–523, 2022.  doi: 10.1177/10943420221102873.
  * ISEC [2004] ISEC.  Benchmark problem 1 tsunami runup onto a plane beach in the third international workshop on long-wave runup models.  Website, 2004.  http://isec.nacse.org/workshop/2004_cornell/bmark1.html.
  * Isozaki and Unoki [1964] I. Isozaki and S. Unoki.  The numerical computation of the tsunami in tokyo bay caused by the chilean earthquake in may, 1960.  _Studies on Oceanography_ , pages 389–402, 1964.
  * Jiang and Shu [1996] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  _Journal of Computational Physics_ , 126(1):202–228, June 1996.  doi: 10.1006/jcph.1996.0130.  URL https://doi.org/10.1006/jcph.1996.0130.
  * Jouppi et al. [2017] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, and A. Borchers.  In-datacenter performance analysis of a tensor processing unit.  In _ACM/IEEE 44th Ann. Int. Symp. on Comp. Architecture (ISCA). IEEE_ , pages 1–12, 2017.
  * Kamiya et al. [2022] M. Kamiya, Y. Igarashi, M. Okada, and T. Baba.  Numerical experiments on tsunami flow depth prediction for clustered areas using regression and machine learning models.  _Earth, Planets and Space_ , 74(127), 2022.
  * Kennedy et al. [2000] A. Kennedy, Q. Chen, J. Kirby, and R. Dalrymple.  Boussinesq modeling of wave transformation, breaking and runup, part I: 1D.  _J. Waterw. Port Coast. Ocean Eng._ , 126:39–47, 2000.
  * Kim and Lynett [2011] D. Kim and P. Lynett.  Turbulent mixing and passive scalar transport in shallow flows.  _Phys. Fluids_ , 23:016603, 2011.
  * Knudson and Bettinardi [2013] B. Knudson and A. Bettinardi.  Estimated economic impact analysis due to failure of the transportation infrastructure in the event of a 9.0 Cascadia Subduction Zone Earthquake, 2013.
  * Komar [1998] P. Komar.  _Beach processes and sedimentation_.  Prantice Hall, 1998.
  * Leschka and Oumeraci [2014] S. Leschka and H. Oumeraci.  Solitary waves and bores passing three cylinders-effect of distance and arrangement.  _Coast. Eng. Proc_ , 39, 2014.
  * LeVeque et al. [2011] R. LeVeque, D. George, and M. Berger.  Tsunami modelling with adaptively refined finite volume methods.  _Acta Numerica_ , 20:211–289, 2011.
  * LeVeque [2011] R. J. LeVeque.  _Finite volume methods for hyperbolic problems_.  Cambridge Univ. Press, 2011.
  * Liu et al. [2021] C. M. Liu, D. Rim, R. Baraldi, and R. J. LeVeque.  Comparison of machine learning approaches for tsunami forecasting from sparse observations.  _Pure and Applied Geophysics_ , 178:5129–5153, 2021.
  * Liu et al. [1994] X. Liu, S. Osher, and T. Chan.  Weighted essentially non-oscillatory schemes.  _Journal of Computational Physics_ , 115(1):200–212, 1994.
  * López-Venegas et al. [2014] A. López-Venegas, J. Horrillo, A. Pampell-Manis, V. Huérfano, and A. Mercado.  Advanced tsunami numerical simulations and energy considerations by use of 3d-2d coupled models: The october 11, 1918, Mona Passage tsunami.  _Pure App. Geophys._ , 171:2863–3174, 2014.
  * Løvholt et al. [2019] F. Løvholt, S. Lorito, J. Macías, M. Volpe, J. Selva, and S. Gibbons.  Urgent tsunami computing.  In _2019 IEEE/ACM HPC for Urgent Decision Making (UrgentHPC)_ , pages 45–50, 2019.  doi: 10.1109/UrgentHPC49580.2019.00011.
  * Lu et al. [2020a] T. Lu, Y. Chen, B. Hechtman, T. Wang, and J. Anderson.  Large-scale discrete fourier transform on tpus.  _arXiv:2002.03260_ , 2020a.
  * Lu et al. [2020b] T. Lu, Y. Marin, T. Zhuo, Y. Chen, and C. Ma.  Accelerating mri reconstruction on tpus.  In _2020 IEEE High Performance Extreme Computing Conference (HPEC)_ , pages 1–9, 2020b.
  * Lunghino et al. [2020] B. Lunghino, A. Santiago Tate, M. Mazereeuw, A. Muhari, F. Giraldo, S. Marras, and J. Suckale.  The protective benefits of tsunami mitigation parks and ramifications for their strategic design.  _Proc. Nat. Acad. Sci. (PNAS)_ , 1911857117, 2020.
  * Lynett et al. [2002] P. Lynett, T. Wu, and L. P.L.F.  Modeling wave runup with depth-integrated equations.  _Coast. Eng._ , 46:89–107, 2002.
  * Lynett [2007] P. J. Lynett.  Effect of a shallow water obstruction on long wave runup and overland flow velocity.  _Journal of Waterway, Port, Coastal, and Ocean Engineering_ , 133(6):455–462, 2007.
  * Ma et al. [2012] G. Ma, F. Shi, and J. Kirby.  Shock-capturing non-hydrostatic model for fully dispersive surface wave processes.  _Ocean Modeling_ , 43-44:22–35, 2012.
  * Macías et al. [2020a] J. Macías, M. Castro, and C. Escalante.  Performance assessment of the tsunami-hysea model for nthmp tsunami currents benchmarking. laboratory data.  _Coastal Engineering_ , 158:103667, 2020a.
  * Macías et al. [2020b] J. Macías, M. Castro, S. Ortega, and J. González-Vida.  Performance assessment of Tsunami-HySEA model for NTHMP tsunami currents benchmarking. field cases.  _Ocean Modelling_ , 152:101645, 2020b.
  * Macías et al. [2017] J. Macías, M. Castro, S. Ortega, C. Escalante, and J. González-Vida.  Performance Benchmarking of Tsunami-HySEA Model for NTHMP’s Inundation Mapping Activities.  _Pure Appl. Geophys._ , 174:3147–3183, 2017.
  * Mandli et al. [2016] K. T. Mandli, A. J. Ahmadia, M. Berger, D. Calhoun, D. L. George, Y. Hadjimichael, D. I. Ketcheson, G. I. Lemoine, and R. J. LeVeque.  Clawpack: building an open source ecosystem for solving hyperbolic pdes.  _PeerJ Computer Science_ , 2:e68, 2016.  doi: 10.7717/peerj-cs.68.
  * Mao et al. [2020] Z. Mao, A. D. Jagtap, and G. E. Karniadakis.  Physics-informed neural networks for high-speed flows.  _Computer Methods in Applied Mechanics and Engineering_ , 360:112789, 2020.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2019.112789.
  * Marras and Mandli [2021] S. Marras and K. T. Mandli.  Modeling and simulation of tsunami impact: A short review of recent advances and future challenges.  _Geosciences_ , 11:5–, 2021.
  * Marsooli and Wu [2014] R. Marsooli and W. Wu.  Numerical investigation of wave attenuation by vegetation using a 3D RANS model.  _Advances in Water Resources_ , pages 245–257, 2014.
  * Maza et al. [2015] M. Maza, J. Lara, and I. Losada.  sunami wave interaction with mangrove forests:a 3-D numerical approach.  _Coast. Eng._ , 98, 2015.
  * Mukherjee et al. [2023] A. Mukherjee, J. Cajas, G. Houzeaux, O. Lehmkuhl, J. Suckale, and S. Marras.  Forest density is more effective than tree rigidity at reducing the onshore energy flux of tsunamis.  _Coastal Engineering (Accepted for publication. In press)_ , 2023\.
  * Nelson et al. [1995] A. R. Nelson, B. F. Atwater, P. T. Bobrowsky, L.-A. Bradley, J. J. Clague, G. A. Carver, M. E. Darienzo, W. C. Grant, H. W. Krueger, R. Sparks, et al.  Radiocarbon evidence for extensive plate-boundary rupture about 300 years ago at the Cascadia subduction zone.  _Nature_ , 378(6555):371–374, 1995.
  * NOAA National Geophysical Data Center [2010] NOAA National Geophysical Data Center.  Crescent city, california 1/3 arc-second navd 88 coastal digital elevation model, 2010.  type : dataset.
  * Oishi et al. [2015] Y. Oishi, F. Imamura, and D. Sugawara.  Near-field tsunami inundation forecast using the parallel TUNAMI-N2 model: Application to the 2011 Tohoku-Oki earthquake combined with source inversions.  _Geo. Res. Letters_ , 42(4):1083–1091, 2015.
  * Park et al. [2013] H. Park, D. T. Cox, P. J. Lynett, D. M. Wiebe, and S. Shin.  Tsunami inundation modeling in constructed environments: A physical and numerical comparison of free-surface elevation, velocity, and momentum flux.  _Coastal Engineering_ , 79:9–21, 2013.
  * Pelties et al. [2012] C. Pelties, P. de la Punte, J.-P. Ampuero, G. Brietzke, and M. Käser.  Three-dimensional dynamic rupture simulation with a high-order discontinuous Galerkin method on unstructured tetrahedral meshes.  _J. Geophys. Res._ , 117:2156–2202, 2012.
  * Petersen et al. [2002] M. D. Petersen, C. H. Cramer, and A. D. Frankel.  Simulations of seismic hazard for the Pacific Northwest of the United States from earthquakes associated with the Cascadia subduction zone.  In _Earthquake Processes: Physical Modelling, Numerical Simulation and Data Analysis Part I_ , pages 2147–2168. Springer, 2002.
  * Peterson and Lowe [2009] M. Peterson and M. Lowe.  Implications of Cumulative Impacts to Estuarine and Marine Habitat Quality for Fish and Invertebrate Resources.  _Rev. Fish. Sci._ , 17:505–523, 2009.
  * Prasetyo et al. [2019] A. Prasetyo, T. Yasuda, T. Miyashita, and N. Mori.  Physical modeling and numerical analysis of tsunami inundation in a coastal city.  _Front. Built Environ._ , 5:46–, 2019.
  * Rasp et al. [2018] S. Rasp, M. S. Pritchard, and P. Gentine.  Deep learning to represent subgrid processes in climate models.  _PNAS_ , 115(39):9684–9689, 2018.
  * Roelvink and Van Banning [1995] J. Roelvink and G. Van Banning.  Design and development of DELFT3D and application to coastal morphodynamics.  _Oceanographic Lit. Review_ , 42:925–, 1995.
  * Ruggiero [2008] P. Ruggiero.  Impacts of climate change on coastal erosion and flood probability in the US Pacific Northwest.  In _Solutions to Coastal Disasters 2008_ , pages 158–169. 2008.
  * Ruggiero [2013] P. Ruggiero.  Is the intensifying wave climate of the US Pacific Northwest increasing flooding and erosion risk faster than sea-level rise?  _Journal of Waterway, Port, Coastal, and Ocean Engineering_ , 139(2):88–97, 2013.
  * Ruggiero et al. [2010] P. Ruggiero, P. D. Komar, and J. C. Allan.  Increasing wave heights and extreme value projections: The wave climate of the US Pacific Northwest.  _Coastal Engineering_ , 57(5):539–552, 2010.
  * Ruggiero et al. [2013] P. Ruggiero, M. G. Kratzmann, E. A. Himmelstoss, D. Reid, J. Allan, and G. Kaminsky.  National assessment of shoreline change: historical shoreline change along the pacific northwest coast.  Technical report, US Geological Survey, 2013.
  * Satake et al. [1996] K. Satake, K. Shimazaki, Y. Tsuji, and K. Ueda.  Time and size of a giant earthquake in Cascadia inferred from Japanese tsunami records of January 1700.  _Nature_ , 379(6562):246–249, 1996.
  * Satria et al. [2012] M. T. Satria, B. Huang, T.-J. Hsieh, Y.-L. Chang, and W.-Y. Liang.  GPU acceleration of tsunami propagation model.  _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , 5(3):1014–1023, 2012.  doi: 10.1109/JSTARS.2012.2199468.
  * Shi et al. [2012] F. Shi, J. Kirby, J. Harris, J. Geiman, and S. Grilli.  A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation.  _Ocean Model._ , 43-44:36–51, 2012.
  * Stoker [1957] J. J. Stoker.  page 333–341.  Interscience Publishers, 1957.
  * Takada and Atwater [2004] K. Takada and B. F. Atwater.  Evidence for liquefaction identified in peeled slices of Holocene deposits along the lower Columbia River, Washington.  _Bulletin of the Seismological Society of America_ , 94(2):550–575, 2004.
  * Thacker [1981] W. C. Thacker.  Some exact solutions to the nonlinear shallow-water wave equations.  _Journal of Fluid Mechanics_ , 107(-1):499, June 1981.  doi: 10.1017/s0022112081001882.  URL https://doi.org/10.1017/s0022112081001882.
  * Titov et al. [2005] V. Titov, F. González, E. Bernard, M. C. Eble, H. Mofjeld, J. C. Newman, and A. Venturato.  Real-time tsunami forecasting: Challenges and solutions.  _Nat. Hazards_ , 35:41–58, 2005.
  * tsunamiTPUlab [2023] tsunamiTPUlab.  https://github.com/smarras79/tsunamitpulab/releases/tag/v1.0.0, 2023.
  * Ueno [1960] T. Ueno.  Numerical computations for the chilean earthquak tsunami.  _Oceanogr. Mag._ , 17(2):87–94, 1960.
  * Ulrich et al. [2019] T. Ulrich, S. Vater, E. Madden, J. Behrens, Y. van Dinther, I. van Zelst, J. Fielding, C. Liang, and A.-A. Gabriel.  Coupled, physics-based modeling reveals earthquake displacements are critical to the 2018 Palu, Sulawesi Tsunami.  _Pure Appl. Geophys._ , 176:4069–4109, 2019.
  * Wang et al. [2013] P.-L. Wang, S. E. Engelhart, K. Wang, A. D. Hawkes, B. P. Horton, A. R. Nelson, and R. C. Witter.  Heterogeneous rupture in the great Cascadia earthquake of 1700 inferred from coastal subsidence estimates.  _Journal of Geophysical Research: Solid Earth_ , 118(5):2460–2473, 2013.
  * Wessels et al. [2020] H. Wessels, C. Wei$\ss$enfels, and P. Wriggers.  The neural particle method – an updated lagrangian physics informed neural network for computational fluid dynamics.  _Computer Methods in Applied Mechanics and Engineering_ , 368:113127, 2020.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2020.113127.
  * Witter et al. [2003] R. C. Witter, H. M. Kelsey, and E. Hemphill-Haley.  Great Cascadia earthquakes and tsunamis of the past 6700 years, Coquille River estuary, southern coastal Oregon.  _Geological Society of America Bulletin_ , 115(10):1289–1306, 2003.
  * Xia and Liang [2018] X. Xia and Q. Liang.  A new efficient implicit scheme for discretising the stiff friction terms in the shallow water equations.  _Advances in Water Resources_ , 117:87–97, July 2018.  doi: 10.1016/j.advwatres.2018.05.004.  URL https://doi.org/10.1016/j.advwatres.2018.05.004.
  * Xing and Shu [2005] Y. Xing and C.-W. Shu.  High order finite difference WENO schemes with the exact conservation property for the shallow water equations.  _Journal of Computational Physics_ , 208(1):206–227, Sept. 2005.  doi: 10.1016/j.jcp.2005.02.006.  URL https://doi.org/10.1016/j.jcp.2005.02.006.

##  5 Appendix

###  5.1 Running the code

Due to the restrictions of the TPU using Google Cloud Storage, Google’s
buckets will need to be used in order to run the notebooks. With a computing
project setup on Google Cloud, users can quickly run any of the example
notebooks or design their own simulation. Any of the example notebooks
available on GitHub (with the exclusion of tpu_tsunami.ipynb, which contains
the full implementation with all of the different scenarios; and
Create_Scenarios.ipynb, which can aid users in generating a custom DEM file)
can be quickly ran by going through the notebook after a few early setup
steps.

  1. 1.

Download the TPU-Tsunami Repository from
https://github.com/smarras79/tsunamiTPUlab/releases/tag/v1.0.0 to your local
machine. Create a project on Google Cloud Platform and associate a publicly
available bucket with the project.

  2. 2.

Modify the user_constants.py file to specify the PROJECT_ID and BUCKET with
the specifics of your Google Cloud Project. If you wish to change some
simulation constants, modify the beginning of the tpu_simulation_utilities.py
file.

  3. 3.

Navigate to https://colab.research.google.com/ and open the example notebook
(or your own notebook) from the TPU-Tsunami Repository using Colab’s open from
Github tool.

  4. 4.

Navigate to Runtime ¿ Change runtime type, and verify that the TPU option is
chosen as the Hardware Accelerator.

  5. 5.

Upload your user_constants.py and tpu_simulation_utilities.py files to your
notebook session using the drag-and-drop feature under Files. Upload any
corresponding DEM files to the session as well.

  6. 6.

Specify a function corresponding to an initial condition for your DEM file (or
use one example initial condition).

  7. 7.

Set initial conditions, boundary conditions as clarified in the bottom of any
example notebook run. Set last simulation parameters defining numerical
resolution (resolution), time step size (dt), output file times, TPU core
configuration (currently only capable of variation of cy), and DEM file name
on bucket (dem_bucket_filename).

  8. 8.

Run the simulation.

  9. 9.

Analyze results.