CALIBRATION AND VALIDATION OF METHODOLOGY THAT INCORPORATES VARIABILITY IN COSEISMIC SLIP SPATIAL DISTRIBUTION

MOTIVATION The estimation of inundation zones is one of the main objectives of tsunami early warning systems but is compounded by several challenges, ranging from the high computational cost to the inherent uncertainties present. The latter stem in part to the incapacity to establish accurately and with certainty the parameters that characterize the coseismic slip during an earthquake. However, here it is assumed that it is possible to incorporate an estimate of the uncertainty and its consequent variability in order to offer a forecast that better portrays the hazard and, eventually, lead to the protection and mitigation of such disasters of natural origin.


MOTIVATION
The estimation of inundation zones is one of the main objectives of tsunami early warning systems but is compounded by several challenges, ranging from the high computational cost to the inherent uncertainties present. The latter stem in part to the incapacity to establish accurately and with certainty the parameters that characterize the coseismic slip during an earthquake. However, here it is assumed that it is possible to incorporate an estimate of the uncertainty and its consequent variability in order to offer a forecast that better portrays the hazard and, eventually, lead to the protection and mitigation of such disasters of natural origin.
The spatial distribution of the slip of an earthquake is a determining factor when estimating a tsunami propagation (Geist, 2002(Geist, , 2013. However, it is not possible to determine with absolute certainty its spatial distribution, even with the benefit of time, owing to epistemic and aleatoric uncertainties (Cienfuegos, 2018). By incorporating variability in the slip distribution in the rupture zone, this uncertainty could be incorporated into the tsunami propagation and inundation.
Stochastic models are used to generate predictive scenarios of slip distribution (D(x,y)). Mai & Beroza (2002) show that these types of models provide a more accurate approach to parameterize the complexity of the slip of an earthquake, in this way epistemic uncertainty can be considered. The slip distribution can be characterized as a spectral power density in the wavenumber domain, where the parameters of the random field are related to variables such as the magnitude and geometry of the fault (Lavallée & Archuleta, 2013). The use of a spectral density model allows the generation of scenarios to preserve this condition, and thus to be able to simulate the movement of the ground with a stochastic realization of the sliding distribution. Typically, the spectral power density is used in conjunction with a random phase to generate scenarios that retain the magnitude but are not constrained regarding the spatial distribution (e.g. Goda 2014Goda , 2015. This is suitable for long term hazard assessment, but is less appealing for early warning hazard assessments where a baseline earthquake should be considered.
Here, a Phase Variation Methodology is introduced to provide a simple method to generate a set of stochastic scenarios, that would have a mean expected value consistent with the baseline earthquake. This is achieved by a single free parameter α controlling the randomness.
However, this parameter has no physical basis, so in the first instance, it is not possible to define its value.
The main objective of this work is to determine whether is it possible to use the Phase Variation Methodology in an early warning context. For this, it is sought to determine the value of α that best reproduces the inundation observed for the Maule event (2010), and that is capable of reproducing the uncertainty that a tsunami prediction process implies. Then, the calibrated methodology is validated for de Illapel event (2015) in a forecast mode.

METHODOLOGY
To evaluate the effects of the tsunami, the inundation to the Central Chile cities of San Antonio, Constitución and Talcahuano that were affected by the earthquake and tsunami of Maule in 2010 have been studied. At least 19 solutions have been published pertaining the source model for that earthquake. Cienfuegos et al. (2018), showed that the median inundation ( 50% ) calculated these 19 source models is optimally adjusted to the flooding observed in 2010. Therefore, this is treated as the target inundation that is sought to be reproduced by the Phase Variation Methodology and it will be used as a basis to validate the results obtained.
To generate the stochastic scenarios, the observed slip distribution is transformed into its power spectral density. This is then modeled by a Von Kármán spectral power function, given by: (1) where k is the wave number defined by = √ 2 2 + 2 2 , and are horizontal and vertical wave number, respectively. and are the correlation lengths that control the absolute level of the power spectrum for small wave numbers, and H is the Hurst number exponent (Mai & Beroza, 2002). P(k) with this the amplitude spectrum model is constructed.
Unlike previous methods, the phase, ϕ, of the observed power spectrum is computed and retained. However, the modeled phase corresponds to a perturbation of the observed phase: where m and n are along-strike and along-dip coordinates in wavenumber space. Each realization of the phase ϕ is calculated as the sum of ̅ , which is the average phase from the reference model, and which represents the perturbation. The latter is given by: Where N(0,1) is a normally distributed random field, and α is the free parameter which determines the range of variability. Then, with the model amplitude spectrum and the phase spectrum the random fields are obtained. Finally, applying Inverse Fast Fourier Transform, the spectral scenarios are converted to the spatial domain, obtaining a set of N stochastic scenarios. By virtue of the central limit theorem, a large number of N would have an expected value that is approximated to the observed scenario.
It is acknowledged that this is an algorithmic approach that does not take into account the underlying physics of the rupture, but its simple implementation is hypothesized is a good compromise for early warning.
However, it is necessary to determine which value of α best reproduces the observed flood and variability.
A set of 10 different values are tested for α ( = ⁄ , = 1 … 10). For each one of these values, 50 slip scenarios are generated using the Phase Variation Methodology.
Subsequently, the representative inundation is obtained in each case, which is given by the median of the 50 scenarios generated, and also is calculated the variability provided by each case, which is characterized as: Where 95% and 5% are the Inverse Cumulative Density Function at 95% and 5%, respectively. Representative hydrodynamic variables are extracted from these results: maximum run-up (z(x,y)) and maximum amplitude (η(x,y)). These results are compared with the observed flooding and with the variability that exists between the inversion models presented in Cienfuegos et al. (2018).
To determine which value is the most optimal, the maximum likelihood method is used (Simmon, et al., 2017). This consists of assigning an adjustment score (BSS) to each prediction model, which depends on the similarity of the results obtained with those to be reproduced. The BSS score is given by: Where corresponds to the modeled inundation and is the observed inundation. This value must be normalized so that the sum of all scores equals 1, then: The value that maximizes the score is selected for both, the inundation and the expected uncertainty. After defining and selecting the value of α that maximize the score, the methodology is validated by generating stochastic scenarios through the Phase Variation Methodology for the Illapel event (2015) using the best value of α. The hypothesis is that this value would be able to reproduce the median inundation and variability given by the inversion models for different events.

RESULTS
To determine the level of variability that the value of α gives to the results, the median inundation ( 50% ) and the uncertainty level (I) are calculated for San Antonio, Constitución and Talcahuano, using Tsunami HySea software for the tsunami modeling, where 4 hours of propagation are calculated with high-resolution numerical grids.
From Figure 1.a) is observed that as the value of α decreases, so does the inundation height. This is explained by high values of α allows the generation of scenarios that differ to a greater extent from the original, being able to produce extreme events inconsistent with the actual source.
Consequently, Figure 1.b) shows that the variability decreases as the value of α also decreases, while for higher values of this parameter the variability is even greater. This is explained because by increasing the value of α the scenarios that are generated differ to a greater extent from each other, while by decreasing the parameter the methodology is restricted to generate scenarios similar to the original, resulting in very little variability between scenarios. The goal is then, to find an adequate balance between these extremes.
Figure 2.a) shows the L BSS scores obtained for each case of α. It is observed that for lower values it is possible to reproduce to a great extent the observed flood, however, it is not possible to reproduce the variability.
Furthermore, using values of α close to π 6 ⁄ , it is possible to satisfactorily reproduce both the observed inundation and the uncertainty level, as it is shown by Figure 2.b).
Finally, the methodology is validated using a value of = 6 ⁄ and = 5 ⁄ , the results using = are also obtained, this with the main of determining the sensitivity of the Coquimbo bay to the parameter. Figure 3.a) and 3.b) shows the median inundation and level uncertainty, respectively, for Coquimbo bay considering = 6 ⁄ . It is observed that both cases the maps obtained are similar to those expected.
Also, as seen in Figure 3.c), the adjustment scores for the Illapel event in 2015 are shown, it is observed that for = 6 ⁄ and = 5 ⁄ an optimal fit is obtained, so it is possible to validate the methodology in the first instance.

CONSLUSIONS
The sensitivity of the Phase Variation Methodology to the parameter α has been identified, which is responsible for assigning variability to the initial deformation that a tsunami will produce. For this reason, it is important to have knowledge of the values that could be assigned to the parameter α in order to obtain uncertainty in the results without moving away from the real inundation. A range that goes from = 6 ⁄ and = 5 ⁄ is proposed, since these values provide an acceptable fit between the ⁄ . The black dotted line represents the inundation and variability obtained from the inversion models. c) Adjustment score for the median flood (blue) and uncertainty level (orange) for the Coquimbo bay, considering α = π 6 ⁄ , α = π 5 ⁄ and α = π. The red dotted line indicates the position for α = π 6 ⁄ . models and the actual inundation, in addition to providing the expected epistemic uncertainty. However, the exact value that will be assigned to the parameter will depend on the context for which it is required, since for reliable sources of initial slip it might not be necessary to assign great variability to the results, while for fast inversion models or that there is few information of the rupture, as occurs in early warning systems, greater uncertainty would be expected, so it would be necessary to assign a higher value to α, however, it is not recommended that it be greater than = 4 ⁄ , since unrealistic scenarios could be generated.