Optimal placement of tsunami sensors with depth constraint

Shallow water equations

The equations in one dimension are obtained by taking a vertical cross-section (in this case, a y-plane cross-section) of some two-dimensional domain and neglecting flows not parallel to the cross-section. Assuming uniform fluid density, incompressible-irrotational flow, and frictionless bed, the set of equations in one dimension is

(1) $$h_t+(hu{)}_x=0,$$

(2) $$u_t+g{\eta }_x+u{u}_x=0,$$where h is the total water depth, u is the horizontal velocity, and g is the gravitational acceleration (throughout this study, the value of g is set to 9.81 m.s⁻²). We also use the variables d for water depth and η for the free surface elevation, giving the relation h = d + η. Here, Eq. (1) represents the conservation of mass, while Eq. (2) represents the momentum balance. Figure 1 illustrates the variables just described.

To extend the model to two dimensions, we will have two equations for momentum balance in x- and y-direction. We would then need one more variable: v to denote the horizontal velocity in the y-direction. We also introduce the advection terms vu_y and uv_x in the momentum-conservation-equations. Thus the 2D governing equations read as:

(3) $$h_t+(hu{)}_x+(hv{)}_y=0,$$

(4) $$u_t+g{\eta }_x+u{u}_x+v{u}_y=0,$$

(5) $$v_t+g{\eta }_y+v{v}_y+u{v}_x=0.$$

All the variables used in the two-dimensional SWE are illustrated in Fig. 2.

Finite volume method on a staggered grid

We apply the finite volume method on a staggered grid to solve the SWE numerically. In a staggered grid arrangement, η is stored in the cell centre, while vectors u and v are stored in the cell faces. We use a structured rectangular grid to partition the spatial domain into intervals of equal size Δx in one dimension, or rectangles of equal size Δx × Δy in two dimensions; integer indices (x_i, y_j) represent the cell centers while half-integer indices $(x_{i+{\displaystyle \frac{1}{2}}},y_{j+{\displaystyle \frac{1}{2}}})$ represent the cell faces. These are visualized in Figs. 3 and 4. We denote by t_n the discretization of the time dimension.

Note that computing the spatial derivatives of p = hu and q = hv is not as straightforward, as they are products of two variables stored using two different indices. We approximate the values using the following upwind scheme to compute h at the cell faces, based on the direction of the flow:

(6) $$\ast h_{i+{\displaystyle \frac{1}{2}}}=(\begin{array}{cc}{h}_i,& u_{i+{\displaystyle \frac{1}{2}}}>0,\\ h_{i+1},& u_{i+{\displaystyle \frac{1}{2}}}\le 0.\\ & \end{array}$$

We then define $\ast p_{i+{\displaystyle \frac{1}{2}}}=\ast h_{i+{\displaystyle \frac{1}{2}}}{u}_{i+{\displaystyle \frac{1}{2}}}$. Now, to compute ${\left(u{u}_x\right)}_{i+{\displaystyle \frac{1}{2}}}^n$, we take note of the following equality:

(7) $$u{u}_x={\displaystyle \frac{p}{h}{u}_x.}$$

This is a modification of the method earlier used by Magdalena, Erwina & Pudjaprasetya (2015) and is preferred here because the calculation of the advection term is simpler and less costly. Next, the value of h at half-integer indices and the value of p and u at integer indices are computed as

(8) $${\overline{h}}_{i+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{1}{2}\left(h_i+h_{i+1}\right),}$$

(9) $${\overline{p}}_i={\displaystyle \frac{1}{2}\left(p_{i+{\displaystyle \frac{1}{2}}}+p_{i-{\displaystyle \frac{1}{2}}}\right),}$$

(10) $$\ast u_i=(\begin{array}{cc}{u}_{i-{\displaystyle \frac{1}{2}}},& {\overline{p}}_i>0\\ u_{i+{\displaystyle \frac{1}{2}}},& {\overline{p}}_i\le 0.\\ & \end{array}$$

These quantities yield the following approximation of $(u{u}_x{)}_{i+{\displaystyle \frac{1}{2}}}^n$:

(11) $${\left(u{u}_x\right)}_{i+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{{\overline{p}}_{i+{\displaystyle \frac{1}{2}}}}{{\overline{h}}_{i+{\displaystyle \frac{1}{2}}}}\left({\displaystyle \frac{\ast u_{i+1}-\ast u_i}{\mathrm{\Delta }x}}\right),}$$

Thus, we obtain this discretization of the one-dimensional SWE:

(12) $${\displaystyle \frac{h_i^{n+1}-h_i^n}{\mathrm{\Delta }t}+{\displaystyle \frac{\ast p_{i+{\displaystyle \frac{1}{2}}}^n-\ast p_{i-{\displaystyle \frac{1}{2}}}^n}{\mathrm{\Delta }x}=0,}}$$

(13) $${\displaystyle \frac{u_{i+{\displaystyle \frac{1}{2}}}^{n+1}-u_{i+{\displaystyle \frac{1}{2}}}^n}{\mathrm{\Delta }t}+g{\displaystyle \frac{{\eta }_{i+1}^{n+1}-{\eta }_i^{n+1}}{\mathrm{\Delta }x}+{\left(u{u}_x\right)}_{i+{\displaystyle \frac{1}{2}}}^n=0.}}$$

Extending the finite volume method on SWE to two dimensions require us to compute ${\left(v{u}_y\right)}_{i+{\displaystyle \frac{1}{2},j}}^n$ ${\left(v{v}_y\right)}_{i,j+{\displaystyle \frac{1}{2}}}^n$ and ${\left(u{v}_x\right)}_{i,j+{\displaystyle \frac{1}{2}}}^n$. We only give the details in approximating the first expression, as the others can be done similarly. This time, we use the equalities

(14) $$v{u}_y={\displaystyle \frac{q}{h}{u}_y,v{v}_y={\displaystyle \frac{q}{h}{v}_y,andu{v}_x={\displaystyle \frac{p}{h}{v}_x.}}}$$

Analogous to the first advection term, we do an approximation for the following variables:

(15) $${\overline{h}}_{i+{\displaystyle \frac{1}{2},j}}={\displaystyle \frac{1}{2}\left(h_{i,j}+h_{i+1,j}\right),}$$

(16) $${\overline{q}}_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}={\displaystyle \frac{1}{2}\left(q_{i+1,j+{\displaystyle \frac{1}{2}}}+q_{i,j+{\displaystyle \frac{1}{2}}}\right),}$$

(17) $$\ast u_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}=(\begin{array}{cc}{u}_{i+{\displaystyle \frac{1}{2},j}},& {\overline{q}}_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}>0\\ u_{i+{\displaystyle \frac{1}{2},j+1}},& {\overline{q}}_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}<0.\end{array}$$

Thus, we have these approximations:

(18) $${\left(v{u}_y\right)}_{i+{\displaystyle \frac{1}{2},j}}={\displaystyle \frac{{\overline{q}}_{i+{\displaystyle \frac{1}{2},j}}}{{\overline{h}}_{i+{\displaystyle \frac{1}{2},j}}}\left({\displaystyle \frac{\ast u_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}-\ast u_{i+{\displaystyle \frac{1}{2},j-{\displaystyle \frac{1}{2}}}}}{\mathrm{\Delta }y}}\right),}$$

(19) $${\left(v{v}_y\right)}_{i,j+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{{\overline{q}}_{i,j+{\displaystyle \frac{1}{2}}}}{{\overline{h}}_{i,j+{\displaystyle \frac{1}{2}}}}\left({\displaystyle \frac{v_{i,j+{\displaystyle \frac{1}{2}}}-v_{i,j-{\displaystyle \frac{1}{2}}}}{\mathrm{\Delta }y}}\right),}$$

(20) $${\left(u{v}_x\right)}_{i,j+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{{\overline{p}}_{i,j+{\displaystyle \frac{1}{2}}}}{{\overline{h}}_{i,j+{\displaystyle \frac{1}{2}}}}\left({\displaystyle \frac{\ast v_{i+{\displaystyle \frac{1}{2},j+{\displaystyle \frac{1}{2}}}}-\ast v_{i+{\displaystyle \frac{1}{2},j-{\displaystyle \frac{1}{2}}}}}{\mathrm{\Delta }x}}\right),}$$and, finally, we obtain this discretization of the two-dimensional SWE:

(21) $${\displaystyle \frac{h_{i,j}^{n+1}-h_{i,j}^n}{\mathrm{\Delta }t}+{\displaystyle \frac{{}^{\ast }{p}_{i+{\displaystyle \frac{1}{2},j}}^n-{}^{\ast }{p}_{i-{\displaystyle \frac{1}{2},j}}^n}{\mathrm{\Delta }x}+{\displaystyle \frac{{}^{\ast }{q}_{i,j+{\displaystyle \frac{1}{2}}}^n-{}^{\ast }{q}_{i,j-{\displaystyle \frac{1}{2}}}^n}{\mathrm{\Delta }y}=0,}}}$$

(22) $${\displaystyle \frac{u_{i+{\displaystyle \frac{1}{2},j}}^{n+1}-u_{i+{\displaystyle \frac{1}{2},j}}^n}{\mathrm{\Delta }t}+g{\displaystyle \frac{{\eta }_{i+1,j}^{n+1}-{\eta }_{i,j}^{n+1}}{\mathrm{\Delta }x}+{\left(u{u}_x\right)}_{i+{\displaystyle \frac{1}{2},j}}^n+{\left(v{u}_y\right)}_{i+{\displaystyle \frac{1}{2},j}}^n=0,}}$$

(23) $${\displaystyle \frac{v_{i,j+{\displaystyle \frac{1}{2}}}^{n+1}-v_{i,j+{\displaystyle \frac{1}{2}}}^n}{\mathrm{\Delta }t}+g{\displaystyle \frac{{\eta }_{i,j+1}^{n+1}-{\eta }_{i,j}^{n+1}}{\mathrm{\Delta }y}+{\left(v{v}_y\right)}_{i,j+{\displaystyle \frac{1}{2}}}^n+{\left(u{v}_x\right)}_{i,j+{\displaystyle \frac{1}{2}}}^n=0.}}$$

Optimization problem

Suppose we have source points $R=\{r_1,\dots ,r_k\}$, where tsunami waves are generated. Let s be an arbitrary point in the spatial domain. We define

(24) $$\begin{array}{c}\tau \left(r_i;s\right):=\text{the}\text{time}\text{it}\text{takes}\text{for}\text{a}\text{wave}\text{generated}\text{from}\text{a}\text{source}\text{point}\\ r_i\text{to}\text{reach}\text{the}\text{point}s.\end{array}$$

To compute $\tau \left(r_i;s\right)$, we use the two-dimensional SWE presented previously. Suppose we wish to place l sensors, $\sigma =\{s_1,\dots ,s_l\}$. Since only one sensor is needed to be triggered for the tsunami to detected, the time it takes for the tsunami to be detected from a source r_i is given by

$$t(r_i;\sigma ):=\underset{1\le j\le l}{min}\tau \left(r_i;s_j\right).$$

Thus, the guaranteed tsunami detection time to reach at least one of the sensors $\sigma =\{s_1,\dots ,s_l\}$ from any of the possible sources is given by

(25) $$T(\sigma )=\underset{1\le i\le k}{max}t(r_i;\sigma ).$$

Suppose the sensors can be placed anywhere in the water domain D, then the goal is to minimize (25). This minimization problem was also used in previous studies (Ferrolino, Lope & Mendoza, 2020, Ferrolino et al., 2020; Astrakova et al., 2009).

In this work, we introduce a constraint by considering the scenario where the sensors can only be placed on locations with depth constraint, as in the case of near-shore tsunami sensors. If we denote by D* the points in D that are within the depth constraint (see Fig. 5), then the minimization problem we consider is given by

(26) $$(\begin{array}{c}\underset{\sigma }{min}T(\sigma ),\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\sigma \in D^{\ast }.\\ \end{array}$$

The constrained domain D* can be defined as D* = {(x, y) ∈ D: d(x, y) ≤ d_max }, where d_max is the maximum depth in which a sensor can be placed. The constrained optimization problem (26) can be reformulated as an unconstrained problem using the penalty method, that is, by considering the problem

(27) $$\underset{\sigma }{min}T(\sigma )+\mu {\chi }_{D\mathrm{\setminus }{D}^{\ast }}(\sigma ),$$where μ is the penalty parameter and ${\chi }_{D\mathrm{\setminus }{D}^{\ast }}$ (σ) is the characteristic function (or indicator function) of D\D*, whose value is 0 if σ ∈ D* and 1 otherwise. We set μ equal to a large number to penalize those sensor locations that do not satisfy the depth constraint. In our simulations, we set μ = 10⁶.

It should be noted that this formulation can be easily modified to consider cases when the tsunami does not originate from a point source (e.g., when the tsunami is due to a landslide). The minimization problem remains the same but the initial conditions of the SWE will have to be modified. Likewise, this approach can also be used to determine the location of deep-ocean tsunami sensors by either setting the depth constraint to a high value or the value of μ to 0.

Particle swarm optimization

As discussed previously, the function τ in (24) is computed by solving the two-dimensional SWE numerically. However, this only gives us the values of τ at points located at the cell centers of the rectangular grid. To calculate the travel time at other locations, we use bilinear interpolation using the values of τ at the surrounding points.

Bilinear interpolation is done by performing linear interpolation in one direction, and then performing it again in the other direction, on every rectangle [x_i, x_{i + 1}] × [y_j, y_{j + 1}]. For example, suppose we want to interpolate τ using some function f(x, y) on the unit square [0, 1] × [0, 1], given its values at the points (0, 0), (1, 0), (0, 1), and (1, 1). Then f(x, y) is given by

(28) $$\begin{array}{c}f(x,y)=f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy,\\ (x,y)\in {[0,1]}^2.\end{array}$$

An illustration of this is provided in Fig. 6.

Despite its name, as seen above, the resulting interpolant is not linear; rather, it is the sum of products of linear functions in x and y, respectively. Hence, the whole function is piecewise quadratic and may not be differentiable at the cell edges. Furthermore, the constrained optimization problem is transformed into an unconstrained problem using a penalty method. The penalty term in (27) includes a characteristic function, which is binary and not differentiable. In view of this, it is fitting to use optimization methods that do not rely on gradients. In this work, we explore the use of derivative-free metaheuristic algorithms in solving the minimization problem.

Metaheuristic algorithms have gained popularity because they only rely on function evaluations and have the capability of obtaining the global minimizer (Yang, 2014). It was shown in Ferrolino et al. (2020) how PSO was effective in solving the tsunami location problem without depth constraint. We adopt the same algorithm for our problem formulation.

The PSO algorithm starts by randomly generating a swarm of solution candidates (particles) with their initial velocity vector. Each particle then move around the search-space according to three parameters: its velocity vector, the vector pointing at its own best-known position, and the vector pointing at the swarm’s best-known solution. For each iteration, the step that each particle takes is a linear combination of these three vectors (see Fig. 7), expecting that the swarm eventually gathers around one point—the optimal solution.

The pseudocode of the algorithm is presented in Algorithm 1. The algorithm starts by generating swarm_size particles randomly in the search space D, and each their own initial velocity vector v. They all then follow a series of instructions based on computations provided in Algorithm 1, where the vectors p (each particle’s best known location) and g (the entire swarm’s best known location) is updated. This process terminates when either the maximum number of iterations (max_iter) is reached, or the change in the objective value f is less than min_func. In this study, the function f is the one presented in Eq. (26); and the value of each hyperparameters are swarm_size = 100, max_iter = 100, min_func = 10⁻⁸, ω = ϕ_p = ϕ_g = 0.5.

More detailed discussions of this method may be found in Kennedy & Eberhart (1995); Shi & Eberhart (1998), Mezura-Montes & Coello (2011). Our implementation of the PSO algorithm made use of the built-in command PySwarms, an extensible research toolkit for particle swarm optimization in Python.

1D standing wave on a flat bottom

A standing wave, simply put, is a wave that does not travel (horizontally). Rather, it is an oscillating wave fixed in space. At any given point in space, the peak amplitude of wave oscillations on a standing wave is constant. This phenomenon can be observed on the surface of a liquid in a vibrating container.

To produce a standing wave simulation in one dimension, we set up a container of length L, constant still water depth d; give it initial conditions $\eta (x,0)=\cos ({\displaystyle \frac{\pi x}{L})}$, u(x, 0) = 0 and wall boundary conditions u(0, t) = u(L, t) = 0. The result of the simulation on a flat bottom with length L = 25 m and depth of 4 m, using Δx = 0.2 m, is shown in Fig. 8. The wave seemingly does not travel horizontally, as a result of interference between two waves traveling in opposite directions. The water level at any point c oscillates in time with amplitude |η(c, 0)|. The analytical solution to this phenomenon is equivalent to that of a string of infinite length with some wave number k and angular frequency ω. In linear and non-dispersive case here, these numbers are just $k={\displaystyle \frac{\pi }{L}}$ and $\omega =k\sqrt{gd}$, which gives

(29) $$\eta (x,t)=\cos \left({\displaystyle \frac{\pi }{L}x}\right)\cos \left({\displaystyle \frac{\pi }{L}\sqrt{gd}t}\right).$$

Simulating this at Δx = 0.2 m, Δt = 0.032 s for 20 s gives a mean absolute error of 0.015 m. To simplify things, all simulations below assume hard wall boundary conditions and no current, unless stated otherwise.

1D wave run-up on a sloping beach

The first nonlinear simulation is wave run-up, which refers to either the phenomenon or the measure of when an ocean wave reaches a beach or sea dike structure and rises above still water level. Obviously, as a wave approaches the beach structure, the water depth gets smaller relative to the amplitude, hence the advection term becomes more significant. Furthermore, in this simulation the wave interacts with the dry region, something that was not observed in the previous simulations.

This simulation is based on an experiment by Synolakis (1987) on a solitary wave that propagates through a single slope. The solitary wave centered at x = x₀ has an initial surface profile $\eta (x,0)={\displaystyle \frac{H}{d}sec{h}^2\left(\gamma \left(x-{\displaystyle \frac{L}{4}}\right)\right)}$ with $\gamma =\sqrt{{\displaystyle \frac{3H}{4d}}}$. This way, similar to the step bottom simulation, the right traveling soliton eventually changes shape, this time as it climbs the sloping beach. This behavior shown in Fig. 9, where the simulation was done with Δx = 0.1 m and Δt = 0.016 s. Several gauges were placed during the experiment to record the wave height every 10 s. We then compared these records to the simulation. As seen in Fig. 9, the model was able to predict the experiment to a pretty high accuracy.

1D dam break over a flat, dry bed

Dam break, as the name suggests, happens when the dam holding a reservoir breaks, causing a sudden rush of water flowing to the other side. The set up of this simulation is on a flat bed with 1-m-deep water on one side of the domain (reservoir), and completely dry on the other side:

(30) $$h(x,0)=(\begin{array}{cc}H,& x<x_0\\ 0,& x>x_0,\\ & \end{array}$$where x₀ represents the location of the barrier. Right after the simulation starts, the previously contained water will stream to the dry side and start to fill it, while the water level on the reservoir goes down. The result of the simulation on a 20 m-long, 1 m-deep bed with Δx = 0.1 m is provided in Fig. 10. To test the robustness of the numerical scheme, we compare it with the analytical solution given by Chanson (2008):

(31) $$h(x,t)=(\begin{array}{cc}H,& x<-ct\\ {\displaystyle \frac{H}{9}{\left({\displaystyle \frac{2-x}{ct}}\right)}^2,}& -ct\le x<(2c-1)t,\\ & \end{array}$$where $c=\sqrt{gH}$. At Δx = 0.1 m, Δt = 0.008 s, the numerical scheme gives a mean absolute error of 0.0034 m.

2D planar surface wave over a parabolic basin

So far the simulations were set on rectangular or trapezoidal shapes. To avoid bias due to the also rectangular shape of the grids, this simulation will be set on a parabolic basin. This way, not only do we have a curved bottom profile, all of the perimeter will be interacting with dry regions. The results of this simulation further confirms the performance of the wet-dry procedure.

The shape of the basin follows the parabolic function $d(x,y)=d_0\left(1-(x^2+y^2)/L^2\right)$ on an 8,000 m × 8,000 m domain, with L² = 8,000 m². The large parabolic shape of the domain is supposed to mimic that of a lake. We would like to see the motion of the water surface initialized by $\eta (x,y,0)={\displaystyle \frac{2a_0{d}_0}{L}({\displaystyle \frac{x}{L}-{\displaystyle \frac{\xi }{2L})}}}$. The analytical solution of the planar surface wave obtained by Thacker (1981) is

(32) $$\eta (x,y,t)={\displaystyle \frac{2\xi d_0}{L}({\displaystyle \frac{x}{L}\cos \omega t-{\displaystyle \frac{y}{L}\sin \omega t-{\displaystyle \frac{\xi }{2L}),\omega ={\displaystyle \frac{1}{L}\sqrt{2g{d}_0}.}}}}}$$

The parameters used in the simulation were d₀ = 1 m, ξ = 400 m using Δx = Δy = 40 m, and Δt = 4.5 s. As expected, the planar surface in the parabolic basin oscillates with period $T=2\pi /\omega =2\pi L/\sqrt{g{d}_0}$, which can be seen in Fig. 11. Quantitatively, the numerical solution agrees with the analytical solution with mean absolute error of 0.029 m.

Tsunami run-up onto a conical island

This experiment was done at the Coastal Engineering Research Center, Vicksburg, Mississippi, as part of a research (Briggs et al., 1995) on a tsunami that struck Babi Island, Indonesia in 1992. The shape of Babi Island is similar to a truncated cone, as seen in Fig. 12, and the experiment may serve as a benchmark test for fluid dynamic models in modeling interaction between waves and such structure. It was conducted on a 25 m × 30 m basin with a conical island at the center, and waves with an initial solitary wave-like profile were generated from one side of the tank. The resulting water surface elevations were measured by 22 gauges. Our numerical simulation was ran with Δx = Δy = 0.2 m. The results of the simulation were compared with data gathered at 4 of the 22 gauges: G6, G9, G16, and G22. As seen in Fig. 13, the model was able to replicate the experiment with great accuracy.

Tsunami run-up onto a complex 3D beach

We then move on to an experiment with a much more complex bathymetry, similar to that found in nature. The experiment was based on a tsunami that struck Okushiri Island, Japan, with a really high current, that resulted in an extreme run up at the tip of a very narrow gulley within a small cove at Monai. It was done at the Central Research Institute for Electric Power Industry (CRIEPI) in Abiko, Japan on a 1:400 laboratory model of Monai (Liu, Yeh & Synolakis, 2008) and the setup in shown in Fig. 14. Just like the previous experiment, a series of waves were generated from the deeper side of the water. The results of our simulation, with Δx = Δy = 0.025 m, were compared with the recorded height on the three gauges placed during the experiment. As can be observed in Fig. 15, the model was able to replicate the experiment with high accuracy.

Step bottom

The first simulation was done on a step bathymetry of size 1,000 m × 1,000 m with d = 20 m for x < 500 m, and d = 10 m elsewhere. Waves were generated from five evenly spaced points along the line x = 200 m, as seen in Fig. 16. The feasible region for this problem is {(x, y): d(x, y) ≤ 10 m}, which is the whole right side of the domain.

The result obtained is pretty intuitive, that is, the optimal location of detection sensor in this case is right in the middle of the line where the domain changes depth. If located far up on the y-axis, it will take long for the sensor to detect waves coming from source points further down (and vice versa). If placed too far to the right, waves from any source point will take a longer time to reach the sensor.

Sloping bottom

The next simulation was done on a sloping bathymetry of size 10,000 m × 10,000 m, with a slope of 0.1 and maximum depth of 1,000 m. This time the waves were generated through a flip-fault motion from six subduction zones in the shape of arcs of the semicircle $x+5000=-\sqrt{{5000}^2-{(y-5000)}^2}$, as seen in Fig. 17. The feasible region for this problem is {(x, y): d(x, y) ≤ d_max, x < −5,000}, with d_max = 900, 800, 700, 600 m. Since in this case the water depth gradually gets smaller as we approach the shoreline on the right side, it is interesting to investigate what happens if the depth constraints are varied.

As seen in Fig. 17, not all sensors go to their respective depth constraint. This is expected because the bottom topography is not flat. As predicted, the sensors obtained with varying depth constraints lie on the bisector of the semicircle. As the water depth increases, the location of the sensors move away from the shore and closer to the subduction zones. However, the sensor cannot get too close to the semicircle to guarantee detection time from all the possible sources. This can be seen on the obtained locations of sensors for depths 800 and 900, i.e., the two locations are almost identical. This is possible because the feasible region for d_max = 800 m is contained in the feasible region for d_max = 900 m.

Step bottom

Using the same setup as the previous one, the PSO algorithm was run to solve the optimal location of two detection sensors. As seen in Fig. 18, the obtained locations of the sensors are on the boundary of the feasible region. This is expected because the sensors must be placed as close as possible to the source points and at the same time must satisfy the depth constraint. It also makes sense that neither sensor was located right in the middle, since one sensor can account for waves coming from the upper side while the other can account for waves coming from the lower side.

Sloping bottom

In this example, the optimal locations of sensors are more spread out to account for waves coming from different angles. The results are illustrated in Fig. 19. Obviously, increasing the number of sensors will always decrease the wave detection time. As shown in Fig. 20, using two sensors instead of one allows the waves to be detected about 37 s earlier or 68% faster. However, using three or more sensors would be less efficient, as this only decreases the detection time by 1 s or even less.

Sloping bottom with an island

This simulation starts with the same setup as the sloping bottom simulation, only this time we introduce a small island in the middle of the domain. In this case, we expect to see interesting results on small depth constraints since there are two “shallow” areas here: the center island and the coastline on the right. It can be seen in Fig. 21 that sensors with depth constraints of 500, 700 m are not really affected by the island. However for maximum depth of 300 and 100 m, the pair of sensors are spread apart, with one approaching the island and the other getting closer to the coastline on the right. Because the domain is symmetric with respect to the line y = 5,000 m, two solutions are obtained. This is illustrated in Fig. 21. This configuration could not be tackled in Ferrolino et al. (2020) because the presence of an island in the simulation requires the use of the nonlinear SWE with wet-dry procedure.

The simulations we have presented so far were formulated so that the optimal solutions can be verified based on the geometry of the subduction zone and the domain. The obtained results using PSO yield solutions that agree with the expected geometric solutions. These tests were done to check that the method works. In a real-world scenario, the optimal placement of sensors is not trivial because the subduction zones do not necessarily follow a geometric pattern. Additionally, the bottom topography of oceans is uneven, making the computation of travel time more complicated. Because we use the 2D nonlinear SWE in simulating tsunami waves, our proposed method can handle complex bathymetric profiles and arbitrary tsunami sources. In our following and last example, we apply our tsunami sensor detection model to one of the more recent tsunami incidents: the 2018 Palu tsunami.