Probabilistic tsunami hazard analysis of the pacific coast of mexico: Case study based on the 1995 colima earthquake tsunami

This study develops a novel computational framework to carry out probabilistic tsunami hazard assessment for the Pacific coast of Mexico. The new approach enables the consideration of stochastic tsunami source scenarios having variable fault geometry and heterogeneous slip that are constrained by an extensive database of rupture models for historical earthquakes around the world. The assessment focuses upon the 1995 Jalisco–Colima Earthquake Tsunami from a retrospective viewpoint. Numerous source scenarios of large subduction earthquakes are generated to assess the sensitivity and variability of tsunami inundation characteristics of the target region. Analyses of nine slip models along the Mexican Pacific coast are performed, and statistical characteristics of slips (e.g., coherent structures of slip spectra) are estimated. The source variability allows exploring a wide range of tsunami scenarios for a moment magnitude ( M w ) 8 subduction earthquake in the Mexican Pacific region to conduct thorough sensitivity analyses and to quantify the tsunami height variability. The numerical results indicate a strong sensitivity of maximum tsunami height to major slip locations in the source and indicate major uncertainty at the first peak of tsunami waves.

An accurate assessment of tsunami hazards and quantification of uncertainty associated with the assessment are essential to mitigate and to control disaster risk exposures from a tsunami risk management point of view.Research on probabilistic tsunami hazard analysis/assessment (PTHA) has been improved after two mega seismic events, the 2004 Indian Ocean Tsunami (e.g., McCloskey et al., 2008) and the 2011 Tohoku Earthquake Tsunami (e.g., Mori et al., 2011).PTHA is a viable approach to evaluate the uncertainty of tsunami sources and related hazard modeling.One of the major challenges for tsunami impact assessment is to predict the earthquake source characteristics of future tsunamigenic events (e.g., location and geometric slip distribution), and then to quantify the uncertainty associated with the variability in earthquake rupture (tsunami generation) and tsunami inundation processes (e.g., Burbidge et al., 2008).In particular, tsunami generation is influenced by the complex and non-linear interaction of earthquake generation properties, while offshore tsunami propagation, affected by changes in sea bathymetry, is generally considered as a solved problem (Geist, 2002;McCloskey et al., 2008;Goda et al., 2014).
There are many scientific studies related to PTHA, which have been conducted worldwide.The earliest study that considered the probabilistic nature of tsunami hazard was by Rikitake and Aida (1988).Although they did not consider a full PTHA methodology, they used historical records and a typical earthquake fault model to estimate the probability of tsunami height exceeding a certain level at the shoreline.After the 2004 Indian Ocean event, there were several PTHA studies for other areas worldwide (e.g., Geist and Parsons, 2006;Annaka et al., 2007;Thio et al., 2007) and it accelerated after the 2011 Tohoku Earthquake Tsunami (e.g., Goda et al., 2015;Park and Cox, 2016).There are three major approaches for PTHA, which are commonly used [see a review article by Mori et al. (2017)].The first approach is to use a combination of many source scenarios based on expert opinion (e.g., González et al., 2009).The second approach is to use a logic-tree based on a combination of slip scenarios and geometric slip parameters with a weight function (e.g., Sánchez and Farreras, 1987;Annaka et al., 2007;Horspool et al., 2014;Fukutani et al., 2015;Lorito et al., 2015;Park and Cox, 2016).The third approach is to generate synthetic slip distributions, which are constrained by seismological theories and models, by slip wavenumber spectra assuming a random phase approximation (e.g., Geist and Oglesby, 2014;Goda et al., 2014Goda et al., , 2015;;Davies et al., 2015).The latter two methods are widely used for PTHA.Moreover, PTHA for landslidetriggered tsunamis is challenging and involves significant uncertainty in the tsunami generation process (e.g., Geist and Lynett, 2014).
The Pacific States of Mexico, i.e., Jalisco, Michoacan, Guerrero, and Oaxaca, are positioned at the subduction interface between the Rivera-Cocos Plates and the North American Plate (Bird, 2003).The slip rate along the plate boundary is in the range between 50 and 70 mm/year (DeMets et al., 1994), and hence, the potential for hosting large earthquakes is high.Historically, many large subduction earthquakes have occurred, causing severe shaking along the coastline and inland areas (e.g., Mexico City), and tsunami damage in coastal areas.Particularly, for this region, one of the major earthquakes that have been studied in the literature is the 28 March 1787 Earthquake (Suárez and Albini, 2009).According to the historical records, tsunami wave heights at specific locations had exceeded 10 m, and extensive inundation had occurred along southern Pacific Mexican coast (Nunez-Cornu et al., 2008).The estimated magnitude for this event ranges from Mw 8.4 to Mw 8.6.More recently, many moderate-to-large earthquakes have occurred in the Mexican subduction zone (Ramírez-Herrera et al., 2012).The large (>Mw 8) earthquakes in this zone include the 1932 Jalisco Earthquake (Farreras and Sánchez, 1991;Ramírez-Herrera et al., 2014), the 1985 Michoacan Earthquake (Mendoza, 1993), and the 1995 Jalisco-Colima Earthquake [denoted by the 1995 Colima Earthquake hereafter; Mendoza and Hartzell (1999)].
The offshore areas near Guerrero have not ruptured since 1911, which are prone to cause large earthquakes.This segment of the Cocos-North American Plate interface is referred to as the Guerrero seismic gap (Kostoglodov et al., 2003).The recent geophysical investigations by Perez-Campos et al. (2008) and Pacheco and Singh (2010) indicate that the Cocos Plate is subducting beneath the North American Plate with a dip angle of about 15°, and the slab reaches a depth of about 25 km at the distance of about 65 km from the Trench.In this region, both interface events and inslab events may be generated.Gradually, the subduction interface becomes flatter and becomes almost horizontal at the distance of about 120 km from the Trench (depth of about 40 km).The low coupling of the interface at the distance range beyond 100 km from the Trench results in low seismicity of inslab earthquakes.At the farther distance of about 300 km from the Trench, the subduction interface falls off sharply into the mantle.In the Guerrero seismic gap, some of the accumulated strain along the plate interface is also released as episodic slow slip events (Pacheco and Singh, 2010).Recently, one example of slow slip event sequences was triggered by the 2010 Maule, Chile Earthquake (Zigone et al., 2012).
Regarding future Earthquake-Tsunami events in the Guerrero region, major concerns are that (i) the segment has a potential to host large mega-thrust subduction events and (ii) no events greater than Mw 8 have not occurred in the segment in the recent history.In case of large subduction events, touristic places, such as Acapulco, will be devastated by earthquake and tsunami.Although Geist and Parsons (2006) showed PTHA results based on the logic-tree approach targeted Acapulco region, it is difficult to setup plausible scenarios in the Guerrero region due to lack of scientific data of historical events.A recent study by Perez-Campos et al. (2013) considered a scenario magnitude of Mw 8.2 for the Guerrero seismic gap, having the fault length and fault width of 210 and 90 km, respectively, and average slip of 4.0 m.Even when the magnitude of the scenario event can be defined for disaster risk mitigation purposes, geometry as well as slip distribution of the earthquake rupture may vary significantly.These rupture characteristics have significant influence on the earthquake ground motions as well as tsunamis.In addition, Jaimes et al. (2016) developed tsunami hazard maps in Mexico using several earthquake source scenarios and highlighted that it is essential to incorporate multiple earthquake source scenarios to produce reliable tsunami hazard maps.Therefore, it is important to take into account various earthquake rupture scenarios as part of disaster risk reduction strategy.Furthermore, there are many PTHA studies but the comparisons with historical records are limited.For benchmarking tsunami hazard predictions with past experience and future development of PTHA research, it is important to compare PTHA results to historical observations.This study develops a stochastic source model for large tsunamigenic earthquakes in the Guerrero region.The stochastic source modeling approach is useful for generating synthetic realizations of realistic source models (Mai and Beroza, 2002).The generated rupture models can be implemented in Monte Carlo tsunami simulation to assess the uncertainty of tsunami propagation and shoaling characteristics (Goda et al., 2014), and such uncertainties can be further propagated in tsunami hazard and damage estimation to promote effective tsunami risk mitigation decisions (Goda and Song, 2016).Recently, new  scaling relationships of various earthquake source parameters, such as fault geometry and slip statistics, have been developed for tsunamigenic events (Goda et al., 2016).The new empirical scaling models are based on an extensive statistical analysis of numerous source inversion models of the past earthquakes obtained from the SRCMOD (Mai and Thingbaijam, 2014; see http://equake-rc.info/SRCMOD/), which is an online database of finite-fault rupture models of past earthquakes.Using these equations, source parameters can be predicted by taking into account uncertainty and dependency of other parameters (note: source parameters are physically inter-related and thus prediction errors of these parameters are correlated).In short, stochastic tsunami simulation is valuable for assessing the regional tsunami impact due to future large earthquakes in the Guerrero region.
The study is organized as follows.First, a summary of the finite-fault rupture models for Mexican subduction earthquakes is presented; nine source models are obtained from the SRCMOD database.Based on the characteristics of the source models, a generic fault model for stochastic source modeling in the Guerrero region is developed.In synthesizing earthquake source models, new scaling relationships are implemented in the Monte Carlo tsunami simulation.Then, a numerical procedure of the stochastic tsunami simulation is described for the Guerrero region.Finally, an application of the stochastic tsunami simulation for the 1995 Colima Earthquake is presented, which facilitates the retrospective investigation and comparison with observed tsunamis during the historical event.

analYsis OF hisTOrical sOUrces analysis of Fault Models
Numerous source models have been developed for historical Mexican subduction earthquakes.In the well-organized source database, such as SRCMOD (Mai and Thingbaijam, 2014), nine finite-fault models are available for the Mexican subduction zone.The locations of the nine finite-fault models are shown in Figure 1 (see also Table 1).The strike angles of the finite-fault models range from 283° to 309°, which are consistent with the boundary between the Cocos Plate and the North American Plate (Bird, 2003).The dip angles of the source models are typically in the range of 12-14°, except for Model 7.
Using these finite-fault model data, Goda et al. (2016)  a Model 6 is not analyzed because this model was judged as unsuitable for spectral analysis based on the criteria set in Goda et al. (2016).The gray shade in the table indicates target of this study.
scaling relationships for source Parameters for guerrero region The stochastic source models are synthesized based on a set of new scaling relationships of source parameters developed by Goda et al. (2016).Equations 1-6 are the relationships for W, L, Da, Dm, Az, and Ax, respectively, and are given as a function of Mw.
The error terms of the source parameters are correlated; The equations for the Box-Cox parameter and the Hurst number are independent of Mw (Goda et al., 2016).The Box-Cox parameter is modeled as a normal random variable with mean equal to 0.312 and SD equal to 0.278.On the other hand, the Hurst number is modeled as a random variable that takes a deterministic value of 0.99 with probability of 0.43 and a sampled (random) value from the normal distribution with mean equal to 0.714 and SD equal to 0.172 with probability of 0.57.See Goda et al. (2016) for further details of the probabilistic models of the source parameters.The prediction errors of the Box-Cox parameter/Hurst number are considered to be uncorrelated with those of other source parameters.
The consistency of the estimated source parameters for the eight finite-fault models of the Mexican subduction earthquakes (Table 1) with the prediction models of the source parameters mentioned above is examined, and the results are shown in Figure 3.In Figure 3, the results for the 1995 Colima Earthquake by Mendoza and Hartzell (1999) are represented by a different symbol.The stochastic tsunami simulation of the 1995 Colima Earthquake is carried out, in comparison with the field observations, in the latter part of this study.The results shown in Figure 3 indicate that the estimated source parameters agree with the scaling relationships; for most cases, the estimated parameters fall within the 16th-84th percentile confidence interval of the prediction equations, Eqs 1-6.Therefore, the use of the developed scaling relationships by Goda et al. (2016) for Mexican subduction events in the Guerrero region can be justified by the historical slip models along the Pacific Mexican subduction zone.

OUTline OF nUMerical MODel numerical Model for Tsunami simulation
A series of numerical simulations for tsunami propagation from the source to coastline is performed by non-linear shallow water equations by Goto et al. (1997).The governing TaBle 2 | Linear correlation coefficients of regression residuals of the scaling relationships for the six earthquake source parameters.individual grid systems with different resolutions.The coarsest 810 m region covers the entire Guerrero region, and the deformation due to fault rupture is computed at this resolution using Okada (1985) and Tanioka and Satake (1996)  equations are evaluated using a leap-frog staggered-grid finite difference scheme.The nesting grid systems that are implemented considering the size of continental shelf for the Guerrero region and nearshore bathymetry have three levels as shown in Figure 2B.
The grid discretization at the coarsest level is 810 m, while the finest level is 90 m; a factor of 3 is considered to connect An accurate inundation calculation is excluded in this study, although inland inundation is simulated during the computation.Therefore, the main results presented in this study focus on the tsunami heights along the coast.Note that the tsunami wave height that is discussed in this study is the height of water flow above mean sea level.It is important to emphasize that GEBCO2014 has bias in very shallow waters and the data also need to match with water depths along the shoreline.Thus, GDEM2 data (2011) are used for onshore topography.The integration of bathymetry data and elevation data are not trivial because the spatial resolutions of these data are very different.For the cases of GEBCO2014 and GDEM2, the resolutions differ by a factor of 30 (i.e., 900 versus 30 m).The effects of the interpolation are expected to be significant at shallow depths near the shoreline.In developing the "depth" data for tsunami simulation (i.e., combined elevation data for a given region), first, three datasets, namely GEBCO2014, GDEM2, and SWBD, are combined without interpolation.The points in the "combined" data are spaced neither regularly nor uniformly.Within the onshore areas, the corresponding GEBCO2014 data between 0 and 200 m in elevation are replaced by the counterparts of GDEM2 data.In addition, the SWBD shoreline data are overlaid as zero elevation data points.Once this composite dataset is developed, linear interpolation is carried out.In future studies, new bathymetry data (probably compiled from local sources) should be incorporated to improve the reliability and accuracy of the tsunami simulation, especially in the very shallow water environment.

stochastic Tsunami simulation
Stochastic tsunami simulation can be conducted by generating multiple stochastic source models for a given earthquake scenario and by performing tsunami forward modeling repeatedly.A computational flowchart of stochastic tsunami simulation is shown in Figure 5.The detail of stochastic tsunami modeling is available in Goda et al. (2014Goda et al. ( , 2016) but it will be explained briefly in the following.
The first step of the method is to define a suitable tsunami source zone model.For the Guerrero region, the model shown in Figure 2A is adopted.The scenario magnitude should be selected according to the objective of the analysis.Within the fault plane, the so-called asperity zone is setup, together with the required slip concentration range.Essentially, the asperity zone works as crude constraints of the generated source model regarding the slip concentration within the fault plane.It requires that a certain amount of slip must be concentrated within the target region.One example is that more than 50% of the total slip should be concentrated in the shallow part of the fault plane (e.g., shallower than 20 km).The asperity zone parameters should reflect the seismological knowledge of earthquake rupture in the target region.
Second, the macro earthquake source parameters, such as W, L, Da, Dm, λ, Az, Ax, and H, are generated using the scaling relationships.Uncertainty as well as correlation associated with the regression models should be taken into account in sampling the values of the source parameters (Eqs 1-6 and Table 2).In the simulation, random variables for these residuals can be sampled from the multivariate normal distribution.In addition, at this stage, consistency among the simulated values of W, L, and Da can be tested by comparing the target seismic moment (as specified by the given scenario magnitude) and the simulated seismic moment (Mo = μWLDa, where μ is the rock rigidity).An inadequate combination of W, L, and Da values is resampled until the scenario magnitude is satisfied.
Third, using the generated spatial slip distribution parameters, a random slip field is generated using a Fourier integral method (Pardo-Iguzquiza and Chica-Olmo, 1993).To achieve slip distribution with realistic positive skewness, the synthesized slip distribution is converted via Box-Cox transformation (Goda et al., 2014).The transformed slip distribution is then adjusted to achieve the target mean slip Da and to avoid very large slip values exceeding the target maximum slip Dm.Subsequently, the position of the synthesized fault plane is determined randomly within the source region.To ensure that the synthesized slip distribution is realistic with respect to the seismotectonic characteristics of the region, two criteria/constraints are implemented to determine the final acceptance of the generated source model.The first constraint requires that the asperity area ratio of the candidate slip distribution falls between 0.2 and 0.3.The second constraint requires that the simulated earthquake slip is more concentrated in the designated asperity region.Multiple slip distributions are simulated repeatedly until an acceptable source model, which has all desirable characteristics, is obtained.
Fourth, for a given acceptable source model, the initial water surface elevation (i.e., initial boundary conditions for tsunami simulation) is evaluated based on formula by Okada (1985) and Tanioka and Satake (1996).Tsunami wave propagation is evaluated by solving non-linear shallow water equations (Goto et al., 1997).The run-up of tsunami is considered by the model but it is incomplete due to coarse grid size in this study.Finally, the above simulation procedure is repeated until a sufficient number of acceptable source models are generated and their tsunami inundation heights/depths at locations of interest are evaluated.The results from the Monte Carlo tsunami simulation are useful for evaluating variability of tsunami simulation results at different locations and for developing stochastic tsunami hazard maps along the coast.synthetic tsunami characteristics.To demonstrate the stochastic tsunami simulation model for the Guerrero region, the 1995 Colima Earthquake is considered, which is one of the most major tsunami events in the northern part of the Guerrero region.More specifically, the selection of the 1995 Colima Earthquake is relevant because the size of the earthquake (Mw 8.0) is sufficiently large to cause tsunami waves and post-event tsunami survey data (e.g., Borerro et al., 1997;Trejo-Gómez et al., 2015) as well as an inverted slip model (Mendoza and Hartzell, 1999) are available for this event.Our aim in setting up the case study is to compare the results of stochastic tsunami simulations with the past survey data of the 1995 Colima Earthquake.
The moment magnitude of the event was Mw 8.0 and occurred near the junction of the Cocos-Rivera-North American Plates (Model 5 in Figure 1).Borerro et al. (1997) reported that the observed run-up heights along the Colima coast line (longitudes between 104 and 105°W) were mostly 2-3 m and was ranged up to 5 m (note: at one particular point, north end of Santiago Bay, the observed run-up height reached 10.9 m; this high run-up was likely to be caused by the very local topographical effect).
To setup the stochastic tsunami simulations for the 1995 Colima Earthquake, a separate source zone model is defined within the whole source zone model for the Guerrero region.The 1995 Colima Earthquake source zone model restricts ranges of slip along the subduction zone shorter than the Guerrero region and is based on the source model by Mendoza and Hartzell (1999), i.e., Model 5 in Figure 1 and Table 1.A zoom-up of the Mendoza-Hartzell fault plane model ( 1999) is shown in Figure 6A.The maximum slip is 4.78 m, and it is located around 19°N in latitude and 106-105.5°W in longitude.To incorporate the variations of the fault plane size, the source zone model for the 1995 Colima Earthquake has the width equal to 150 km and the length equal to 310 km, which is larger than the original fault plane by Mendoza and Hartzell (1999) (note: the top edge of the fault plane model is identical).It is also noteworthy that the Mendoza-Hartzell model agrees with the scaling relationships developed by Goda et al. (2016) shown as the square marks in Figure 3.

synthetic Tsunami simulation
To carry out the stochastic tsunami simulation, two sets of 100 stochastic source models are generated for the Mw 8.0 scenario.The first set takes into account uncertainty of the scaling relationships of the source models, whereas the second set does not.The practical consequences of not incorporating the uncertainty of the scaling relationships are that the generated source models have identical dimensions and slip heterogeneity parameters, while the slip distribution and location of the source model are varied.On the other hand, when the uncertainty of the scaling relationships is considered, both dimensions and slip heterogeneity are varied.Note that such uncertainty sometimes is ignored, and scaling relationships are used as deterministic equations relating two source parameters.In addition to the two sets of 100 stochastic source models, a hindcast simulation for the 1995 Colima Earthquake Tsunami is performed by using the original Mendoza-Hartzell fault plane model (Figure 6A) as a benchmark.Figure 6B shows examples for four realizations of the stochastic source models for the cases where the uncertainty of the scaling relationships is accounted for.The synthetic slips in Figure 6B demonstrate that realistic locations of maximum slip along the Trench and coherent structure of sub-slips around the maximum slip given by spectral decomposition.Note that each realization has the similar tail of slip spectra, which are controlled by the spectral decomposition and Box-Cox parameter.The generation of arbitrary number of synthetic slips is one of advantages of using the random phase approach in PTHA, which is different from the logic-tree approach.Therefore, it is possible to discuss probabilistic tsunami heights along the coast.Although the examples of stochastic source models without the uncertainty of the scaling relationships are not shown due to limitation of space, the synthetic source models considering the uncertainty of the scaling relationships shown in Figure 6B change fault size, i.e., L and W, and the other slip characteristics due to the uncertainty terms ε in Eqs 1-6.Consequently, the rupture aspect ratios (i.e., L/W) become variable even when the same  magnitude is considered.To discuss the tsunami heights along the coast probabilistically and to investigate the uncertainty effects on earthquake source modeling, the tsunami simulation results are presented in two sections, i.e., sensitivity of tsunami simulated heights and effects of accounting for parameter uncertainty in earthquake source generation.

Sensitivity of Tsunami Simulated Heights
Figure 7 shows the maximum wave height of the 1995 Colima Earthquake Tsunami simulated by the Mendoza-Hartzell model.The large tsunami heights are located in the northern part of the domain due to the large slip concentration as shown in Figure 7.There is less significant tsunami amplification inside bays; however, the amplification of tsunami heights along the northern coast facing the Pacific due to shallow water shoaling effects is remarkable.The tsunami height in the southern part of the domain is less than 1 m, noting that this should be regarded one of realizations from ensemble events.The mean and SD of maximum tsunami heights based on the 100 stochastic models considering the uncertainty of scaling relationships are shown in Figure 8.Besides the southern end of the computational domain, the mean maximum tsunami height by the stochastic models are more uniformly distributed along the coast compared with other sites, although the wave deformation from the offshore to onshore are different especially in the middle of the domain.The SD of maximum tsunami heights along the coast is large outside of bays facing the Pacific but is not significant inside.
The comparison of the mean maximum tsunami heights from the stochastic tsunami simulations and the maximum tsunami height of the hindcast tsunami simulation is conducted by calculating the ratio between them.The results are shown in Figure 9.The ratio of the maximum tsunami heights, shown in Figure 9, is mostly larger than 1.Thus, the predicted tsunami hazards based on the inverted source model developed using actual observed geophysical data for the 1995 Colima Earthquake Tsunami are less than the averaged realization of the stochastic tsunami simulation.There are three major factors to amplify tsunami heights at particular coastal locations.The first is source characteristics (e.g., location of maximum slip), the second is wave shoaling and focusing due to large-scale bathymetry features, and the third is convergence of energy due to the shape of a bay.The second effects, i.e., wave shoaling, are proportional to h −1/4 where h is water depth, if wave nonlinearity is negligible (i.e., Green's law). Figure 10 shows the spatial distribution of h −1/4 of GEBCO2014 data around the target area shown in Figure 7.The value of h −1/4 constantly increases from offshore to onshore lower than 19.3°N in latitude.Therefore, the amplification of tsunami heights shown in Figure 9 is not influenced by offshore bathymetry.It can be concluded that the source characteristics are the main cause of the observed spatial inhomogeneity of mean maximum tsunami height rather than the tsunami propagation processes.
The maximum tsunami heights along the coastal line (elevation is between -1.0 and 1.0 m) are then extracted for every 450 m (i.e., five grids) in the longitudinal direction to examine the spatial variation of the tsunami heights along the coast; the results can be compared with the tsunami survey results.The locations of the extracted points are shown in Figure 11.In total, there are 284 sites along the coast.In addition, tsunami waveforms (i.e., temporal profile) at three recording locations along the Colima coast are also focused on to investigate the temporal characteristics of the tsunami profiles for different source models.In relation to the surveyed locations by Borerro et al. (1997), the recording location 2 is in Tenacatita Bay, whereas the recording location 3 is in Manzanillo Bay.The depths at the recording locations 1-3 are 1.3, 10.1, and 36.2 m, respectively.
First, the maximum tsunami heights along the coast line from north to south are shown in Figure 12.The mean and the upper and lower 16th percentiles of the simulated tsunami heights with and without the uncertainty of the scaling relationships (red solid and dashed line) are shown in the figure.These results are also compared with the hindcast simulation of the 1995 Colima Earthquake Tsunami (blue solid line) and observed run-up heights (dots) by Trejo-Gómez et al. (2015).The past tsunami height profile of the 1995 Colima Earthquake shows that the wave heights for sites 1-100 are higher than others, which are consistent with the source model by Mendoza and Hartzell (1999), having large asperities in the north-western segment of its fault plane.The maximum tsunami profiles for the stochastic source models vary significantly along the coast; for instance, the median curve varies between 2.0 and 5.5 m, while the 84th percentile curve varies between 3.0 and 7.5 m.The spatial distributions of mean maximum tsunami heights from sites 0 to 70 are generally similar to the historical run as well as the measured tsunami run-up heights by Trejo-Gómez et al. (2015).On the other hand, the maximum tsunami height profile for the stochastic source models from the sites 100 and more (toward south-east) differs from the historical run and measured tsunami run-up heights, noting that the locations shown in Figure 2 of Borerro et al. (1997) approximately correspond to sites 100-200.It emphasizes that the Mendoza-Hartzell source model is not based on the tsunami data; therefore, disagreement between the observed and the simulated tsunami results is not unexpected.

Effects of Accounting for Parameter Uncertainty in Earthquake Source Generation
Regarding the effects of the uncertainty of the scaling relationships (i.e., Figure 12A versus Figure 12B), it can be seen that the variability of the simulation results is significantly increased, especially for the maximum tsunami height profiles (although general trends of the results are similar for both cases).The mean and upper 16th percentile values of the simulated results with the uncertainty of the scaling relationships are larger than the results without the uncertainty of the scaling relationships.Especially, when the uncertainty of the scaling relationships is taken into account, the upper 16th percentile and upper bound of the maximum tsunami height, shown as background gray lines in Figure 12, become quite large and are related to the increase of the maximum slip due to the uncertainty by Eq. ( 4).Thereby, the estimation of scaling relationship for the maximum slip and its uncertainty by analyzing inversion slip models are important for tsunami hazard assessment.The impact of stochastic modeling on time series of tsunami profile at the recording locations 1-3 (from north to south) are illustrated in Figure 13.The tsunami waveforms at the recording locations 1-3 show that the wave amplitudes are generally higher at the recording location 1 than those at the recording stations 2 and 3.The simulated tsunami waveforms for the recording locations 1-3 also exhibit large variations in the temporal tsunami profiles, indicating that tsunamis having amplitudes up to 5 m may be expected at offshore locations for the Mw 8.0 earthquake scenario.The largest variability due to different stochastic source models is found at the first peak of tsunami waves; such large variability is also observed for a few subsequent waves after the first one.However, the signs of the first wave do not change over the stochastic simulations, and the first waves always begin with positive change.Moreover, edge waves following continental shelf are mainly affected by large-scale nearshore bathymetry, which is independent of slip source modeling.Consequently, there is less influence due to the stochastic source models on the later parts of the simulated waves.

cOnclUsiOn
PTHA for the Guerrero region, Mexico has been carried using the novel stochastic tsunami simulation method.The method takes into account uncertainty of the key source parameters and randomness of slip heterogeneity over the fault plane and, hence, is capable of quantifying the tsunami hazard probabilistically.Such a methodology has not been implemented in the previous PTHA studies for Mexico.The scaling relationships used in the stochastic earthquake source generation have been developed based on extensive statistical analyses of the source models parameters estimated from the SRCMOD database.The bathymetry and elevation data for the region were compiled based on the GEBCO2014 and GDEM2 to develop the nesting grid systems that are suitable for regional tsunami simulation studies.Finally, the developed stochastic tsunami simulation method was applied to the 1995 Colima Earthquake scenario.The results indicated that the effects of the source model characteristics on the simulation results are important.It was also found that the tsunami simulation results using the stochastic source models exhibit significant variability of tsunami profiles, while the results overall agree with the tsunami run-up survey results for the 1995 Colima event.
The extension of the source zone model to the Guerrero region by varying earthquake scenario magnitudes will be the focus of our future study.Such investigations have been carried out for Japan using the similar stochastic tsunami simulation method (Goda et al., 2017).It is also important to simulate tsunami inundation and run-up and to assess tsunami damage to structures based on the stochastic tsunami simulations along the Pacific Mexican coast.

FigUre 1 |
FigUre 1 | Finite-fault models, along the Mexican Pacific coastline, available in the SRCMOD database (see details of parameters in Table1).
evaluated the macro source parameters: fault length (L), fault width (W), mean slip (Da), maximum slip (Dm), Box-Cox parameter (λ), correlation length along strike direction (Az), correlation length along dip direction (Ax), and Hurst number (H) as a function of moment magnitude.The fault width and length, together with strike and dip, define the geometry of the fault plane.The mean slip, maximum slip, and Box-Cox parameter characterize the probability distribution of the slip values.The correlation lengths and Hurst number are used to model the spatial heterogeneity of the slip values.The obtained values of the source parameters for the Mexican subduction earthquakes are listed in Table1, except for Model 6.Because the majority of the slip values of Model 6 are 0 (Figure1), Model 6 was regarded as unsuitable for spectral analysis of the source model and thus excluded from further investigations in this study.Based on the geometry of the Mexican finite-fault models, a generic fault model for the Guerrero region is defined for the synthetic source generation (Figure2).It covers the offshore region of the Pacific Mexican coast.The length and width of the Guerrero source zone are 930 and 170 km, respectively.The top edge of the fault plane is positioned at a depth of 3 km.The fault plane has a constant strike of 293° and a constant dip of 13°.For stochastic source modeling and Monte Carlo tsunami simulation, the Guerrero source zone is discretized into 10/10 km sub-faults.In stochastic source modeling, slip values that are consistent with the considered spatial slip distribution characteristics are generated.TaBle 1 | Summary of the finite-fault source parameters for the Mexican subduction earthquakes.

FigUre 2 |
FigUre 2 | Tsunami source zone model and tsunami propagation domain for the Guerrero region.(a) Tsunami source zone model.(B) Tsunami computational domain (810 m-270 m-90 m).
equations for the region shown in Figure 2A.There are two 270 m regions shown in Figure 2B; within the 270 m −1 and 270 m −2 regions, five 90 m regions are defined, respectively.The nested domains of 270 and 90 m resolutions are shown in Figure 2B.It is noted that at the same resolution, regions overlap with the nearby (same-resolution) regions.This is to ensure that the solutions of tsunami waves, especially edge wave, in simulation are propagated across different regions properly.The duration of each simulation is set to 2 h and the time step for the simulation is 0.25 s to satisfy the Courant-Friedrichs-Lewys condition, which is a necessary condition of time and space discretization for convergence when partial differential equations are numerically solved by the finite difference method.No tidal variation is taken into account.

FigUre 3 |
FigUre 3 | Comparison of the estimated source parameters for the eight finite-fault models of the Mexican subduction earthquakes with the corresponding scaling relationships.
resUlTs anD DiscUssiOn: case sTUDY OF The 1995 cOliMa earThQUaKe synthetic slips A brief discussion of generated synthetic/stochastic tsunami sources for the Guerrero region is given before analyzing

FigUre
FigUre 6 | Continued

FigUre 6 |
FigUre 6 | Stochastic source models for the 1995 Colima Earthquake scenario by considering prediction errors of the scaling relationships.(a) Mendoza-Hartzell model for the 1995 Colima Earthquake.(B) Synthetic source models.

FigUre 7 |
FigUre 7 | Maximum surface elevation, in meters, of the 1995 Colima Earthquake Tsunami based on the Mendoza-Hartzell model.

FigUre 8 |
FigUre 8 | Mean and SD of maximum surface elevation, in meters, by the stochastic models considering prediction errors of the scaling relationships.(a) Mean.(B) SD.

FigUre 9 |
FigUre 9 | Ratio of average maximum surface elevation by the stochastic model considering prediction errors of the scaling relationships to maximum surface elevation calculated based on the Mendoza-Hartzell model.

FigUre 11 |
FigUre 11 | Site locations (1-284) and wave recording locations (1-3) for the analysis of time series of tsunami surface elevation.

FigUre 12 |
FigUre 12 | Maximum tsunami wave height profile along the coastal line for the stochastic models with and without considering prediction errors of the scaling relationships [red solid line: mean, red dashed line: upper and lower 16%, gray line: individual stochastic run, blue solid line: historical run, circle: observed run-up height by Trejo-Gómez et al. (2015)].(a) Considering prediction errors of the scaling relationships.(B) Without considering prediction errors of the scaling relationships.

Table 2
lists the linear correlation coefficients of the prediction errors for Eqs 1-6.