Probabilistic quantification of tsunami current hazard using statistical emulation

In this paper, statistical emulation is shown to be an essential tool for the end-to-end physical and numerical modelling of local tsunami impact, i.e. from the earthquake source to tsunami velocities and heights. In order to surmount the prohibitive computational cost of running a large number of simulations, the emulator, constructed using 300 training simulations from a validated tsunami code, yields 1 million predictions. This constitutes a record for any realistic tsunami code to date, and is a leap in tsunami science since high risk but low probability hazard thresholds can be quantified. For illustrating the efficacy of emulation, we map probabilistic representations of maximum tsunami velocities and heights at around 200 locations about Karachi port. The 1 million predictions comprehensively sweep through a range of possible future tsunamis originating from the Makran Subduction Zone (MSZ). We rigorously model each step in the tsunami life cycle: first use of the three-dimensional subduction geometry Slab2 in MSZ, most refined fault segmentation in MSZ, first sediment enhancements of seabed deformation (up to 60% locally) and bespoke unstructured meshing algorithm. Owing to the synthesis of emulation and meticulous numerical modelling, we also discover substantial local variations of currents and heights.


Introduction
Following the unexpected damage incurred at ports from the tsunamis of 2004 (Indian Ocean), 2010 (Chile) and 2011 (Japan) [1,2], it is paramount to investigate the associated hazard. Despite recent studies [1,[3][4][5] and advances in high-fidelity modelling [6], probabilistic methods for quantification of future tsunami hazard due to strong flows in harbours are sparse [7,8]. The need for such a quantification is further accentuated by certain peculiarities related with the phenomena of large tsunami currents in ports, e.g. the drifting of the 285 m ship Maersk Mandraki on 26 December 2004 at the Omani port of Salalah [2], despite small wave heights. It is deceptive to associate high wave amplitudes with high velocities. The treacherous nature of the currents was aggravated by the fact that strong currents in harbours continued for hours after the waves with maximum amplitude had arrived. This is all the more consequential since conventional tsunami warnings may be lifted after visibly perceptible signs of the tsunami (i.e. vertical displacement) have disappeared, whereas the strong currents may manifest later on.
Probabilistic scenario-based tsunami hazard assessment (PTHA) delivers a priori critical data to buttress tsunami disaster planning and practice. Scenario-based assessment scores over its catalogue-based counterpart through a more comprehensive exploration of plausible scenarios. Probabilistic scenario-based assessment surpasses deterministic scenario-based assessment in its assignment of probabilities and weighed integration of the different plausible scenarios. There exist variants in the probabilistic methodologies employed in PTHA: Monte Carlo [9], logictree [10] and Bayesian [11]. For an in-depth discussion on PTHA, the reader is referred to the recent review of Grezio et al. [12]. However, apart from the difficulties in assigning probabilities to scenarios, the computational burden expended for simulating each scenario prohibits an exhaustive sweep over the entire range of plausible scenarios. In this work, we pursue another probabilistic route via statistical emulation to quantify uncertainties in future tsunamis due to the uncertain earthquake sources (see workflow in figure 1a). Since the probability of large magnitude events is small, a comprehensive coverage of the Gutenberg-Richter (G-R) relation requires a large number of runs for the diversity of plausible events to be well represented across magnitudes and source locations (at least thousands for a coarse quantification and orders more for realistic assessments). Due to the considerable computational complexity of each simulation of coastal tsunami currents, such a probabilistic endeavour is made feasible by essentially replacing the numerical tsunami model by a statistical surrogate: the emulator. To our knowledge, this is the first time that Gaussian process (GP) emulation has been marshalled to generate future earthquake-generated tsunami currents; it has been employed only once in the past for currents, for a single source of landslide-generated tsunamis with considerable benefits in terms of computational costs and hazard assessment [13]. Here, with a design of only 300 full-fledged training runs, we fit an emulator to rapidly predict the impact of 1 million plausible tsunamis at prescribed locations. These emulated runs enable us to characterize uncertainty in future tsunami currents. A recent work by Kotani et al. [14] adopts a similar strategy of approximating the input-output response surface, albeit using nonlinear regression. Zhang & Niu [15] showcase a comparable 1.38 million scenarios, although using linear combination of waves from unit sources. Another recent strategy for reduction of the number of tsunami simulations employs an event tree coupled with cluster filtering of sources [16,17].
Additionally, formidable computational challenges must be addressed in order to accurately represent both the actual geophysical processes and their uncertainties. Despite possible issues arising from handling fine resolutions, our main challenge lies in encapsulating a large number of these high-definition simulations within a statistical framework. This is an essential requirement for PTHA and stretches the limit of current high-performance computing (HPC) facilities, even with the latest graphics processing unit (GPU) acceleration [18]. Often, the trade-off between capability and capacity in HPC is left unresolved by either radically simplifying the physics (e.g. a linear tsunami propagation till say 100 m depth with the use of an empirical relationship thereafter), or running very few fine resolution simulations as scenarios. Given a validated tsunami model, we argue that our emulation framework, in this context of currents that are nonlinear and very sensitive to near shore bathymetry, attempts a solution to this trade-off between precision and coverage of uncertainties. It requires manipulation of very large datasets on HPC, as well as complex post-processing on diverse software and data platforms. Overall, this work pushes the boundaries of current state-of-the-art in quantifying port hazard-with multithreaded emulation platform for large-scale (1 million) predictions, built on 300 high-definition simulations on smart unstructured meshes (10 m), using massively parallel multi-GPU-enabled simulations of validated tsunami model VOLNA, and hierarchical file formats-all integrated in an overarching workflow. We illustrate the emulation framework for the Karachi port in the Makran Subduction Zone (MSZ).
The MSZ has given rise to tsunamis in 1524 [19], 1945 [20,21] and 2013 [22]. Recent studies estimate the mega-thrust potential for the eastern part of the MSZ (blue rectangle in figure 1b) to be M w 8.8-9.0 [23]. Thus, here is a pressing need for a comprehensive quantification of tsunami hazard, especially port velocities and associated uncertainties. However, the accurate simulation of tsunami currents at shallow depths requires accurate coastline definition and bathymetry, with adequately refined meshes over a long duration to capture the maximum. Thus, in this study, we employ spatial resolutions of 10 m for the computational mesh, 30 m for bathymetry and 10 m for coastline, locally in the vicinity of Karachi port, for a total simulation time of 12 h. Furthermore, we employ here an earthquake source designed with segments of size 5 × 5 km with carefully constructed positive slip kernels to preserve fidelity to both magnitude scaling [24] and slip scaling relations [25]. The presence of a considerable sediment layer over the MSZ demands incorporation of its influence on the seabed deformation, since an appreciable amplification of up to 60% can be generated [26]. Section 2 describes the models and methods, §3 details the emulation framework, §4 discusses the results and conclusions are drawn in §5.

Models, data and methods
In this section, we describe in §2a the MSZ, in §2b the finite fault (FF) apparatus and slip profile, in §2c integration of the sediment amplification over the slips for generating seabed deformation (or uplift) and in §2d tsunami propagation.

(a) MSZ
The MSZ is formed by the subduction of the Arabian plate under the overriding Eurasian plate. It extends approximately 900 km from the Ornach Nal fault (approx. 67 • E) in the east to the Minab-Zendan-Palami fault (approx. 52 • E) in the west [20,27,28]. The mega-thrust potential of the entire MSZ is estimated at M w 9.07-9.22 [23]. Constraints imposed by GPS data resulted in three major segments and an estimated approximately 58% mean coupling ratio between the plates [27]. The subduction interface is divided into the eastern and western MSZ, with the eastern half being more seismically active. Given the scope of this work, we limit ourselves to the eastern MSZ, since tsunamis from western MSZ would have less appreciable effects on Karachi port than those arising from the western MSZ. Furthermore, paleoseismic accounts hypothesize that the western MSZ is seismically inactive compared with the eastern MSZ [29,30].
Here, the probability distribution function (pdf) for the G-R relation is modelled as the doubly truncated exponential distribution [31] where β = b log e 10, and the lower M m w and upper M M w limits of truncation are 4 and 8.8, respectively. The upper limit of M w 8.8 derives from the mega-thrust potential of eastern MSZ [23]. The rate parameter b of 0.92 is taken from the recent Earthquake Model of Middle East database (see electronic supplementary material, table S2 in [32]), and refers to the whole MSZ. For the scope of this work, we assume it as representative of the eastern MSZ. The complementary cumulative distribution function (ccdf), also called probability of exceedance or survival function, is then: Two cases of the truncated G-R distributions are plotted in figure 2a, i.e. for maximum magnitudes M M w of 8.8 and 8.6. Figure 2b shows histograms of actual samples from the distribution (used later in this work).

(b) Finite fault and slip profile
A FF on the eastern section of MSZ (blue rectangle, figure 1b) is constructed using a total number (n F ) of 2295 rectangular segments. The overall dimension of the FF model is 420 × 129 km 2 (L max × W max ). The slip on a segment is denoted by S i , where i varies from 1 to 2295. Okada  superposition of vertical displacements due to all the activated fault segments. Among the FF parameters, the dip angle and fault depth (d f ) are sourced from the recent plate boundary model, Slab2 [34,35]. The strike and rake angles are kept constant at 270 • and 90 • . A segment size (h 2 s ) is approximately 5 × 5 km 2 (l i × w i ), and the segments are arrayed in an 85 × 27 grid. This segment size is not chosen arbitrarily. It is selected based on a numerical study of the fidelity of the segmentation viz. 5 × 5 km 2 , 10 × 10 km 2 and 20 × 20 km 2 (figure 15a) to the earthquake dimension-magnitude scaling relation [24] (figure 1c). The discrepancy to the scaling relation appears as discontinuities in the realizable fault lengths (L) and widths (W) (figure 15a, inset). The size of the discontinuities are ∼h s .
We use the definitions of the seismic moment M 0 = n F i=1 μl i w i S i and moment magnitude M w = (2/3)(log 10 M 0 − 9.1), with μ = 3 × 10 10 N m −2 being the modulus of rigidity. Our implementation of the Okada suite is adapted from the dMODELS 1 code [36,37]. Slips are usually modelled to be uniform on the FF segments, even though inversions of seismic sources evidence localized concentrations of high slips (asperities) over a backdrop of lower slips [12]. Appendix A details the construction of the non-uniform slip profile used in this work.

(c) Influence of sediment amplification on seabed deformation
Incorporation of the effect of sediments influences tsunami modelling mainly in two ways. First, the interplay of sediment transport and tsunami flow gives rise to enhanced coupled morph-and hydro-dynamics [38,39]. Second, the Okada deformation model [33], with the assumptions of an elastic, homogeneous, isotropic medium in a semi-infinite domain, can be improved by sediment models that exhibit nonlinear, non-homogeneous and an-isotropic behaviour. Considerable amplification (up to 60% locally) of crustal deformation due to the presence of layers of sediments on the seafloor can occur [26]. In this section, we limit the incorporation of the effect of sediments to the deformation model by making use of a sediment amplification curve (figure 3c), extracted from elastodynamic simulations of layered sediment-rock seabed [26]. The curve uses the relative  where d i s is the sediment thickness over the segment interpolated from GlobSed 2 [40], and d i f is the down-dip fault depth of the segment taken from Slab2 [34] (figure 3a). Given d i r , the sediment amplification curve supplies the sediment amplification factor (S i a ) on the segment (figure 3d). The amplification due to the sediments is incorporated by multiplying the slip S i with the sediment amplification factor S i a resulting in an effective slip S e i (figure 3e) Okada's closed-form equations transform the effective slips S e i into the effective vertical displacement U e (figure 4b) [33]. The influence of sediments not only increases the slips effectively but also modifies the profile, as evident in the emergence of a double-lobed profile (figure 3e). The effect is more conspicuous in the associated deformations (compare figure 4a,b). The amplification factor (S a ) peaks at a relative depth of approximately 0.13, after which it decreases. Given the geometry of the fault and overlying sediment profile, a significant number of segments have an amplification factor between 0.4 and 0.6 (or, equivalently 40-60% amplification) (figure 3c inset and d). Furthermore, the sediment amplification factor is strongly dominated by the fault depth rather than the sediment thickness, which is near-uniform. The sediment amplification curve is defined only till a relative depth of 0.23 [26]. We linearly extrapolate the curve in order to be as conservative as possible in the region where it is not defined as well as to smoothly transition from regions of higher to lower fault depths. The counterparts of average slip S avg and maximum slip S max of S (without sediments) are defined as average effective slip S e avg and maximum effective slip S e max of S e (with sediments). Similarly, effective moment magnitude M e w is defined, by replacing S i with S e i in the expression of M w . The effect of sediments on slips are compared in figure 5a. Here, the increased scatter of S e max compared with S e avg is due to the spatial distribution of S a , which significantly amplifies S e max depending on the epicentre (X o , Y o ). Also, the increase in scatter of S e max as M w decreases is due to the decrease in fault dimensions that allow many earthquake scenarios to be situated in areas of lower S a . This aspect is pronounced in a similar comparison of M e w with M w in figure 5b.

(d) Tsunami propagation
Analysing wave heights requires few hours of simulation, while investigating the velocities needs longer simulation times. Thus, each scenario is run for 12 h of simulation time T s to obtain the maximum tsunami velocity and wave height, and is therefore computationally expensive. It is not only imperative that the numerical algorithms in the computer code for tsunami simulations run efficiently at fine mesh resolutions (10 m) needed to capture the currents, but also that the code is amenable to adequate parallelization, e.g. [41,42]. Thus, to run 300 such scenarios, we employ VOLNA-OP2 3 that runs efficiently for unstructured meshes on parallel GPUs [18]. The number of full-fledged scenarios (i.e. 300) is considerably higher than in existing studies related to MSZ [43][44][45]. Usual simulations employ the Green's functions approach to superpose the tsunami wave heights from a multi-segment FF source. Here, the nonlinear shallow water equations model not only the propagation of the tsunami but also the run-up/down process at the coast [46]. The finite volume (FV) cell-centred method for tessellation of control volume is used in VOLNA, and the barycentres of the cells are associated with the degrees of freedom. Details of numerical implementation, validation against standard benchmarks and comprehensive error analysis are available [18,47]. An important factor affecting the fidelity of long-lasting simulation of currents is numerical dissipation. Giles et al. [48] studied the numerical errors in VOLNA-OP2, wherein they are analysed by decomposing them into dispersion and dissipation components. Furthermore, an inter-model benchmarking of different numerical models highlighted the pitfalls in high-resolution current simulations [6]. In line with the scope of this work, we limit our numerical studies using VOLNA-OP2. It may be noted that although the emulation framework is independent of the specific numerical model employed, the accuracy of the emulator is limited by the accuracy of the underlying numerical model. VOLNA models the tsunami life cycle with where H(x, t) = b + η is the total water depth defined as the sum of free surface elevation η(x, t), and time-dependent bathymetry b(x, t). The two horizontal components of the depth-averaged fluid velocity are in v(x, t), g is the standard gravity and I 2 is the 2 × 2 identity matrix. The maximum tsunami velocity v max and wave height η max at location x at time t are computed as  and The dynamic bathymetry b(x, t) is the sum of static bathymetry b s (x) and U e , the effective deformation due to the influence of sediments. Here, an instantaneous fault is assumed, i.e. U e is supplied once at the beginning of the simulation. Furthermore, to reduce the computational burden of calculating deformations from 300 events, U e is computed only within a uplift calculation box covering the fault (green box in figure 4).
Accurate bathymetry, precise coastline and good quality computational mesh are vital for a proper modelling of velocities and currents in shallow water and near the coast. Thus, b s uses GEBCO 2019 (15 resolution) [49] complemented with hydrographic charts for Karachi port (approx. 30 m resolution), and SRTM v3 topography (1 resolution) [  port structures and breakwaters along the coastline, Google Earth's satellite imagery (approx. 10 m resolution) is used. The merging is described in the electronic supplementary material. The non-uniform unstructured mesh is designed in three stages corresponding to three regions, viz. offshore, onshore and near the port. This three-pronged strategy strikes a balance between having a fine mesh resolution (10 m) near Karachi port and reducing the overall computational cost with approximately 2.64 × 10 6 triangles in total. The mesh is generated using Gmsh 4 [51]. The construction of the mesh is described in appendix B.
The outputs v max and η max for two training samples (nos. 1 and 129) are plotted in figures 6 and 7, respectively, alongside snapshots taken at various time instants during the simulation.

(a) Emulator construction
The numerical simulation of the tsunami life cycle, i.e. its generation, propagation and inundation at fine mesh resolutions is computationally expensive due to model nonlinearity, and typically consumes hours on supercomputers. This is all the more prohibitive for a probabilistic quantification where thousands of runs of the tsunami code are required to exhaust the range of plausible scenarios. Statistical surrogates (or emulators) provide a computationally cheap approximation of the complex tsunami solvers, together with estimates of uncertainties in the predictions. In this study, the three input model parameters are moment magnitude (M w ) and epicentre coordinates (X o , Y o ) (figure 1c, inset). The coordinates have their origin as the southwest corner of the MSZ. The inputs are transformed into effective seafloor deformation. The consequent tsunamis are propagated till Karachi port. The outputs of interest in our case are the maximum wave height (η max ) and maximum wave velocity (v max ) at n G (193) virtual gauge locations around the port.      Thus, the computer code (denoted by M) simulates a multi-physics two-stage physical model, i.e. from the input parameters (M w , X o , Y o ) to deformation U e , then from U e to tsunami outputs v max and η max . An essential stage is the creation of an informative dataset for constructing the emulator. This is also called the design of computer experiments and the dataset is termed as the training set. The specific purpose of the design stage is to capture the functional relationship between the input parameters (M w , X o , Y o ) and output quantities (η max , v max ) at a location. The Latin Hypercube Design (LHD) generates a set of points that are nearly uniformly spread to cover the input parameter space. Specifically, it maximizes the minimum distance between points in the set, a feature that explores the functional relationship better than a random scatter. In a physical sense, this spread of points endeavours to capture the information inherent in the input-output relationship as much as possible. The model is evaluated by computer runs of M at the training points. Here, we employ an LHD of size 300 for three parameters ( figure 8). This is large enough to capture complex nonlinear combined sensitivities to the input parameters (e.g. the influence of size and location in relatively small and mid-size events closer to Karachi, or large regional variations in spatial distributions of slips), but still fits within our computational budget. The GP emulator (denoted by M) interpolates across the input-output points in the training set. In other words, the constructed emulator works as an approximation of M, and can be used to generate predictions (or, evaluated) at any point in the space of input parameters. The predictions will be exact at the training points, but uncertain elsewhere. This uncertainty is modelled by a normal distribution whose mean and standard deviation are calculated using the Kriging formula (mean quantities denoted byv max andη max ) explicitly accounting for the design. This structure allows for any nonlinear relationship to be modelled with uncertainties dependent on the location of the design points, unlike in more standard linear or even nonlinear regressions where the structure is fixed a priori. Derivations and exact equations can be found in Beck & Guillas [52]. GP emulation has been instrumental in successfully quantifying uncertainties in tsunami heights generated by landslides over the North Atlantic and the Western Indian Ocean as well as earthquakes over Cascadia [13,[53][54][55]. We use the efficiently implemented multiple-output Gaussian process emulator (MOGP) 5 for emulation.
The covariance kernel is a key component in the construction of the emulator. Here, we use the Matern 5/2 kernel that is smooth enough to avoid a rough GP, but not extremely smooth thus being suitable for modelling the physics. The piecewise polynomial, rational quadratic, exponential and squared exponential functions are other candidates [56]. The parameters (or length scales) in the kernels and other hyperparameters are found via nonlinear optimization  Maximum velocity magnitudes (and heights) are positive. In order to respect this physical constraint and not predict negative velocities (and heights), we feed the logarithm of v max (and η max ) into the construction of the emulator. Since the constructed emulator is now in the logarithmic scale, we transform the predicted quantities back to the original scale by accounting for the lognormal nature of the predicted distributions. Hence, the confidence intervals for the predictions, representing uncertainties, are all rendered positive, and naturally skewed in that direction. Once the emulator is constructed, it needs to be validated before employing it for predictions.

(b) Emulator diagnostics
In order to validate the quality of the emulation, we provide Leave-one-out (L-O-O) diagnostics here. Our training set consists of 300 pairs of input-output quantities. In L-O-O, a reduced training set of 299 pairs is employed to build an emulator, which is then used to predict the output at inputs of the one pair that was left out. The predicted output (and its uncertainty) is compared with the actual output of the left out pair. This procedure is repeated 300 times to cover all the pairs in the training set. These tests are passed by the emulator, as seen for predicted v max in figure 9 andη max in figure 10. The comparison between the mean of predictions from the emulator M and the training data from the tsunami simulator M shows that the emulator approximates well the simulator. The vertical line segments connect the predicted mean with its counterpart in the training data. More importantly, the uncertainties in the predicted mean, quantified in the form of 90% prediction intervals (green bars in figures 9 and 10), represent well the uncertainties about these predictions (or are even slightly conservative), since around    In these cases, the location with respect to the port is such that negligible wave energy is radiated to the port. Conversely, the low M w events that do show appreciable velocities are located such that considerable wave energy reaches the port. The L-O-O also shows a decrease in this positional dependence as M w increases, due to an accompanying increase in fault area and energy. Additionally, numerical dissipation in the model does play a role here, and numerical schemes tailored for reducing numerical dissipation would increase the accuracy [48]. To be able to generate 1 million predictions, we employ MOGP. Once the predictions are finished, we are left with two histograms (one each forv max andη max ) at every virtual gauge, each made up of 1 million samples of predicted quantity. The histograms are processed to extract P e (I(x) ≥ I 0 ), the probability of exceedance. P e is the probability of the tsunami having I(x) ≥ I 0 at a gauge x. The intensity I is the measure of hazard, i.e. eitherv max orη max , and I 0 is the intensity threshold for the hazard quantity under consideration.

Results and discussion
We first plot the raw output from the 1 million predictions, i.e. the histograms at 193 gauges in figure 11a,b. At each gauge, two histograms are superimposed on each other. These correspond to the two G-R relations with varying maximum moment magnitude assumptions, i.e. M M w 8.6 and M M w 8.8 ( figure 2). The histograms also act as visual indicators for the measure of the hazard at the gauge, and will be cast as hazard maps in figures 13 and 14. Near the tip of breakwaters and the mouth of the harbour, we observe relatively higher velocities than in other regions. We also observe a complementary relation between the histograms of velocities and wave heights: the gauges having thicker histograms for velocity have thinner histograms for wave heights and vice versa. These phenomena can also be observed in the snapshots (compare figure 6b with 7b).
As expected, there is a clear reduction of hazard when the maximum moment magnitude is reduced. For closer inspection, we enlarge the normalized histograms at gauge no. 91 in figure 11c,d. Gauge no. 91 is located in the centre of the map near the mouth of the port and is chosen since there is substantial spread of both maximum velocities and wave heights in its histograms. In figure 11c, the normalized histograms for maximum velocity are plotted. The range of velocities for M w 8.8 extends till approximately 16 ms −1 , while it extends to only approximately 6.2 ms −1 for M w 8.6. Thus, we observe approximately 61% reduction in maximum velocity hazard for a M w 0.2 reduction in maximum moment magnitude. By comparison, for the same reduction in maximum moment magnitude, the reduction in hazard from maximum wave height is only approximately 38% (from approx. 4.5 to 2.8 m in figure 11d). The probability of exceedance P e that is extracted from the histograms is plotted in the inset of the respective figure. Figure 12a,b compare normalized histograms for 1 million (1 M) and 10 000 (10 k) samples of input parameters (figure 2b). The corresponding probability of exceedance P e plots with their 99% confidence intervals can be seen in the inset. In figure 12a, we observe that the histogram corresponding to 10 000 predictions is curtailed around 7.5 ms −1 and becomes very sparse for higher velocities. This is due to a deficit of samples that results in the isolated bars for higher velocities. This behaviour also translates into larger uncertainties (or wider confidence intervals) for estimates of low probabilities of P e . By contrast, 1 million predictions adequately sweep through the entire range of velocities resulting in lower uncertainties (  intervals) for the tail probabilities. It may be noted that tail probabilities in the P e curve correspond to extreme events with higher velocities. Similar behaviour is seen in figure 12b, where the deficit of samples is observed for maximum wave heights higher than 2.7 ms −1 for the case of 10 000 predictions.
In figure 12c,d, we plot the probability of exceedance curves extracted from the histograms of 1 million predictions for the 193 gauges. Superimposed on top are the P e curves for 10 000 predictions. The horizontal lines in the plots are the chosen values of probability of exceedance, 10 −1 , 10 −2 and 10 −3 , progressively decreasing by an order of magnitude. The vertical lines in figure 12c denote maximum velocities of 1.5, 3.1 and 4.6 ms −1 (or 3, 6 and 9 knots, respectively), values that demarcate categories of damage [5].  although the lower probabilities (around 10 −4 ) have been made accessible by 10 000 events, they require 1 million events for accurate resolution: with only 10 000 samples, both probabilities and quantities are overestimated between 10 −3 and 10 −4 . Hence, being able to produce a very large number of predictions is crucial to hazard assessment. Only with the utilization of the emulatorneeding only 300 simulations-are we able to afford realistic predictions of velocities and wave heights at high resolution. Port hazard is represented on maps by velocity zonations, a time-threshold metric and safe depths for vessel evacuation [3,5]. In this work, the probability of exceedance curves in figure 12 are cast as hazard maps [7,8].   velocities is concentrated at the tip of breakwaters and along the dredged channel leading into the port (seen in port bathymetry, electronic supplementary material), as also observed in Lynett et al. [4]. This is also supported by the patterns of localized higher maximum velocities in figure 6a,c. By contrast, the spatial distribution of P e for maximum wave height shows a complementary behaviour and is more spread out. Conversely, for chosen probabilities of exceedance, the corresponding hazard thresholds at the gauges are plotted in figure 14. As expected, the overall intensity thresholds increase with decrease in probability of exceedance. Again, the bulk of the maximum velocity threshold is concentrated at the tip of breakwaters and along the dredged channel ( figure 14a). Here too, we see a complementary behaviour for maximum wave height in figure 14b.
Velocities have more spatial variation than heights [57], and show increased sensitivity to port configurations, compared with wave heights [58]. The larger spatial variation of velocities in figure 12c compared with wave heights in figure 12d is evident in the probability of exceedance plotted for all the gauges. This can be attested in figure 11a,b, where the bulkiness of velocity histograms varies spatially much more than that of the heights. Additionally, at a given gauge, we observe that the spread of velocities is much more than those of the heights for the same set of earthquake scenarios, e.g. compare figure 11a,b for gauge no. 91. These behaviours can also be observed for individual runs from the spatial variations of maximum velocity and wave height, compare (a) and (c) in figure 6 to those of figure 7.
The probability of exceedance extracted in this work acts as the basic input for common hazard outputs of probability of occurrence (and return periods), especially the approximately 2475 year mean return period for the maximum considered tsunami as laid out in ch. 6 of ASCE 7-16 [59]. It also feeds into loss estimation functions [60]. Although a full/complete probabilistic description of hazard may remain elusive, a realistic goal of 'fullness' will be to carefully define and perform each step in the PTHA. In these terms, a 'full' probabilistic assessment would ideally need to include further sources of uncertainties, including a thorough analysis of the source uncertainties in its seismic and tectonic setting. These include layers of uncertainties that are either epistemic or aleatoric in nature. Epistemic uncertainties include the scaling relation, and the G-R approximation of the occurrence-magnitude relationship [61], i.e. both the maximum moment magnitude and the b-value. For MSZ, the major influence of the maximum magnitude was illustrated in an initial work [62], with a simplified tsunami modelling strategy. Here, we only assess two cases, for M M w 8.6 and M M w 8.8. Uncertainties in the near shore bathymetry also have a large influence on near shore hazard [63]. Furthermore, the entire MSZ needs to be modelled for an area-wide assessment of hazard at the major ports in Pakistan, Iran, Oman and India, while accounting for crustal, outer-rise and imbricate faults. Secondary tsunamigenic effects from earthquakes in the continental crust (submarine slumps and slides) need additional parameters, e.g. 27 November 1945 M w 8.1 [64], and 24 September 2013 M w 7.7 [22] events. Similarly, with appropriate additional parameters, outer-rise and splay faults can be incorporated into the source, e.g. barrier models. Although a large increase in the number of parameters (especially for spatial fields of parameters) presents a challenge to emulation, a solution presents itself in the combination of dimension reduction and emulation [63].
Aleatoric uncertainties in the variations of the geometry in the seafloor uplift and subsidence can be readily incorporated. An alternative to our slip profile generation is to directly parameterize the co-seismic deformation profile using three parameters (or more) [55]. The Okada model that transforms the slips to the vertical deformation is then bypassed. This route is quite attractive since it allows the creation of very realistic deformation patterns with a fixed number of parameters, and does away with the dependency of the deformation/slip on the resolution of the segmentation (shown in figure 15a, inset). Our work uniformly samples the 1 million samples for epicentre coordinates (another aleatoric uncertainty). However, a recent spatial distribution of locking has been made available for the MSZ [27]. It would be even more realistic to sample the epicentre coordinates using the locking distribution, since zones of high locking act as a major cause for earthquake recurrence, as recently hypothesized [65]. The locations could be further distributed based on the depth-dependent rigidity [66]. Randomness in tide levels at the time of impact (consequent changes of up to 25% reported [67]) could be included. A better approximation of the currents would be through threedimensional modelling that accounts for fluid behaviour of the vertical water column and variable vertical flow [6,68]. Better designs of computer experiments than the LHD could be employed to reduce uncertainties in the emulator's approximation, such as sequential design [52]. Instead of investigating a range of scenarios, if one only wants to examine the maximum wave height in order to build defences for instance, a recent surrogate-based optimization could be pursued whereby the design of the experiment is combined with a search for the maximum, saving large quantities of computational time and increasing accuracy due to the focus on the optimization [69]. To be able to emulate a sequence of multiple models of seabed deformation and tsunami propagation, and possibly a three-dimensional model of currents locally, a new approach, called integrated emulation, allows even better designs [70]. The most influential models are run more times where it matters, and the integrated emulator propagates uncertainties with higher fidelity by taking into account the intermediate models in the system of simulators. This approach has the potential to enable fully realistic end-to-end coupling of three-dimensional earthquake sources models with tsunami models [71].

Conclusion
In this paper, we provide a novel end-to-end quantification of uncertainties of future earthquakegenerated tsunami heights and currents in the MSZ: (i) We replace the complex, expensive high-resolution tsunami simulator by a functionally simple, cheap statistical emulator trained using 300 tsunami simulations at 10 m mesh resolution in the vicinity of the port. We propagate uncertainties from the G-R relation to tsunami impacts of maximum velocities and wave heights in the port area of Karachi, Pakistan. We observe maximum (extreme event) velocities and wave heights of up to 16 ms −1 and 8 m, respectively, for the range M w 7.5-8.8 ( figure 11). (ii) We perform the largest emulation using 1 million predictions/source scenarios. To our knowledge, this is the first large-scale uncertainty quantification of earthquake-generated tsunami current hazard. We are able to display the necessity of this very large number of predictions for resolving very low probabilities of exceedance (less than 10 −3 )-very high impact extreme events (v max > 7.5 ms −1 and η max > 3 m) with tighter uncertainties ( figure 12). (iii) We observe that reduction in hazard due to a reduction in maximum moment magnitude is more for velocities than wave heights. Near the mouth of the harbour, the reduction in hazard is approximately 61% for maximum velocity, but only approximately 38% for maximum wave height (corresponding to a reduction in maximum moment magnitude from 8.8 to 8.6) (figure 12c). (iv) We generate the first area-wide probabilistic hazard maps of tsunami currents from 1 million predicted scenarios at the Karachi port (figures 13a and 14a). It shows patterns that are geophysically meaningful and important for the next steps of disaster risk reduction. We identify concentrations of high probability of exceedance around the port for given intensity threshold (a maximum of approx. 18%, 10% and 4% for 3, 6 and 9 knots, respectively) ( figure 13a). Conversely, the same regions also have high intensity thresholds given the probability of exceedance (a maximum of approx. 3 15a). An h s ∼ 5 km gives 2295 segments for the overall FF dimensions of L max ∼ 420 km and W max ∼ 129 km. To resolve the slip profile adequately, we require a fault to span a minimum of four segments along both the length and width directions. Using the scaling relation for h s ∼ 5 km, this requirement gives a minimum M w 6.32 that can be accommodated on the FF. This is sufficient as our region of investigation starts at M min w = 7.5. The scaling relation also limits the maximum M w that can be accommodated on the FF area of L max × W max , giving M sat w = 8.65 (figure 1c  [25]) S max and S avg for a = 0.5 S max and S avg for a = 1 S max and S avg for a = 2 proportionately increase the slip on the maximum dimensions. Now, to generate the slip profile, a positive kernel function φ is used (figure 15b, inset): where the gamma function Γ enters the normalization constant, length scale r defines where φ is non-zero and α adjusts the steepness of φ. With φ as the core, the bi-lobed kernel Φ is defined as where r l and r r are the length scales of the left and right lobes, their values depending on the position of epicentre (X o , Y o ) with respect to fault length (L) and width (W). The tensor product of two bi-lobed kernels, one along the length and another along the width of the fault, yields the surface Φ ⊗ (figure 16): where [−r W , r E ] × [−r S , r N ] denotes the domain of the fault and r ⊗ = {r W , r E , r S , r N }, the distances of western, eastern, southern and northern sides of the fault rectangle from (X o , Y o ). A normalization of Φ ⊗ with the required moment magnitude yields the final slip profile S (e.g. figure 4). We select α = 1 by varying α to mirror the maximum slip S max and average slip S avg curves from empirical scaling relations (see table 2   training sample #1 f (x; r E , a = 1) f (x; r W , a = 1) F (x; r W , r E , a = 1) Here, λ o is 60 km, which is approximately 60% of the diagonal in the smallest fault, i.e. of size approximately 94 km × 34 km for a M w 7.5 event (sample no. 300). Assuming the time period to be the same everywhere, the wavelength λ n at depth b s (x) is found as [73] λ n b s (x) This relates the characteristic length of mesh triangle (or mesh size) h λ (b s ) at depth b s (x) as  The mesh sizing function h(π ) is based on the coast proximity π (x) (figure 17b inset), which is defined as the minimum distance of point x from the coastline C of the merged bathymetry b s π (x) = min The construction of h(π ) is split into three regions, viz. inundation, stretch and blow-up (figure 17b). In the inundation region, which is defined to extend inland for a distance π I (2.5 km) from the coast, the mesh size is prescribed as the minimum mesh size h m (500 m). Thus, the inundation region facilitates smooth transition between the onshore and offshore meshes. Further inland, we require the triangle sizes to explode quickly to the maximum mesh size h M (25 km). This region is called the blow-up region (from π S to π B in figure 17b). We introduce the stretch region between the end of the inundation region and the beginning of the blow-up region (i.e. from π I to π S in figure 17b), for a gradual transition of corresponding mesh sizes, i.e. from h m to h S (10 km). This gradual change is achieved by setting the size ratio ρ, which is the ratio of characteristic lengths of adjacent triangles (or grading gauge [74]) to 1.3. The stretch distance π S − π I is calculated as π S − π I = h m + ρh m + ρ 2 h m + . . . + ρ n S h m . ( B 7 ) Equation (B 7) is a geometric series that approximates the distance by summing up the sizes of n S + 1 triangles, lined up end-to-end in a straight line, monotonically increasing in size by a factor