Natural Hazards and Earth System Sciences High resolution tsunami inversion for 2010 Chile earthquake

Abstract. We investigate the feasibility of inverting high-resolution vertical seafloor displacement from tsunami waveforms. An inversion method named "SUTIM" (small unit tsunami inversion method) is developed to meet this goal. In addition to utilizing the conventional least-square inversion, this paper also enhances the inversion resolution by Grid-Shifting method. A smooth constraint is adopted to gain stability. After a series of validation and performance tests, SUTIM is used to study the 2010 Chile earthquake. Based upon data quality and azimuthal distribution, we select tsunami waveforms from 6 GLOSS stations and 1 DART buoy record. In total, 157 sub-faults are utilized for the high-resolution inversion. The resolution reaches 10 sub-faults per wavelength. The result is compared with the distribution of the aftershocks and waveforms at each gauge location with very good agreement. The inversion result shows that the source profile features a non-uniform distribution of the seafloor displacement. The highly elevated vertical seafloor is mainly concentrated in two areas: one is located in the northern part of the epicentre, between 34° S and 36° S; the other is in the southern part, between 37° S and 38° S.


Introduction
The 2010 Chile earthquake occurred off the coast of the Maule Region on 27 February 2010, registering a magnitude of 8.8 on the moment magnitude scale (USGS).In this event, the segment of the fault zone was estimated to be between 500 km to 700 km long with an average slip about 10 m (USGS, BGS).A precise GPS measurement (http://researchnews.osu.edu/archive/chilemoves.htm) showed that several cities south of Cobquecura were raised by up to 3 m at least.Lorito et al. (2011) and Fritz et al. (2011) further estimated that the maximum surface Correspondence to: T.-R.Wu (tsoren@ncu.edu.tw)deformation was about 4 to 5 m.On the other hand, the tsunami signal recorded by the buoys and tidal gauges around the epicenter might open us a window to understand the seafloor displacement and earthquake characteristics.
Tsunami is a gravity wave generated mainly by earthquakes and landslides.A strong earthquake occurring in the ocean is often accompanied by tsunamis.Because the wavelength is much longer than the water depth, tsunami is categorized as a pure linear wave which obeys the shallow water equation (SWE) in wave propagation.In addition, the wave propagation speed is much slower than the seafloor movement, therefore the initial tsunami profile, or the tsunami source, can be treated as the seafloor vertical displacement (Wu et al., 2009;Megawati et al., 2009;Huang et al., 2009;Saito and Furumura 2009).Once the tsunami initial freesurface profile is inverted from the wave gauges, some fault parameters can be determined.Compared to the conventional seismic inversion method, the tsunami inversion requires only wave signals and bathymetry data, which can be accessed relatively easily.This idea is not new; Satake (1987) proposed a tsunami inversion method to infer the seismic parameters from the tidal gauge data.His method was further used on the tsunami forecast (Wei et al., 2003;Yamazaki et al., 2005) and exploration of earthquake events with sea surface deformation (Pires and Miranda 2001; Baba et al., 2002Baba et al., , 2005;;Satake et al., 2005;Piatanesi and Lorito 2007).However, those implementations used very large unit sources such that the details of the tsunami profiles could not be described.With the advance of fast computation and modern tidal gauge facilities, the small unit inversion method becomes feasible.Saito et al. (2010) used small unit sources to study the 2004 earthquake off Kii Peninsula with a reasonably good result.In their paper, they emphasized the importance of dispersion effect to the smallscale tsunamis.However, no performance tests and model validations were presented in their studies.In this paper, we introduce the small unit source tsunami inversion method (SUTIM) to rebuild the fine structure of tsunami source.Instead of using the Gaussian distribution (Saito et al., 2010)  the top-hat small unit sources are adopted to maintain the mass conservation.To save computational costs, moreover, a grid shifting method is proposed.To ensure the stability of the inversion result, the smoothing constraint is applied.In the present model, the wave dispersion effect is assumed to be small and negligible.To evaluate the performance of SUTIM, a fictitious earthquake event is created as the benchmark problem.Last, the SUTIM is utilized to study a real event: the 2010 Chile Earthquake.The result is presented and discussed in this paper.

Shallow water equation and COMCOT Numerical model
The wave motion can be described by the potential flow theory.By integrating the equations along the water depth, the shallow water equation can be derived (Dean and Dalrymple, 1991).Equations ( 1) and ( 2) are the nonlinear SWE in a non-dimensional form.
Equation ( 1) is the continuity equation and Eq. ( 2) is the momentum equation, in which h represents the local water depth, η is wave height, t is time, u is the horizontal velocity vector, and ε is the nonlinearity defined as ε = a/ h, a is the wave height.If ε < 1/20, the wave is categorized as a linear shallow water wave and the terms with ε can be omitted.Another important parameter used to define the dispersion effect is dispersion coefficient, µ, defined as µ = h/λ, in which λ is wavelength.If µ < 0.05, the wave is a nondispersive shallow water wave and can be described by SWE (Dean and Dalrymple, 1991).In general, an earthquake tsunami propagating in the deep ocean is a typical linear and nondispersive wave.Both nonlinear and dispersion effects can be ignored.
To calculate the tsunami propagation, a SWE code named COMCOT is utilized to perform the calculation.COMCOT is able to solve both linear and nonlinear shallow water equations in Cartesian or Spherical coordinate systems.It has been widely used to study many historical tsunami events, 18  such as the 1992 Flores Islands tsunami (Liu et al., 1994(Liu et al., , 1995)), 2004 Indian Ocean tsunami (Wang andLiu, 2006, 2007), and 2006 Pingtung tsunami (Wu et al., 2008;Chen et al., 2008).The accuracy reaches the 2nd order in both spatial and time domains.COMCOT also supports the nested grid system that allows the finer grid to be placed on a coarser grid to increase local resolution.

Tsunami inversion method
The energy of a tsunami wave can be accumulated linearly because of its linear characteristic.In other words, the wave energy can be divided into many small units and superposed back on the original one, i.e. each element of vertical displacement is deduced numerically without referring to any dislocation formulation.Upon the linear superposition theory and the least square method, SUTIM is developed.The algorithm used in SUTIM is shown in Fig. 1 and below: 1. Select waveforms from buoys or tidal gauges.Judging by the data quality and azimuthal distribution, we choose the proper wave gauges to perform the inversion.The gauges close to the source region with deeper local water depth and less site effects are preferred.2. Decide the inverse region.Although this region has to cover the entire possible tsunami source area, using smaller coverage can save computational time and constrain the uncertainty of inversion.
3. Place a series of unit sources over the inverse region.Each unit source shares the same wave height, which is a unit.
4. Calculate the wave propagation of each unit source, and record the computed waveforms at each gauge.
5. Adjust the initial wave height by multiplying the weighting factor, and sum the adjusted gauge results together to have the best fit to the real gauge record.We adopt the least square method to obtain the best-fit result.
6. Apply the grid shifting method to increase the resolution of the tsunami source.
In the following sections, the detailed algorithm of SUTIM will be described.

Least square method
In SUTIM, the numerical gauges record the computed waveforms from each unit source.By comparing the results with field data, a best-fit weighting matrix can be determined.The weighting matrix is a factor matrix for the unit source array and the multivariable least square method is utilized.We set matrix a as the factor matrix of the unit source array.The linear matrix-vector equation is shown in Eq. ( 3).The first lower-index indicates the data number.The second index indicates the number of the unit source.Vector b is the observation gauge data.Vector x is the weighting factor, which will be solved by the least square method.Constant m is the total number of the unit source, and n is the total number of the gauge data.
Equation ( 3) can be rewritten in a matrix form: Based on the least square method, the weighting vector x can be acquired by solving Eq. ( 5): Equation ( 5) is the standard formula of the least square method used to invert the equation.However, this method has to be modified for the multi-gauge and multi-event cases.
Because the least square method suffers from the unstable and non-physical features, in which case a smoothing constraint is required in the following sections.

Multi-gauge inversion
In most of the real tsunami events, more than one wave gauge or buoy record will be used to perform the tsunami inversion.
When multiple gauges are employed in the inversion method, Eq. ( 4) can be extended as follows: in which m represents the total number of combined data.

Multi-event inverse
In some earthquake events, multiple earthquakes happened in a short period.The 2006 Pingtung earthquake doublet is one example.Two earthquakes with similar magnitude occurred within 8 min.In such case, both events shall be included in the tsunami inversion.In the multi-event approach, we create unit sources with multiple events, and then calculate all the wave propagation of each unit source.
In the case of the multi-event, the observation data b has the combined influence from all events, and therefore the calculated wave height shall be the summation of each multiunit source: After composing all the unit source results, the modified linear matrix-vector equation can be derived.Equation (7) shows the example of two events.
in which the elements from a 11 to a 1n represent the gauge data of the first event, and elements from a 1 n + 1 to a 1n+n represent the gauge data of the second event.
x 1 to x n are the weighting factors of the first event.x n+1 to x n+n are the weighting factors of the second event.b 1 to b M are the observed data containing signals of both two events.

Smoothing constraint
Based on the seismic theory and field observations, the surface of the ground movement should be smooth (Hartzell and Heaton, 1983).However, the seafloor profile inverted directly from the tsunami inversion might be bumpy and non-physical, a smoothing constraint is required.Adding a smoothing constraint as a damping function has been widely used in the seismic inversions (Yen et al., 2008;Lee et al., 2008) and some tsunami inversions (Baba et al., 2005;Saito et al., 2010).The smoothing constraint constrains each tsunami unit source to have a relationship with its neighboring unit source as: in which x s is the weighting factor of one specific unit source, x i represents the weighting factor of neighboring units of the specific unit.k is the total number of neighboring unit source on x s , and swt is the smoothing weight scalar.Thus, each unit has defined relationship with its neighboring units.For example, considering one specific tsunami unit source 33 with the neighboring units 32, 34, 23, and 43, the equation of the smoothing constraint can be written as: in which x 33 represents the weighting factor of unit 33, and unit 32, 34, 23, 43 are x 32 , x 34 ,x 23 , and x 43 .The swt is the smoothing strength.Equation ( 9) can be added to the matrix a as shown in Eq. ( 10): in which i and j are the IDs of the unit source, M is the length of data, and swt is the smoothing strength.With this constraint, each unit source builds up the relationship with its neighboring units, and the result will be smoothed.

Grid Shifting Method (GSM)
In order to reproduce the detail of the source distribution, a large number of unit sources have to be applied to the possible tsunami source region.Each unit source requires 20 grid points to resolve the wave propagation.Therefore, double resolution requires quadruple unit sources, not to mention using a finer grid size and smaller time steps.The calculation is very time consuming even on modern computers.In   this study, we propose the grid shifting method to enhance the source resolution with limited extra CPU load.
The grid shifting method is described as follows: 1. Create the first layer of unit source array as previously described.
2. Create the second layer of unit source array identical to the first one.However, we slightly shift the second unit source array with 1/2 grid spacing (Fig. 2).As the result, there is an overlapped area where the resolution is 1/2 of the original one.
3. Calculate the wave propagation of both unit source arrays, and perform the least square method to find the weighting factors on both arrays.
4. Average both weighting factors on each overlapped unit source.
The idea of grid shifting method is that the two-layer unit sources shall be able to capture more detailed information than just one-layer ones.Shifting the unit sources by half grid spacing indicates that the sampling resolution is doubled.Because the size of the unit source is larger than the sampling resolution, we average both layers to obtain the final solution.This might raise an issue of over-smoothing problem.However, the grid-shifting method indeed provides detailed information of tsunami source profile.The good performance of this method can be seen in the following sections.
As the length of averaged result on each unit source is 1/2 of the original unit source, the inversed tsunami source has finer resolution and the calculation requires only the double of the computational load.We shall demonstrate the performance of GMS in the following sections.

Restriction in SUTIM
Several assumptions are made for SUTIM; we shall list all the limitations while using it.
1.The tsunami has to be a linear water wave, which obeys a/ h < 1/20.
2. The tsunami has to be a long water wave, which follows h/λ < 0.05.
3. The tsunami source is large and the slip duration is short.One of the important assumptions used in SU-TIM is that the vertical component of seismic motion is the same as the tsunami source.However, it is not always true.Saito and Furmura (2009) proposed the limitation to this assumption.The criteria are λ 13h and T λ/(8c), in which T is slip duration, and c is the phase velocity of the tsunami.In other words, the length of source has to be much longer than the water depth, and the duration of seismic motion has to be very short to ensure this assumption is correct.

Performance Test of SUTIM
In order to examine the performance of SUTIM, a fictitious tsunami event is created as shown in Fig. 3a.The size of the fictitious source region is 172.7 km by 155.4 km, located at 100 km off Taiwan.The source profile is generated by the half-space semi-elastic model (Okada, 1985) with the fault parameters listed on Table 1. Figure 3 shows the bathymetry profile as well as the gauge location.The wave gauges are all placed in the southern part of Taiwan.The ill-deployed gauges and complex bathymetry raise the difficulty for the inversion.The maximum surface elevation of the source is 0.049 m and −0.243 m with local water depth about 2500 m (Fig. 4a).The fictitious source profile is smooth, yet the result from SUTIM will be a stair-like profile, which adds some uncertainties to the inversion.Figure 4b is the inversion result directly from the least square method.Although the negative and positive tsunami sources as well as the strike angle can be roughly identified, the coarse resolution makes recognition difficult and the surface bumpy.
Figure 4c is the result after the smoothing constraint is applied.The inversed tsunami source is smoother compared to Fig. 4b.However, the low-resolution issue still exists.Figure 4d shows the results after the adoption of both smoothing constraint and GSM.The source profile is much cleaner than that in Fig. 4b and c.The surface is smooth and the resolution is high enough to easily identify the elevation and the strike angle.The inversion result is very close to the fictitious source.Figure 5 shows the comparisons on the cross-section in Fig. 4a.The cross-section is located at 21.5 • N and extends from 118.5 • E to 120.4 • E. As seen in Fig. 5b, although SU-TIM fails to predict the flat surface between 118.5 • E and 119.0 • E, the steep slope near 119.4 • E can still be easily identified.This results from the low resolution of the unit source.Figure 5c shows the smoothed result, and Fig. 5d shows the final result with smoothing constraint and GSM.The result indicates that, with the smoothing constraint and GSM, SUTIM can improve the quality of the inversion results.

2010 Chile earthquake
Before applying the SUTIM to the event of 2010 Chile earthquake, we shall inspect the restrictions that require parameters, a, h, λ, and T .The possible wave height a is assumed to be 3 m estimated from the GPS measurement in the south of Cobquecura.The local water depth h varies from 500 m to 3000 m.To ensure the inversion quality, h = 500 m is adopted.The wavelength or the length of the source λ can be deduced from the distribution of the aftershocks.In this case, it is estimated to be 100 km.The source duration or the rise time T can be estimated from slip divided by slip velocity.The slip is about 10 m from BGS, and the slip velocity is about 1 m s −1 in most of the earthquake events.The source duration T is therefore estimated to be 10 s.
Based on the above mentioned parameters, the characteristics of the tsunami source can be determined: a/ h = 0.006, smaller than 1/20; h/λ = 0.005, which is smaller than 0.05; λ = 100 m h = 500 m; and T = 10 s λ/(8c) = 178 s.We conclude that SUTIM is suitable for studying the 2010 Chile earthquake.

Parameter predetermination
To employ SUTIM, the possible source region has to be determined a priori.Based on the seismic report from USGS website and aftershocks reported from Global CMT, the length is estimated to be 500 km.The rupture width is assumed to be 100 km upon the announcement of USGS.According to the half-space semi-elastic model (Okada, 1985), the scaling law, and fault parameters from Global CMT, the tsunami source can be verified (Wu et al., 2008;Wu and Huang, 2009).Table 1 lists the fault parameters used to generate the tsunami source for the 2010 Chile Earthquake.We refer this method to the traditional source-determination method.Figures 6 and 7 show the tsunami source and the gauge comparisons, respectively.We can see that the tsunami source from the traditional source determination is smooth and uniform spatially.However, the source profile does not seem natural.The calculated waveforms shown in Fig. 7 are simulated by COMCOT.In Fig. 7, although the arrival time of the first peak is predicted correctly, the detailed waveforms deviate from the field data.This indicates that the source location is correct, yet the source profile is not.

Source inversion using SUTIM
We apply forty-five unit sources to the first layer and fiftythree unit sources to the second layer for GSM.The size of each unit is a 40 × 40 square.Six gauges are chosen from GLOSS (The Global Sea Level Observing System) and one is from NDBC DART station as shown in Fig. 6. Figure 6 also shows the computation domain and unit sources.Figure 8b shows the inversion result with GSM and smooth constraint.In the figure we can see the seafloor displacement distributes non-uniformly.Two primary clusters of positive elevation can be observed.One cluster distributes along the coastline between 35 • S and 36 • S, which is the asperity of this event.The maximum elevation is 3.63 m, which is very close to the GPS measurement.The other is weaker and distributed in the south part between around 37.2 • S with local maximum elevation 1.87 m.
Our sea-floor vertical distribution is very similar to the slip distribution proposed by Tong et al. (2010), using a GPS and InSAR data set.Our result is also close to Lay et al. ( 2010), Delouis et al. (2010), andLorito et al. (2011) in terms of the spatial pattern, especially for the locations of both primary slip maximum near 35 • S, 72.5 • W and the secondary slip maximum near 37.2 • S, 73.5 • W. We as well present the inversion result without GSM in Fig. 8a to point out that a better inversion result can be expected with GSM In addition to comparing the spatial distribution to the results from seismic data, we also compare the inverted profile with the distribution of aftershocks (Fig. 9a) and the tsunami source determined form the traditional method (Fig. 9b).The profile of the seafloor elevation agrees well with the distribution of aftershocks (Fig. 9a).The strike angle and length from SUTIM are similar to those obtained from the traditional method of source determination (Fig. 9b).
Because the exact source area cannot be determined before obtaining the final inversion result, we tend to choose a large source zone to cover most of the tsunami source area.However, it is obvious that the source zone in the 2010 Chile case extends too much in the northwest direction and reducing it would have allowed better constraining in the inversion.
We further examine the calculated waveforms with the filed data.Figure 10 shows the time-history free-surface at each gauge location.Compared to Fig. 7, the result obviously has better fit to the field data in terms of the arrival time, wave period, and waveform.This indicates that the high-resolution result from SUTIM accurately predicts the detail of the tsunami source.
In the current study, the SUTIM does not take into account a triggering sequence along the rupture area for megaearthquakes.However, for a large rupture like the 1200 km along-strike of the 2004 Sumatra event, rising tsunami time from the epicenter to the fault might be crucial to recover the entire tsunami phase.The initial tsunami is not one quasistatic wave but a wave sequence.This can be improved by using a multi-event procedure in the future.
As far as direct inversion with tsunami information-only is concerned, a further validation of the procedure is needed for testing the tsunami modeling itself and the robustness of the simulation.One possibility is through the runup observations.This might be an excellent opportunity to validate our derived source.It would be interesting to state that the source we have is unique.However, such validation might be difficult with the lack of an accurate bathymetric field of the area.Another possibility is a cross validation of available sources constructed with geophysical data.This can be done in the near future.

Conclusions
Without considering earthquake parameters, tsunami waveform inversion serves as an alternative way to reconstruct the seafloor vertical displacement.With high-resolution inversion, the strike angle and the detailed source profile can be depicted.In this study, we introduce a tsunami inversion method named small unit source tsunami inversion method (SUTIM).A complete inversion procedure is presented.This procedure includes the least square inversion method, multigauge method, and multi-event method.We propose the grid shifting method (GSM) to increase resolution with limited extra CPU effort.We also introduce the smoothing constraint to enhance accuracy.We demonstrate the performance of SUTIM by simulating a fictitious earthquake event in which the seafloor motion is known.Good SUTIM performance can be seen.
We then perform SUTIM to study the 2010 Chile earthquake event.By examining the resections of SUTIM, we discover that SUTIM is suitable for studying the 2010 Chile earthquake event.The inversion result shows a stronger elevation appears along the coastline, which agrees well with the GPS measurement and the distribution of aftershocks.
From the performance test and validation, GSM indeed raises the resolution of the inversion.However, the mathematic derivation fails to be derived.This shall be discussed in the future studies.
Figure 1.The procedure of SUTIM.

Figure 2 .
Figure 2. The sketch of grid shifting.The 1 st layer of unit source array is presented in solid lines.The 2 nd layer is presented in dashed lines.The shaded area is the final resolution by applying the grid shifting method.

Fig. 2 .
Fig. 2. The sketch of grid shifting.The 1st layer of unit source array is presented in solid lines.The 2nd layer is presented in dashed lines.The shaded area is the final resolution by applying the grid shifting method.

Figure 4 .Fig. 4 .
Figure 4.The tsunami source profiles.(a) Fictitious tsunami source.(b) by Least square 3 inversion only.(c) with smoothing constraint.(d) with smoothing constraint and GSM.The 4 white line in (a) shows the cross-section used in Figure 5. 5 6

Figure 5 .
Figure 5.The cross-section profiles of the tsunami source.(a) Fictitious tsunami source.(b) by Least square inversion.(c) with smoothing constraint.(d) with smoothing constraint and GSM.

Fig. 5 .
Fig. 5.The cross-section profiles of the tsunami source.(a) Fictitious tsunami source.(b) by Least square inversion.(c) with smoothing constraint.(d) with smoothing constraint and GSM.

Figure 6 . 7 Fig. 6 .
Figure 6.The computational domain, bathymetry, gauge location (marked in white triangle), 3 epicenter of the 2010 Chile earthquake (marked in red asterisk), tsunami source from Okada 4 model (drawn in contour lines with 1m interval), and the unit sources (drawn in small red 5 squares).The color bar shows the bathymetry in meters.6 7

Figure 8 .
Figure 8. Inversion result by SUTIM of 2010 Chile earthquake.(a) by Least square inversion and smoothing constraint.(b) with GSM.Color bar indicates the vertical displacement of the seafloor in meters.

Fig. 8 .
Fig. 8. Inversion result by SUTIM of 2010 Chile earthquake.(a) by Least square inversion and smoothing constraint.(b) with GSM.Color bar indicates the vertical displacement of the seafloor in meters.
Figure 9. (a) The inversion result with SUTIM of vertical displacement and aftershocks of the 3 February 27, 2010 event.(b) the inversion results with SUTIM in gray and from traditional 4 tsunami source determination method presented in counter lines.The blue circle is the 5 epicenter location reported by Global CMT. 6 7

Table 1 .
The fault parameters of fictitious earthquake and 2010 Chile earthquake.