Consistent Estimates of Tsunami Energy Show Promise for Improved Early Warning

Abstract

Early tsunami warning critically hinges on rapid determination of the tsunami hazard potential in real-time, before waves inundate critical coastlines. Tsunami energy can quickly characterize the destructive potential of generated waves. Traditional seismic analysis is inadequate to accurately predict a tsunami’s energy. Recently, two independent approaches have been proposed to determine tsunami source energy: one inverted from the Deep-ocean Assessment and Reporting of Tsunamis (DART) data during the tsunami propagation, and the other derived from the land-based coastal global positioning system (GPS) during tsunami generation. Here, we focus on assessing these two approaches with data from the March 11, 2011 Japanese tsunami. While the GPS approach takes into consideration the dynamic earthquake process, the DART inversion approach provides the actual tsunami energy estimation of the propagating tsunami waves; both approaches lead to consistent energy scales for previously studied tsunamis. Encouraged by these promising results, we examined a real-time approach to determine tsunami source energy by combining these two methods: first, determine the tsunami source from the globally expanding GPS network immediately after an earthquake for near-field early warnings; and then to refine the tsunami energy estimate from nearby DART measurements for improving forecast accuracy and early cancelations. The combination of these two real-time networks may offer an appealing opportunity for: early determination of the tsunami threat for the purpose of saving more lives, and early cancelation of tsunami warnings to avoid unnecessary false alarms.

Appendices

Appendix 1: GPS Method

The Global Navigation Satellite System (GNSS) (including American GPS, Russian GLONASS, Chinese Beidou, European Galileo, Japanese QZSS, and Indian IRNSS) will soon emerge as a practical technology for monitoring, forecasting and rescue management and future planning. By 2020, there will be over 160 GNSS satellites broadcasting over 400 signals across the L-band, nearly double the number today at any location. Currently, there are increasing numbers of dense GNSS networks with hundreds of stations that can and do provide real-time data in coastal regions that have a history of great tsunamigenic earthquakes, including the Cascadia subduction zone, the Japanese archipelago, and the Southern American trench. This represents an untapped new valuable data source for tsunami and earthquake early warning.

NASA has funded the early development of a prototype real-time GPS-based rapid earthquake and tsunami characterization system (http://www.nasa.gov/topics/earth/features/tsunami-20080117.html). Real-time GPS data, in various forms and processing stages, were primarily from the NASA Global Differential GPS System (GDPS—http://www.gdgps.net). With 150+  real-time GPS receivers, it is the world’s largest globally distributed, centrally managed real-time GPS network. The GDGPS system is a high accuracy GPS augmentation system, developed by JPL, to support the real-time positioning, timing, and orbit determination requirements of NASA science missions. Its real-time products are also used for GPS situational assessment, natural hazard monitoring, emergency geolocation, and other civil and defense applications.

The utility of GPS data for earthquake magnitude determination and for tsunami monitoring was demonstrated by Blewitt et al. (2006) who showed how the correct magnitude of the 2004 Sumatra earthquake could have been determined in near real time by complementing the seismic data with measurements from relatively few GPS tracking sites. Song (2007) demonstrated that coastal GPS stations are able to detect continental slope displacements of faulting due to large earthquakes and the GPS-detected displacements are able to estimate the disturbed oceanic energy and tsunami scales within a few minutes after the quake. This innovative method has been successfully tested in a NASA pilot program—the GPS-aided Real-Time Earthquake and Tsunami (GREAT) Alert System, during the events of the 2010 Chilean M 8.8 earthquake and associated tsunami (NASA release 2010: http://www.nasa.gov/topics/earth/features/tsunami_prediction.html; Naranjo 2013). The basic steps are the following: (1) locate an earthquake epicenter from seismometers (a few minutes after an initial earthquake); (2) collect near-field GPS-derived land velocities and infer the seafloor motions (a few more minutes of latency are possible); and (3) calculate the tsunami energy based on the GPS-predicted seafloor motions and local topography. For locations with real-time GPS coastal station and an earthquake model that includes local topography, the time to estimate the tsunami energy could be as little as 5 min after earthquake origen time. This approach was further tested after the 2011 Japan tsunami (Song et al. 2012; Xu and Song 2013).

For completeness, here we briefly describe the GPS approach. Notice that GPS stations only measure near-source ground motions of faulting, not the undersea ground motions that generate tsunamis. Therefore, the seafloor displacements due to the earthquake have to be derived from near-source GPS measurements. Although sophisticated GPS-inversion models have been explored using the shear and tensile fault theory in a half-space (Okada 1985), for tsunami prediction, we only need the seafloor (surface) displacements. Because great earthquakes cause large-scale ground motions, the coastal GPS-measured displacements are part of the major motions at the fault. Therefore, the seafloor displacements can be projected or extrapolated from the near-source GPS measurements. For tsunami early warnings, time is of the essence. To get the tsunami source from the real-time GPS measurements without delay, Song (2007) formulated an empirical profile extrapolation to project the seafloor displacements.

The formulation is as follows. Let (∆E _j , ∆N _j ) be the near-field GPS offsets (horizontal displacements) in the eastward and northward direction, respectively, where ∆ represents a small increment. The subscript j labels the GPS stations alongshore (parallel to the fault line), but only for those closer stations. For each j station, a cross-shelf section through the station, usually perpendicular to the fault, is identified. The seafloor displacements along the cross-shelf section are modeled to follow an empirical profile until the frontal edge of the fault (the trench):

$$\Delta E(r) = \Delta E_{j} \exp (r^{2}_{j} - r^{2} ) + \Delta e_{j2}$$

(1)

$$\Delta N(r) = \Delta N_{j} \exp (r^{2}_{j} - r^{2} ) + \Delta n_{j2}$$

(2)

$$\Delta U(r) = \alpha \sqrt {\Delta E(0)^{2} + \Delta N(0)^{2} }\times \left\{ {\exp ( - ar^{2} ) - \sqrt {\frac{\pi }{4a}} \exp ( - r)} \right\}$$

(3)

Here, r = d/W is the normalized distance such that r = 0 on the fault line and r = r _j at the j-th GPS station and, d is the physical distance from the fault line perpendicularly, and W is the mean distance of the GPS stations from the fault. An empirical value of W = 320 km is used for the 2004 Sumatra and 1964 Alaska earthquakes and 250 km for the 2011 Tohoku earthquake. If more GPS stations are available (besides the j labeled stations), a least-square fitting between the model results and the data is used for (∆e _j2, ∆n _j2) to correct the projection, where the j2 represents the extra GPS stations. The vertical uplift/subsidence is determined from the horizontal displacements by conserving mass such that the uplift volume equals the subsidence volume, in which the non-dimensional constant a = 1.5 is used. Notice that r represents a dimension in the cross-shore direction, while j gives the other dimension in the along shore direction. Once the seafloor displacements are derived from the GPS measurements, they are then interpolated onto a grid over the fault area, which is defined by the source model.

As seawater is almost incompressible, the vertical displacement of seafloor in each increment would result in a similar motion on the sea-surface because water particles over the fault area cannot flow away before being pushed up. In addition, the horizontal displacements of slope also contribute vertically to the sea-surface deformation (Tanioka and Satake 1996). The total sea-surface perturbation due to the vertical motion (the vertical acceleration of water particles) in each increment is

$$\Delta \eta \approx \Delta h = \Delta U + \Delta E \times h_{x} + \Delta N \times h_{y} ,$$

(4)

where η is the sea-surface elevation from the mean sea level, h is the water depth, Δ represents a small increment, and h _x and h _y are the eastward and northward slopes of the subfault surface, respectively. The tsunami potential energy gained by the ocean in each increment is

$$\Delta PE = g \times \rho \times \frac{{\Delta \eta^{2} }}{2} \times \Delta x \times \Delta y,$$

(5)

where g is the gravity acceleration and ρ is the water density. The total accumulated potential energy is the integration of Eq. (5) over the whole faulting area and within the period of rupture time.

To determine the kinetic energy transferred from a moving slope, we need the displacement velocity of the seafloor. Since each increment of the displacements has been achieved within a time interval of Δt s, the three-dimensional motion (velocity) of a subfault can be obtained by (ΔE, ΔN, ΔU)/Δt. The vertical acceleration of water particles does not contribute to the tsunami propagation; however, the resultant sea-surface perturbation would give the ocean potential energy, which has been accounted for in Eq. (5). The accelerated horizontal velocity of water particles in the vicinity of moving slopes can be written as:

$$\Delta u_{b} \left( z \right) = \left\{ {\begin{array}{ll} {\Delta E/\Delta t,} & { - h \le z \le - ~~R_{x} = h - \min \left\{ {h,~L\left| {h_{x} } \right|} \right\}} \\ {0,} & {{\text{otherwise}}} \\ \end{array} } \right.$$

(6)

$$\Delta v_{b} \left( z \right) = \left\{ {\begin{array}{ll} {\Delta N/\Delta t,} & { - h \le z \le - ~~R_{y} = h - \min \left\{ {h,~L\left| {h_{y} } \right|} \right\}} \\ {0,} & {{\text{otherwise}}} \\ \end{array} } \right.$$

(7)

Here, t is the rise-time of the subfault, z is the vertical coordinate, and L _H is the effective range of the horizontal motion (Song et al. 2008). The tsunami kinetic energy gained by the ocean in each second due to the horizontal motion is

$$\Delta {\text{KE}} = \rho \times \frac{1}{2}\left( {\Delta u_{b}^{2} + \Delta v_{b}^{2} } \right) \times \Delta z \times \Delta x \times \Delta y$$

(8)

where Δz is the vertical grid size in the bottom layer. The total accumulated kinetic energy is the integration of Eq. (8) over the whole faulting area and within the period of rupture time. The sum of integrating Eqs. (5) and (8) gives the total tsunami source energy (E _T):

$$E_{\text{T}} = \mathop \sum \limits_{i,j} \Delta {\text{PE}} + \mathop \sum \limits_{i,j,k} \Delta {\text{KE}}$$

(9)

where the i, j, and k indicate the horizontal and vertical grid, respectively.

Appendix 2: DART Method

The National Oceanic and Atmospheric Administration (NOAA) has developed a tsunami forecast system assimilating real-time tsunami measurements from Deep-ocean Assessment and Reporting of Tsunami (DART) stations (http://www.ndbc.noaa.gov/dart.shtml/) for more accurate tsunami warnings (Titov et al. 2005; Titov 2009; Bernard and Titov 2015). The forecast system integrates three key components, the DART data, the pre-computed tsunami propagation database, and high-resolution, site-specific flooding forecast models. They correspond to the three distinct stages of tsunami evaluation: generation, propagation and coastal run-up. Together, these components provide a full forecast capability. Presently, the DART network is composed of 63 stations, and is globally owned by nine countries (Fig. 1b). The propagation database has 1990 unit scenarios covering the most active subduction zones in the Pacific, Indian and Atlantic Oceans. High-resolution flooding forecast models have been developed for 75 US coastal communities (54 in the Pacific, 21 in the Atlantic/Caribbean).

The tsunami forecast system has been tested with 40 past tsunamis and has shown good accuracy when compared with observed tsunami amplitude time series at tide gauges, particularly for destructive tsunami waves (Titov et al. 2005; Titov 2009; Wei et al. 2008, 2012; Tang et al. 2008, 2009, 2012, 2015; Bernard et al. 2014). The NOAA system is designed to predict tsunami flooding, but not all tsunamis produce flooding. However, during the 2011 Japanese tsunami the NOAA forecast system correctly predicted flooding at Kahalui, Hawaii 4 h before the tsunami arrived. Another feature of the NOAA Center for Tsunami Research (NCTR) forecast system is to quickly estimate the tsunami energy. The total energy transmitted by tsunami waves is one of the most fundamental macroscopic quantities for quickly estimating the potential impact of a tsunami. As shown in Fig. 1b, DART stations are strategically positioned close to subduction zones throughout the world’s oceans, where they detect tsunamis within 1 h of most tsunamis generated. In 2012, Japan installed three additional DART stations (Fig. 1) close to the Japanese coastline. With these additional stations, the December 7, 2012 Japanese tsunami energy could have been estimated within 20 min after the earthquake origen time (Bernard et al. 2014).

The DART inversion approach exploits the simple linearity assumption of tsunami propagation to circumvent the large uncertainties in earthquake sources. Wave dynamics of tsunami propagation in deep-ocean are assumed to be linear (Mei et al. 2005; Kanoglu and Synolakis 2006; Liu 2009). Therefore, a propagation scenario (η) can be decomposed into a superposition of a subset of source units (h):

$$h = \mathop \sum \limits_{i = 1}^{N} \alpha_{i} h_{i}$$

(10)

h is the tsunami amplitude, h _i is the tsunami amplitude of the i-th source unit, α _i is the i-th coefficient, N is the total number of source units used. A source unit h _i is a unit propagation scenario that could be generated by an earthquake, a landslide or other type of sources. Currently, each source unit h _i in the forecast system represents a unit propagation scenario generated by a typical Mw 7.5 subduction earthquake. It is modeled as the instantaneous rupture of a single rectangular fault plane characterized by pre-defined parameters described the location, orientation and rupture of the plane (Gusiakov 1978; Okada 1985). Each source unit has predefined parameters of length 100 km, width 50 km, slip 1 m and rake angle 90°. Other parameters, including depth, dip angle, and strike angle, are location-specific (Titov et al. 1999; Gica et al. 2008).

The DART-inverted tsunami source refers to the source scenario or the synthetic propagation scenario (η) from a superposition of source units (h) that gives a best fit to the observed tsunami amplitude time series (η _obs) in deep ocean:

$$\eta_{\text{obs}} = \eta + {\text{res}}_{t} = \mathop \sum \limits_{i = 1}^{N} \alpha_{i} h_{i} + {\text{res}}_{t}$$

(11)

$$\mathop \sum \nolimits {\text{res}}_{t}^{2} = \mathop \sum \nolimits \left( {\eta_{\text{obs}} - \eta } \right)^{2} = \mathop \sum \nolimits \left( {\eta_{\text{obs}} - \mathop \sum \limits_{i = 1}^{N} \alpha_{i} h_{i} } \right)^{2}$$

(12)

η _obs is the observed tsunami amplitude in deep ocean, res_t is a residual term between the observation and model. α _i was estimated using the least squares method (Eq. 12), e.g., minimizing the sum of the squared residuals (Percival et al. 2011).

The Method Of Splitting Tsunami (MOST) model, a suite of finite difference numerical codes based on the long-wave approximations (Titov and Synolakis 1998; Titov and Gonzalez 1997; Synolakis et al. 2008), is used to compute the tsunami propagation caused by source units:

$$\frac{\partial h}{\partial t} + \frac{\partial }{\partial x}\left[ {\left( {h + d} \right)} \right]u + \frac{\partial }{\partial y}\left[ {\left( {h + d} \right)} \right]v = 0$$

(13)

$$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = - g\frac{\partial h}{\partial x} - \frac{{C_{f} u\left| u \right|}}{h + d}$$

(14)

$$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = - g\frac{\partial h}{\partial y} - \frac{{C_{f} v\left| u \right|}}{h + d}$$

(15)

where d is the undisturbed water depth, u = {u, v}, u and v are the depth-averaged velocities in the x and y directions, respectively, g is the gravity acceleration, f is the Coriolis parameter, C _f  = gn ²/(h + d)^1/3, and n is Manning coefficient (assumed to be small during deep-ocean propagation). In the MOST model, these equations are solved numerically using the splitting method (Titov and Synolakis 1998).

At a given time, the instantaneous tsunami wave energy E is the sum of the kinematic energy, KE, and the potential energy PE. It is computed as:

$$E = {\text{KE}} + {\text{PE}}$$

(16)

where

$${\text{KE}} = 0.5\rho \mathop \int \limits_{s}^{{}} \left( {\eta + d} \right)\left( {u^{2} + v^{2} } \right){\text{d}}s$$

$${\text{PE}} = 0.5\rho g\mathop \int \limits_{s}^{{}} \eta^{2} {\text{d}}s$$

ρ is the ocean water density and g is the gravity acceleration coefficient. The integration is over the propagation domain S.

Based on the zero initial velocities assumptions according to the adopted elastic deformation model, the tsunami energy, E _T (the total energy transmitted by the tsunami waves), is then computed as the difference between the potential energy from the initial deformed body of water and that of the quiescent state in the same basin:

$$E_{\text{T}} = {\text{KE}}_{t = 0} + {\text{PE}}_{t = 0} = 0 + 0.5\rho g\int\limits_{s} {\eta_{0}^{2} {\text{d}}s}$$

(17)

where η ₀ is the sea-surface elevation at the initial stage (e.g., the initial ocean surface deformation). Substituting Eq. (10) into Eq. (17), we have the numerical formula:

$$E_{\text{T}} = 0.5\rho g\mathop \sum \limits_{j}^{{}} s_{j} \left( {\mathop \sum \limits_{i = 1}^{N} \alpha_{i} h_{0,i,j} } \right)^{2}$$

(18)

where j represents j-th wet grid node, s _j is a cell area corresponding to the j-th node, h _0,i,j is the initial sea-surface displacement from the i-th tsunami source function at j-th node. More discussion on tsunami energy can be found in Tang et al. (2012).