From offshore to onshore probabilistic tsunami hazard assessment via efficient Monte Carlo sampling

SUMMARY Offshore Probabilistic Tsunami Hazard Assessments (offshore PTHAs) provide large-scale analyses of earthquake-tsunami frequencies and uncertainties in the deep ocean, but do not provide high-resolution onshore tsunami hazard information as required for many risk-management applications. To understand the implications of an offshore PTHA for the onshore hazard at any site, in principle the tsunami inundation should be simulated locally for every earthquake scenario in the offshore PTHA. In practice this is rarely feasible due to the computational expense of inundation models, and the large number of scenarios in offshore PTHAs. Monte Carlo methods offer a practical and rigorous alternative for approximating the onshore hazard, using a random subset of scenarios. The resulting Monte Carlo errors can be quantified and controlled, enabling high-resolution onshore PTHAs to be implemented at a fraction of the computational cost. This study develops efficient Monte Carlo approaches for offshore-to-onshore PTHA. Modelled offshore PTHA wave heights are used to preferentially sample scenarios that have large offshore waves near an onshore site of interest. By appropriately weighting the scenarios, the Monte Carlo errors are reduced without introducing bias. The techniques are demonstrated in a high-resolution onshore PTHA for the island of Tongatapu in Tonga, using the 2018 Australian PTHA as the offshore PTHA, while considering only thrust earthquake sources on the Kermadec-Tonga trench. The efficiency improvements are equivalent to using 4–18 times more random scenarios, as compared with stratified-sampling by magnitude, which is commonly used for onshore PTHA. The greatest efficiency improvements are for rare, large tsunamis, and for calculations that represent epistemic uncertainties in the tsunami hazard. To facilitate the control of Monte Carlo errors in practical applications, this study also provides analytical techniques for estimating the errors both before and after inundation simulations are conducted. Before inundation simulation, this enables a proposed Monte Carlo sampling scheme to be checked, and potentially improved, at minimal computational cost. After inundation simulation, it enables the remaining Monte Carlo errors to be quantified at onshore sites, without additional inundation simulations. In combination these techniques enable offshore PTHAs to be rigorously transformed into onshore PTHAs, with quantification of epistemic uncertainties, while controlling Monte Carlo errors.


A.1 Proof that Equation 9 is the variance of Equation 7
As noted in the Section 3.2, Equation 7 can be written as: where B p b,i,T , N (M w,b ) is a binomial random variable with probability p b,i,T and size N (M w,b ). The binomial term has mean p b,i,T N (M w,b ) and variance N (M w,b )p b,i,T (1 − p b,i,T ) (e.g., Bolker, 2008). We wish to prove that Equation 9 is the variance of Equation A.1. Equation A.1 is a constant multiple of a binomial random variable. For any random variable X and constant k, if the variance of X is σ 2 (X) then the variance of kX is k 2 σ 2 (X). This follows directly from the definition of the variance (e.g. Bolker, 2008).
Applying these results to calculate the variance of Equation A.1 gives: which gives Equation 9.

A.2 Justification of Equation 10
Equation 10

B.1 Preliminary results
The following standard results are used in the subsequent arguments.
For a discrete random variable X with probability mass function ρ(x), the mean of X is: and the variance of X is: Bolker, 2008). Basic importance-sampling (e.g. Owen & Zhou, 2000;Tokdar & Kass, 2009;Lie & Quer, 2017) can be used to approximate Equation A.5 using a weighted random sample drawn from X with replacement, where the sampling weights correspond to a chosen probability mass function w τ (x), which can differ from ρ(x). Denote this sample as X τ , and assume it has sample size N . The basic importance-sampling estimate of µ(X) is: In our context, the offshore PTHA gives the within-magnitude-bin exceedance-rate as: gives the within-magnitude-bin conditional probability of scenario e (note for stratifiedsampling, w SS b,i is the same as the sampling weights, see Equation 6). If the scenarios are instead sampled with alternative weights w SIS b,i (e), then from Equation A.7 the basic importance-sampling estimate of Equation A.8 is: which is the same as Equation 17:

B.2 Proof that Equation 19 is the variance of Equation 17
Below we directly compute the variance of Equation and its variance is (using Equation A.6): Next suppose we sample N (M w,b ) scenarios. Consider the summation term on the RHS of Equation A.10. The variance of this sum is equal to the sum of the variances of the individual terms (because each sample is independent, Bolker, 2008): Finally, using the fact that for any random variable X and constant k it is true that σ 2 (kX) = k 2 σ 2 (X), the variance of Equation A.10 is: The last step gives Equation 19.

B.3 Justification of Equation 20
Equation 20 is an estimate of Equation 19 derived from basic importance-sampling. The summation term on the RHS of Equation A.21 can be approximated using the basic importance-sampling estimator (Equation A.7) as follows: • ρ(x) in Equations A.5 and A.7 is replaced with w SIS b,i (e) • x in Equations A.5 and A.7 would ideally correspond to [φ SIS b,i (e)1 (Q(e)>Q T ) − p b,i,T ] 2 . But it involves the term p b,i,T which cannot be computed using only the sampled scenarios.
-To work around this issue, in place of p b,i,T we use q b,i,T , which estimates the former using only the sampled scenarios. So , and X τ is replaced with E SIS b,i . These substitutions lead directly to Equation 20.

C Further information on the inundation model setup
The model domain (Figure 9) uses three levels of two-way grid nesting. The outer domain has 1850 m cell-size, while high resolution areas around Tongatapu have 7.5 m cell size. Elevation data was derived from a mixture of lidar, gridded bathymetric surveys near Tongatapu (Damlamian et al., 2013), and global-scale data derived from GA250 (Whiteway, 2009) and GEBCO-2014 Grid (version 20150318).
Tsunamis were simulated using the open-source SWALS code, which provides a number of explicit shallow water solvers including the linear solver used for PTHA18 (Davies & Griffin, 2018. This code has a test-suite of more than 20 analytical, laboratory or field problems, including well known tsunami problems from NTHMP (2012) (excluding landslides) and other recent NTHMP problems (Park et al., 2013;Lynett et al., 2017;Macías et al., 2020;Gao et al., 2020). The model setup herein mirrors that used by  to compare modelled and observed tsunamis at 16 Australian tide-gauges. On all nested grids a second-order finite-volume scheme is used to solve the full nonlinear shallow water equations with uniform Manning-friction of 0.03 (details similar to Davies & Roberts, 2015). On the outer domain the linear shallow water equations are combined with Manning-friction and solved using a classical scheme (Goto & Ogawa, 1997), with details matching .
To check the model performance, a number of convergence tests were conducted by factor-of-two coarsening of the model resolution. This had a small impact on the hydrodynamic results, suggesting the model resolution is fine enough. We also simulated five historic earthquake-tsunamis measured on the Nuku'alofa tide-gauge ( For each historic event a similar tsunami scenario was selected from among PTHA18 scenarios that fit DART buoy observations relatively well (as identified in Davies, 2019). To simulate the Tohoku and Chile events the outer domain was extended to cover the Pacific Ocean and it was evolved for 24 hours. For the Tonga events the model was evolved for 5 hours with an outer-domain matching Figure 9, identical to that used to simulate the random scenarios.