Rapid prediction of alongshore run-up distribution from near-field tsunamis

Rapid prediction of the spatial distribution of the run-up from near-field tsunamis is critically important for tsunami hazard characterization. Even though significant advances have been made over the last decade, physics-based numerical models are still computationally intensive. Here, we present a response surface methodology (RSM)-based model called the tsunami run-up response function (TRRF). Derived from a discrete set of tsunami simulations, TRRF can produce a rapid prediction of a near-field tsunami run-up distribution that takes into account the influence of variable local topographic and bathymetric characteristics in a given region. This new method reduces the number of simulations required to build an RSM model by separately modeling the leading order contribution and the residual part of the tsunami run-up distribution. Using the northern region of Puerto Rico as a case study, we investigated the performance (accuracy, computational time) of the TRRF. The results reveal that the TRRF achieves reliable prediction while reducing the prediction time by six orders of magnitude (computational time: <1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$< 1$$\end{document} s per earthquake).


Introduction
Tsunamis are some of the most destructive and costly natural hazards for coastal areas around the world. The 2004 Indian Ocean tsunami and the 2011 Tohoku tsunami are prime examples of how tsunamis can cause extensive damage to coastal communities, especially in near-field areas (Titov et al., 2005;Wei et al., 2013). To mitigate damage and build resilient coastal communities, it is critically important to develop rapid prediction capacities for a near-field tsunami run-up distribution along the coastlines. Physics-based numerical simulation is currently the most accurate method for predicting a tsunami run-up distribution. Though significant advances have been made over the last decade (LeVeque et al., 2011;Lin et al., 2015;Popinet, 2015;Shi et al., 2012), these physics-based numerical models still remain time consuming. For example, robust probabilistic tsunami hazard assessment (PTHA) requires tsunami run-up estimates for a large number of scenarios to allow for accurate quantification of the hazard and related uncertainty (Mori et al., 2018). However, due to the computational burden associated with physics-based numerical simulation, a logic-tree approach is typically employed: it limits the number of scenarios based on historical earthquake characteristics (e.g., magnitude, recurrence interval) used to evaluate uncertainty in tsunami hazard (Annaka et al., 2007;Park and Cox, 2016;Park et al., 2018). The issue with the logic-tree approach is that it relies on expert judgment, which is difficult to justify. Also, the computational burden associated with physics-based numerical simulation-especially for near-field tsunami forecasting-is a major obstacle. For this reason, pre-computed simulation databases are widely used. These databases can provide fast prediction by selecting the best-matched simulation or by interpolating between simulations immediately after the source mechanism is known (Kamigaichi, 2011;Mulia et al., 2018;Setiyono et al., 2017). A problem with the database approach is that it can have substantial errors in real-world scenarios that do not exist in the selected databases.
The response surface methodology (RSM) is an effective statistical-based approach for establishing a relationship between a set of input variables and the output of a system (Box and Wilson, 1951;Myers et al., 2016). Once the RSM model is built, output can be rapidly estimated across the continuum of input spaces. However, because high-dimensional input requires a large number of simulations-which is prohibitively expensive-the RSM has not been used to predict a tsunami run-up distribution. For example, a tsunamigenicearthquake (the input in an RSM model) is usually represented by nine fault parameters (Fig. 1). A full factorial design is one of the most widely employed designs of experiments (DoE) used to measure the response of every possible combination of independent variables. If we design the synthetic tsunami scenarios using a three-level full factorial approach with nine fault parameters, 19,683 (= 3 9 ) simulations are required. Moreover, if the input/output relationship shows large nonlinearity, a higher level of DoE may be needed, which would necessitate exponentially more simulations. Here, we present a new methodology to rapidly predict the near-field tsunami run-up distribution: the tsunami run-up response function (TRRF). It is based on RSM but requires only 729 (= 3 6 ) simulations through reducing input dimensionality. Input dimensionality is reduced through a decomposition of the leading order tsunami run-up contribution and the residual part of the run-up distribution. We demonstrated the TRRF approach in northern Puerto Rico, where a significant tsunami generated by an earthquake along the Puerto Rico Trench could devastate coastal communities on the northern shore (Grilli et al., 2010;López-Venegas et al., 2015;Reid and Taber, 1919) (Fig. 2).

Tsunami run-up response function (TRRF)
The main concept of TRRF is to decompose the tsunami run-up distribution R(x) into source run-up S(x) and topographic run-up T (x) (Fig. 3): where the x-axis is parallel to the coastline. The source run-up S(x) is a leading order contribution that can be represented by Okal and Synolakis (2004)'s empirical formula (hereafter OS formula): where the coefficient a is related to the width of the source run-up, b is the maximum source run-up, and c is the distance from the x-axis origin to the location of the maximum source run-up.   (Mw ≥ 4.5, (USGS, 2017)). The filled black circles and dashed black lines represent the epicenters of NOAA's pre-defined unit sources and fault orientation, respectively (Gica, 2008). The dashed red line represents the contour line where the water depth is 8 km. The blue dashed square represents the region where the National Geophysical Data Center (NGDC)'s 3-second topographic grid (NGDC, 2005) is used for numerical simulation.
The topographic run-up T (x) is the residual run-up remaining after subtracting S(x) from R(x). It represents the local (de)amplification of the incoming tsunami wave and the resulting run-up arising from topographic variation. The T (x) can be normalized as follows, hereafter called normalized topographic run-up N T (x).
The TRRF can predict a tsunami run-up distribution R p (x) by putting the source run-up S p (x) and the normalized topography run-up N T p (x) (where superscript p represents prediction) to the following equation: where N T p (x) is the 50th percentile (or median) of N T (x) among all the simulations used to build the TRRF. The source run-up S p (x) can be estimated by inputting the OS formula coefficients a, b, and c into Eq. 2. The RSM approach is applied to estimate the OS formula coefficients a and b from six fault parameters (hereafter RSM fault parameters): where f a and f b are the best-fitting curves to these coefficients; hereafter, these curves are called RSM functions. Since the RSM function inputs consist of six fault parameters, 729 (= 3 6 ) simulations are required to derive RSM functions following a three-level full factorial design (See Appendix A). Epicenter longitude LON can be considered separately because northern Puerto Rico's coastline runs parallel to the latitude lines. In this condition, the coefficients a and b are independent of epicenter longitude, and the coefficient c can be represented by a function of epicenter longitude: Step 1 Step 2 Step 3 Fig. 4 Computational flow of TRRF development.
where LON 0 is the longitude of the origin. The strike angle ST R and rake angle RAK are also not included in the RSM function inputs because the OS formula (Okal and Synolakis, 2004) is only applicable to an earthquake fault oriented in shore-parallel strike direction (= 90 • for northern Puerto Rico) with 90 • rake angle. Since this is not the only case that occurs in nature, we developed a method that can represent a fault where strike angle and/or rake angle are not 90 • as a series of hypothetical faults where both the strike angle and the rake angle are 90 • , hereafter called the angle projection (AP) method. Section 2.1 will describe the procedures of building a TRRF. Section 2.2 will explain how to apply the AP method to predict a tsunami run-up distribution for the case where strike angle and/or rake angle are not 90 • . Lastly, Section 2.3 will describe the procedures of predicting a tsunami run-up distribution once the TRRF is built. Fig. 4 shows the procedure of TRRF development. The first step is to simulate 729 tsunamigenic-earthquake scenarios using a physics-based numerical model. The second step is to extract the run-up and apply the OS formula (Eq. 2) to obtain the normalized topographic run-up N T p (x). The last step is to fit the earthquake fault parameters and the OS formula coefficients (a and b) North STR RAK East Fig. 5 Schematic sketch of step 1 of AP method. The yellow-filled rectangle is the original fault where ST R is strike angle and RAK is rake angle. The red rectangle represents the adjusted fault where θ is the adjusted strike angle and λ is the adjusted rake angle. The arrows represent the slip direction.

TRRF development
to the 2nd order polynomial model to obtain the RSM functions. Once the N T p (x) and the RSM functions are derived, this procedure does not have to be repeated to predict the tsunami run-up distribution.

Angle projection (AP) method
The AP method comprises three steps: adjustment of strike angle and rake angle, fault rotation, and decomposition of slip. Note that the units of LON and LAT in this section are in kilometers.

Adjustment of strike angle and rake angle
The direction of near-field tsunami propagation is related to the interaction between the strike angle ST R and the rake angle RAK. To consider the interaction between ST R and RAK, the first step involves adjusting the ST R and the RAK as follows (Fig. 5): where θ is the adjusted strike angle and λ is the adjusted rake angle. The α, β, and gamma are the site-specific coefficients (see Section 4).

Fault rotation
If the adjusted strike direction is parallel to the coastline (θ = 90 • ) and the adjusted rake angle λ is 90 • , the maximum source run-up will be located at a shore perpendicular from the longitude epicenter LON . On the other hand, if the adjusted strike direction is not parallel to the coastline (θ = 90 • ), the location of the maximum source run-up will be shifted to a direction perpendicular to the adjusted strike direction. To consider the location of the maximum source run-up depending on the adjusted strike angle, the second step involves rotating the adjusted fault (θ = 90 • ) until θ becomes 90 • (Fig.  6). The epicenter of the rotated fault (LAT p 1 , LON p 1 ) can be calculated as follows:

Decomposition of slip
If the adjusted rake angle λ is not 90 • , the run-up will be spread in the slip direction. To consider the spread of run-ups depending on the adjusted rake angle, the third step involves representing the rotated fault (λ = 90 • ) as a series of hypothetical faults having slips perpendicular to the coastline (λ = 90 • ) (Fig. 7). Since the tsunami energy is proportional to SLP , we assume that the source run-up will be spread proportionally to a component of SLP parallel to the coastline. Based on this assumption, while the LEN , W ID, DEP , and DIP are identical to the original fault, the epicenter (LAT , LON ) and SLP of the ith hypothetical fault (i = 1, 2...n) are defined as follows: where n is the total number of hypothetical faults (see Section 4). The SLP p 1 and SLP p n are the slips of the first and last hypothetical faults, respectively, defined as follows: The LAT p n and LON p n are the epicenter latitude and longitude of the last hypothetical fault and can be calculated based on a geometric setup (see two red lines in Fig. 7): Following this procedure, a fault where strike angle and/or rake angle are not 90 • can be converted into a series of hypothetical faults where both strike angle and rake angle are 90 • .

Run-up
Schematic sketch of step 3 of AP method. The green rectangle represents the first hypothetical fault where the epicenter is (LAT p 1 , LON p 1 ) and the slip is SLP p 1 . The blue rectangle represents the last hypothetical fault among a series of hypothetical faults where the epicenter is (LAT p n , LON p n ) and slip is SLP p n . Gray circles and arrows represent the epicenters and the slips of the hypothetical faults, respectively; these are linearly distributed between the first hypothetical fault and the last hypothetical fault. Two red lines are of the same length. One line is parallel to the slip direction spanning from the epicenter of the rotated fault to the point where it meets the coastline. The other line is the vertical distance from the epicenter of the last hypothetical fault to the coastline. Fig. 8 shows the procedure for how the TRRF predicts a tsunami run-up distribution once the TRRF is built. The first step is to convert the earthquake fault into a series of hypothetical faults using the AP method. The second step is to estimate the OS formula coefficients a p i , b p i (i = 1, 2...n) of hypothetical faults using the RSM functions. The third step is to estimate the OS formula coefficient c p i (i = 1, 2...n) by inputting the epicenter longitudes into Eq. 7. Step 3

TRRF application for prediction
Step 2 Step 4 Step 5 Fig. 8 Computational flow of TRRF application for prediction.
The fourth step is to estimate the final source run-up S p (x) by inputting the OS formula coefficients into Eq. 2 and taking the maximum values of the estimated source run-ups for all hypothetical faults. Finally, the tsunami runup distribution R p (x) can be estimated by inputting the source run-up S p (x) and the normalized topographic run-up N T p (x) to Eq. 4.

Numerical simulation
The 729 tsunamigenic-earthquake scenarios were simulated based on the numerical model Basilisk, which solves the Green-Naghdi equations and employs both Adaptive Mesh Refinement (AMR) and parallelization to facilitate efficient computation.The Basilisk model has not only been validated with several benchmark problems but also been applied to simulate the 2011 Tohoku tsunami (Popinet, 2015). The 729 scenarios were designed as shown in Table  1. The range of the epicenter latitude LAT was determined based on National Oceanic and Atmospheric Administration (NOAA)'s pre-defined unit sources and historical earthquake records in northern Puerto Rico (Fig. 2). The range of the fault length LEN , fault width W ID, and slip SLP was set based on the assumption that the moment magnitude (M w ) should be larger than 7.0 for a tsunami to occur. We used the empirical regression of Hanks and Kanamori (1979) and fundamental equation of Aki (1966) to calculate the moment mag-nitude: where M 0 is a seismic moment (N m), µ is rigidity modulus of the Earth's crust (N m −2 ), and the units of LEN , W ID, and SLP are in meters. We assumed that the rigidity modulus µ is 4.2 × 10 10 N m −2 in northern Puerto Rico following Grilli et al. (2010). We limited the maximum moment magnitude to 8.0 considering the historical seismic events that led to tsunamis in Puerto Rico (Nealon and Dillon, 2001). We assumed that the LEN should be longer than the W ID, and the range of the LEN and W ID should follow the scaling laws introduced by Blaser et al. (2010). The range of the dip angle DIP and the depth of the top edge DEP were determined based on the characteristics of a subduction-interface earthquake that usually causes a tsunami. According to Thingbaijam et al. (2017), subduction-interface earthquakes occur between 10 • and 30 • dip angles and within a slip-centroid depth of 50 km. We assumed that the fault rupture occurred instantaneously, where the initial free surface displacement was calculated using the Okada equations (Okada, 1985). Nearshore bathymetry and onshore topography in the inundation zone were from the 3 arc-second National Geophysical Data Center data set (NGDC, 2005), while the 1 arc-minute ETOPO1 data set (Amante and Eakins, 2009) was used for the entire region (Fig. 2). Considering the grid size, the minimum and maximum AMR levels were set to 5 and 11, respectively. The bottom friction was parameterized using a quadratic drag law in which the bottom drag coefficient C f was set to 10 −4 . The numerical model was used to simulate two hours of tsunami propagation to ensure that complete inundation of the onshore areas was captured. The maximum envelope of the water level was interpolated bilinearly onto a regular grid (0.001 • interval). We excluded four simulations, which failed to finish the simulations because of instability issue, to build the TRRF. We obtained tsunami runup distribution R(x) by extracting the maximum inundation height along the coastline ranging from 67.100 • W to 65.620 • W . The tsunami simulations were conducted in a spherical coordinate system, but the TRRF was defined based on a Cartesian coordinate system. To align the coordinate systems, Vincenty's formulae (Vincenty, 1975) were used to convert the unit of the geometric point from degrees to kilometers. We set the origin at (18.450 • N , 66.400 • W ), and thus the positive LON (LAT ) refers to the longitude(latitude) of the epicenter that is farther east (north) than the origin.

RSM functions and N T p (x)
The RSM functions and the normalized topography run-up N T p (x) were derived as follows. We calculated the OS formula coefficients a and b by fitting the tsunami run-up distribution R(x) to the OS formula for each simulation (Eq. 2). Here, the OS formula coefficient c was fixed to zero because we set the longitude of the origin and the epicenter longitude of simulations identically (Eq. 7). We derived the RSM functions by fitting the RSM fault parameters to the OS formula coefficients a and b using second order polynomial models. The normalized topographic run-up N T (x) was calculated for each simulation following Eqs. 1-3. We derived N T p (x) by selecting the 50th percentile of N T (x) among all simulations (Fig. 9).

Fault parameter range for TRRF prediction
We set the fault parameter range for TRRF prediction as shown in Table 1. The range of six fault parameters (LAT , LEN , W ID, DIP , SLP , DEP ) was set to the same range as the tsunamigenic-earthquake scenarios used in the TRRF development. In order to avoid an extrapolation beyond the inference space of the RSM functions, we only considered cases where all LAT of the hypothetical faults fell within the range for LAT . The range of LON was set to the extent that the fault does not fall outside the region used in the numerical simulation. The strike angle is usually set in the direction tangential to the subduction zone (Gica et al., 2008), and thus we set the range of ST R to be from 50 • to 130 • . Even though some tsunamis are generated by strike-slip earthquakes (Heidarzadeh et al., 2017), most tsunamis are caused by thrust earthquakes. Following this characteristic of RAK, we set the range of RAK to be from 50 • to 130 • .

TRRF calibration
To apply the AP method, (1) the site-specific coefficients (α, β, and γ) of Eqs. 8 and 9 and (2) the number of hypothetical faults (n) must be defined in advance.
To determine the coefficient α, we simulated 80 additional cases (hereafter called STR cases) that were not used in building the TRRF. These additional cases had a fixed longitude of 66.400 • W , where 10 sets varying the RSM fault parameters were randomly selected. For each set, eight different strike angles between 50 • and 130 • were selected, at 10 • intervals, except 90 • . The rake angle was fixed to 90 • so that θ could be independent of β and λ could be fixed to 90 • (see Eqs. 8 and 9). The coefficient α was selected by minimizing TRRF error as represented by normalized root mean square error (N RM SE): where R p (x) is the tsunami run-up distribution predicted by the TRRF,R p (x) is the numerically simulated tsunami run-up distribution, and N is the total number of alongshore locations. For each case, we found the θ value that shows the minimum NRMSE in the range of 45 • and 135 • . We fixed the number of hypothetical faults (n) to 100, which was large enough to provide a convergent prediction. We set the coefficient α to 0.585 by fitting the ST R and the θ in Eq. 8 (Fig. 10).
To determine the coefficients (β, and γ), we simulated 80 additional cases (hereafter called RAK cases) where all fault parameters but the rake angle were set in the same way as the STR cases. Unlike the STR cases, the rake angle was set to the same value as the strike angle. For each case, we found the θ value that shows the minimum NRMSE in the range of 45 • and 135 • . At the same time, we found the λ value that shows the minimum NRMSE in the range of 90 • and 179 • (if RAK < 90 • ) or the λ value that shows the minimum NRMSE in the range of 1 • and 90 • (if RAK ≥ 90 • ). We set the coefficient β to -0.284 by fitting the RAK and the θ to Eq. 8 ( Fig. 11(a)). We set the coefficient γ to -0.754 by fitting the RAK and the λ in Eq. 9 ( Fig. 11(b)).
To determine the number of hypothetical faults (n), we revisited the RAK cases. For each case, we decreased the number of hypothetical faults (n) from 100 to 2 (Fig. 12). Then, we found the minimum value needed for convergence since the computational time increases as n increases. In this study, we set the n to 30, which shows less than 0.1% difference in NRMSE.

Accuracy
The accuracy of the TRRF was investigated by comparing TRRF predictions against the direct numerical simulations. Since the TRRF considers the earthquake fault parameters in three groups ((1) RSM fault parameters, (2) epicenter longitude, and (3) strike angle and rake angle), we systematically tested the accuracy of the TRRF as follows: -Test 1: RSM Fault Parameters. We simulated 100 additional cases in which the RSM fault parameters were randomly selected, while the epicenter longitude was fixed to 66.400 • W and both the strike and rake angles were fixed to 90 • . -Test 2: Epicenter Longitude. We simulated 100 additional cases in which the fault parameters were selected based on the following conditions. While both the strike and rake angles were fixed to 90 • , 10 sets of the RSM fault parameters were randomly selected. For each set, 10 longitudes were selected at a uniform interval in the range of 65.800 • W and 67.000 • W . -Test 3: Strike Angle and Rake angle. We investigated the RAK cases defined in Section 4. -Test 4: Overall Accuracy. We simulated 100 additional cases in which all fault parameters were randomly selected. Fig. 13 shows the comparison of the OS formula coefficients based on the 100 cases of Test 1. Note that these 100 cases were never used to derive the RSM functions. The x-axis is the OS formula coefficient obtained by fitting the numerical simulation result to the OS formula. The y-axis is the OS formula coefficient obtained by putting the six fault parameters (LAT , SLP , LEN , W ID, DIP , DEP ) to the RSM functions. The high-correlated results confirm that the RSM functions can predict the OS formula coefficients well. Fig. 14 shows the selected alongshore tsunami run-up predictions for each test. Comprehensive summary statistics of test results are presented in Table  2. Fig. 14(a) and (b) show the best case (minimum N RM SE) and the worst case (maximum N RM SE) of Test 1, respectively. In both cases, the TRRF prediction followed the overall trend of the numerical simulation result well. However, there are a few localities where the TRRF did not predict the runup well such as the run-ups near 66.2 • W in the worst case. Fig. 14(c) and (d) display the Test 2 results where all the fault parameter conditions were the same except the epicenter longitude. The results show that the TRRF can effectively capture the influence of the epicenter longitude. Summary statistics, which have only increased slightly compared to Test 1, confirm that Eq. 7 is valid (Table 2). Fig. 14(e) and (f) present the Test 3 results in which all fault parameter conditions are the same except for the strike and rake angles. Note that the TRRF predictions strongly align with the numerical simulation results while capturing the asymmetrical shape of the tsunami run-up distribution. Summary statistics, which increased only slightly from Test 1, confirm the  [66.400, 19.058, 90,18, 90, 62.549, 28.405, 2.57, 23.686] [66.400, 19.518, 90, 10, 90, 88.456, 30.107, 2.96, 17.456] [65.800, 19.264, 90, 14, 90, 82.580, 27.544, 2.71, 17.185] [67.000, 19.264, 90, 14, 90, 82.580, 27.544, 2.71, 17.185] [66.400, 19.149, 50, 20, 50, 103.693, 51.944, 2.56, 27.802] [66.400, 19.149, 130, 20, 130, 103.693, 51.944, 2.56, 27.802] [66.856, 19.378, 120, 13, 111, 77.682, 39.220, 3.37, 23.018] [66.203, 19.244, 126, 16, 127, 119.959, 38.144, 3.99, 10.339]  performance of the AP method (Table 2). Fig. 14(g) and (h) show the best case (minimum N RM SE) and worst case (maximum N RM SE) of Test 4, respectively. Both examples have a few localities where the TRRF did not predict the run-ups well, but the overall trend of the TRRF predictions agrees well with the numerical simulation results.

Computational time
The efficiency of the TRRF was investigated by comparing the computational time between the physics-based numerical model and the TRRF. When the physics-based numerical model (Basilisk) was used to predict the tsunami runup distribution, the computational time was about one hour (24 CPU hours) on average for each scenario (24 cores, OpenMP, Intel Xeon E5-2680v3). On the other hand, when the TRRF was used to predict the tsunami run-up distribution, the computational time was only 0.01 CPU second per scenario (desktop computer with one core, Intel I7-7700). The TRRF's CPU time is nearly 9 million times shorter than that of numerical simulation. The difference in computational time between the TRRF and the numerical model would be even greater given higher resolution grids and larger geographic areas than those used in this study.

Discussion
The performance of TRRF was investigated based on total 380 additional simulations in Section 5. When the TRRF predictions are compared against the direct numerical simulations, it is clear that the TRRF can produce reliable run-up predictions over real topography, given the computational time. However, as shown in Fig. 14, even though the TRRF predicts the leading order of tsunami run-up distribution well, there are a few localities where the difference of the run-up is more than two-fold. We found that these localities are correlated to the uncertainty (or the range of percentiles) of N T p (x) (Fig. 9). For example, a large uncertainty was commonly found in places with complex topography, such as areas surrounded by mountains (e.g., 65.735 • W ), areas containing a river (e.g., 66.955 • W ), steep cliffs (e.g., 66.444 • W ), and coastal dunes (e.g., 66.239 • W ). Even though it is difficult to fully interpret the physics behind the normalized topographic run-up N T p (x), we think that this high uncertainty may be attributed to the nonlinear behavior of the tsunami wave as it propagates and inundates complex topography. In its present form, the TRRF does not directly consider potential nonlinearities between the source and topographic run-up components in the hypothesis that the tsunami run-up distribution can be expressed as a sum of the source and topographic run-ups (Eq. 1). Future studies should investigate ways to account for the uncertainty to improve the accuracy of the TRRF approach. Moreover, future studies should expand the applicability of the TRRF by considering the following limitations. One is that the TRRF is only applicable to uniform slip distribution. Several studies have shown that tsunami prediction can vary depending on heterogeneous slip models even when the earthquake magnitude is the same (Geist, 2002;Li et al., 2016;Ruiz et al., 2015). The other limitation of the TRRF is that it is only applicable to tsunamis generated by seafloor displacements associated with earthquakes. After earthquakes, landslides are the second most common cause of tsunamis (Harbitz et al., 2014). Moderate earthquakes do not always cause tsunamis themselves, but they can, in some instances, trigger large landslides that result in tsunamis (Uri et al., 2009). Though landslide-generated tsunamis are rare, a single occurrence can cause substantial damage and loss of life. For example, in 2017, a landslide-generated tsunami off the western coast of Greenland flooded several villages and resulted in casualties (Paris et al., 2019). A recent study also revealed that the 2018 Indonesian tsunami, which claimed more than 2,000 lives and severely damaged coastal communities, was caused by the combination of an earthquake and a landslide (Sassa and Takagawa, 2019). Several other key elements would merit attention in future studies. For example, the arrival time and inundation distance are as important to consider as the run-up. A high tide could enhance tsunami inundation, while a receding tide could dissipate tsunami energy (Tolkova et al., 2015;Zhang et al., 2011). Likewise, a modest amount of sea-level rise could dramatically impact the tsunami run-up distribution (Li et al., 2018).

Conclusions
In the present study, we presented a new methodology, called TRRF, that can predict the alongshore run-up distribution from a near-field tsunami. We adopted the OS formula and developed what we call the AP method to reduce the number of simulations to build the TRRF. The tsunami run-up distribution was decomposed into source run-up and topographic run-up, that source run-up can be modeled by earthquake fault parameters, and that normalized topographic run-up is associated with local topographic characteristics. Using the northern region of Puerto Rico as a case study, the performance of the TRRF was investigated based on total 380 additional simulations. The results showed that the TRRF can produce rapid near-field tsunami run-up predictions over real topography. We expect that future applications of the TRRF will have the potential to save lives and promote resiliency of coastal communities. data, is enough to fit the fault parameters to the OS formula's coefficients a and b (Fig. 15).