Instant tsunami early warning based on real-time GPS – Tohoku 2011 case study

Taking the 2011 Tohoku earthquake as an example, we demonstrate the ability of real-time GPS to provide qualified tsunami early warning within minutes. While in earlier studies we demonstrated the power of the so-called GPS shield concept based on synthetic data, we here present a complete processing chain starting from actual GPS raw data and fully simulate the situation as it would be in a warning center. The procedure includes processing of GPS observations with predicted high precision orbits, inversion for slip and computation of the tsunami propagation and coastal warning levels. We show that in case of the Tohoku earthquake, it would be feasible to provide accurate tsunami warning as soon as 3 min after the beginning of the earthquake.


The Tohoku earthquake and tsunami
The magnitude M W = 9.0 Tohoku earthquake on 11 March 2011, with epicenter at 142.9 • E 38.1 • N (JMA) was the largest event ever instrumentally observed in Japan, and the fourth-largest world wide.In the last century, 22 earthquakes with magnitude between 7.5 and 8.3 occurred in that region (USGS, 2012), but the last event of comparable size was probably the Jogan Sanriku earthquake in 869, which has been assigned an upper magnitude estimate of 8.9 according to associated tsunami deposits reaching as far as 4 km inland.Two additionally identified sediment layers indicate a possible recurrence interval of about 1000 yr for these giant events (Minoura et al., 2001).
The earthquake-generated tsunami, which affected the Sendai, Miyagi and Iwate prefectures within less than half an hour, demonstrated run-ups exceeding 40 m (Mori et al., 2012) and locally inundated up to 5 km inland (Fujii et al., 2011).It caused more than 15 000 fatalities and massive devastation, including damage of the Fukushima nuclear power plant.
Due to numerous and various sensor networks deployed in Japan and in the Pacific, this event is now the most extensively recorded megathrust event ever.Coseismic deformation was recorded by the GEONET GPS array (GSI, 2011) and deduced by sea floor geodesy (Sato et al., 2011).Source parameters were the scope of many investigations including inversion of ocean data, seismic waveforms, as well as joint inversions (Romano et al., 2012;Yokota et al., 2011).Most studies reveal a compact slip distribution (Pollitz et al., 2011), very high maximum slip and average slip about twice as high as for the M W = 9.3 Sumatra earthquake (Tajima et al., 2013).

GPS in tsunami early warning
Near-field tsunami early warning systems (TEWS) for coastal regions exposed to subduction earthquakes should be able to provide a warning as early as 5 to 10 min after an earthquake.Working within this time limit, traditional seismic methods tend to underestimate the total moment magnitude in case of giant (Sumatra 2004, Tohoku 2011) as well as slow, "tsunami" earthquakes (Java 2006, Mentawai 2010).
In case of the Tohoku event, the first tsunami warning bulletin was issued by the Japan Meteorological Agency (JMA) already 3 min after origin time and was based on the M JMA = 7.9 magnitude estimate of the earthquake early warning (EEW) system.Based on this earthquake magnitude and location, a major tsunami warning was declared for the three nearby prefectures.While being rated as "major", the warning still significantly underestimated the true coastal tsunami impact due to large magnitude underestimation: M JMA saturates for events with M W >8.0 (Hoshiba et al., 2011).The JMA magnitude was revised to M W = 8.4 after one hour and twenty minutes, and then, after 3 h, further upgraded to M W = 8.8.Only after several days, the earthquake was found to have M W = 9.0 (Ohta et al., 2012;Ozaki, 2011).Thus, the earthquake was dramatically underestimated during the period of highest tsunami activity.In fact, real seismic moment results to be more than 30 times larger than reported in the beginning.
The Great Sumatra earthquake and tsunami in 2004 triggered active development and operational implementation of seismic methods capable of providing full moment magnitude in much shorter time than traditional methods.For example, the W-phase moment tensor inversion method was demonstrated to be able to provide a moment solution with M W = 9.0 as soon as 20 min after origin time (Duputel et al., 2011).Despite this progress, approximation of a tsunami source as a point source, i.e. epicenter plus magnitude, is usually insufficient for reliable near-field tsunami early warning, since for large subduction events happening not far from the coast, rupture propagation direction and source extent cannot be neglected (Babeyko et al., 2010).Buoys together with ocean bottom sensors help to gain additional source information in case of transoceanic tsunamis, but would not be effective in the near-field due to time constraints.
Recent progress in high-precision real-time GPS processing makes GPS arrays a valuable component of near-field TEWS (Blewitt et al., 2006(Blewitt et al., , 2009;;Ohta et al., 2012;Sobolev et al., 2006Sobolev et al., , 2007)).GPS receivers placed close to the epicenter minimize response time while not suffering from clipping as broadband seismometers do, or tilting artifacts as accelerometers after double integration for displacement.They directly measure surface deformation, which is crucial for tsunami source inversion.Started after the event of Sumatra 2004, the concept of the "GPS shield" was proposed (Sobolev et al., 2006(Sobolev et al., , 2007)), where we suggested incorporating GPS-based tsunami early warning in all tsunami-prone areas worldwide, including Japan, and which is now becoming an operative part of the German Indonesian Tsunami Early Warning System GITEWS (Babeyko et al., 2010;Behrens et al., 2010;Falck et al., 2010).
The Geospatial Information Authority of Japan operates the unique GEONET array, comprising more than 1200 realtime GPS stations.
While not yet integrated into the tsunami early warning procedure, this network has great potential to contribute.We demonstrate this by replaying the Tohoku 2011 event.To do so, we combine real-time GPS Precise Point Positioning (PPP) processing software developed at the German Research Centre for Geosciences in Potsdam (Ge et al., 2009), with fast slip inversion and tsunami computation.This allows observing of the earthquake growth in almost real-time and simultaneously computing tsunami warning levels.
It is important to note that our primary goal in this study is not to provide the most accurate source model, but to show that even a fast inversion of real-time quality data from a limited number of GPS stations provides a source model which is adequate for tsunami early warning purposes.

Geodetic data
We use GPS raw data from the GEONET (2011) array operated by the Geospatial Information Authority of Japan from one day before to one day after the earthquake, with a sampling interval of 30 s. Data of the previous day are processed in static mode for obtaining precise positions serving as reference for the displacement.Of the more than 1200 stations we select a subset, as illustrated in the upper left panel of Fig. 1: the orange dots are on epicentral rays up to 550 km distance at increasing spacing.The first two GPS stations on each ray are selected.This procedure results in 50 stations located close to the coast around the earthquake and is motivated by the findings of Sobolev et al. (2007), which show that a set of two stations perpendicular to the coast line, and thus capturing the gradient of displacement, is much better in locating slip in depth than a single station, and the spacing along the coast is sufficient to map the slip distribution along the trench.Performing the procedure using only 20 stations, or using all of the stations in Japan does not significantly change the slip model.

GPS real-time processing
We employ the PPP (Precise Point Positioning) approach (Kouba, 2003;Larson et al., 2003;Shi et al., 2010), which includes two parts: server-end and client-end.At the serverend, precise satellite clocks and orbits are estimated using data from a global reference network and made available for the monitoring stations, for example the International GNSS services (IGS) real-time products, so that the coordinates of the client-end stations can be estimated station-by-station in the reference frame defined by the orbits and clocks.Usually, dual-frequency geodetic receivers are used in order to eliminate ionospheric delays, while receiver clock bias, tropospheric delays and carrier-phase ambiguities are estimated together with station positions.Ambiguity-fixing can be performed at a single station, which improves the horizontal accuracy significantly.However, orbital and clock errors are propagated directly into the estimated position, as they are fixed in the estimation.The station-by-station processing is very suitable, especially for large-scale dense GNSS networks as the GEONET array with several hundreds of realtime stations.
GFZ has been working on real-time GNSS data processing software since 2007 for various applications in geosciences.The EPOS-RT software (Ge et al., 2009) is now providing real-time global PPP service with ambiguities fixing of 2-3 cm-level accuracy (Ge et al., 2011).The results are shown online at http://kg6-dmz.gfz-potsdam.de/rtgnss/.
In this study we use exactly the same processing strategy to simulate the real-time PPP service to assess the capability of such a service in a tsunami early warning system.We use data from about 90 IGS reference stations to estimate satellite clocks epoch-by-epoch with the predicted part of the IGS Ultra-Rapid orbits fixed.With the estimated clocks and the fixed orbits kept in a file, PPP is carried out in simulated real-time kinematic mode with an inter-epoch constraint on station coordinates of 0.05 m/ √ s.Therefore the whole data processing is a replay of the real-time PPP service running at GFZ now.In order to get an accuracy of 1 cm in horizontal components and 2-4 cm in vertical (Geng et al., 2010), a convergence phase of 1 to 2 h is needed, which is in general not a problem for the continuously operational monitoring stations.
The lower panels of Fig. 1 show obtained displacement time series with 30 s sampling interval for the four stations depicted in the upper left panel.Stations close to the source (middle left panel) show ramp-like behavior and large static displacements.Far away stations (lower right panel) show more wave-like signals and small static displacements.

Slip inversion
The inversion procedure follows a straightforward static approach.An inversion of the 3 component displacements of the selected stations is performed every 30 s using Green's functions for elastic deformation of a layered medium computed using code PSGRN/PSCMP (Wang et al., 2006).Fault geometry consists of 216 rectangular subfaults on a curved surface with varying dip angle based on RUM slab geometry (Gudmundsson and Sambridge, 1998).The resulting slip coefficients minimize a cost function consisting of chi-squarereduced (misfit), smoothing terms for slip magnitude and slip azimuth including boundary conditions.Calibration of the smoothing weight factors is done in advance using synthetic forward models comprising checkerboard, homogenous and realistic slip distributions by adding noise to the resulting displacement vectors and matching the inverted to the forward models.The whole procedure is described in more detail by Hoechner et al. (2008).Resolution analyses for several hypothetical GPS sites constellations and assumed GPS realtime accuracies for Sumatra and Japan were presented earlier (Babeyko and Hoechner, 2012).

Tsunami modeling
The coseismic vertical sea floor deformation associated with the different slip models is used as initial condition for tsunami modeling.Tsunami wave propagation is computed using our code "easyWave" on a 2 arc minute grid resampled from ETOPO1 global bathymetry (Amante and Eakins, 2009).It integrates shallow-water equations in linear approximation on spherical coordinates.The numerical leapfrog scheme on a staggered finite-difference grid essentially follows the well-established and widely used TUNAMI-F1 algorithm (Imamura et al., 1997).Boundary conditions impose full reflection on land and free transmission at the open sea boundary, and the computational time step is 2 s, satisfying the Courant stability criterium.Since inundation modeling is not performed directly, run-ups are extrapolated from offshore positions (typically, in 50-100 m water depth) using Green's law accounting for wave shoaling, an approach which is also employed by the Japanese tsunami early warning system (Kamigaichi, 2009).
Organizing the computation on an expanding grid makes it very efficient, especially at the early stages of tsunami propagation.Even using a serial version of easyWave it is possible to estimate tsunami impact at the nearby Japanese coast within 10-15 s (for 2 h model time).This actually eliminates the need to use pre-computed tsunami Green's functions as done in other studies, e.g.Ohta et al. (2012).Computations for much larger distances and thus longer model times (like for the whole Pacific, 24 h model time) require, of course, longer simulation runs with the serial version of easyWave.However, our first experiences with an adapted version for GPU (Graphical Processing Units), show that it is still possible to provide trans-oceanic simulations in a few minutes.Anyway, timing is not critical for long-distance forecasting.Real-time GPS displacements (left panels) are inverted for slip at the fault (middle panels).These source models are used to compute estimated coastal tsunami wave heights (right panels).Lower right panel: observed wave heights.The table contains key characteristics of the three succeeding earthquake and tsunami models.M w : moment magnitude, H max : maximum observed horizontal displacement at the GPS stations, S max : maximum inverted slip at the fault, U max : maximum modeled uplift at the sea floor, V t : initial tsunami volume, E t : initial tsunami potential energy, W max : maximum estimated wave height at the coastal points of interest.

Results
Figure 2 shows real-time horizontal displacements at the GPS stations, inverted source models and predicted coastal wave heights for 60, 90 and 180 s after earthquake origin time.Here, time denotes the lower limit imposed by physics, that is rupture propagation and seismic wave travel time.Processing time is discussed in Sect. 5.The coloring scheme roughly follows the JMA thresholds, additionally showing in magenta peak coastal amplitudes exceeding 10 m.One minute after origin time, inverted maximum slip is around 6 m (see table in Fig. 2).Nevertheless, the estimated magnitude already corresponds to M W = 8.4 and warning levels show major tsunami warning for the Sendai prefecture.Thirty seconds later, the estimate reaches M W = 8.7.Large coastal areas are already expected to have more than 10 m peak wave height.The slip model stabilizes about three minutes after origin time.Magnitude reaches its final value of M W = 9.0.Peak wave height prediction is 29 m.The lowermost panel of Fig. 2 contains field measurements of coastal wave heights (Mori et al., 2012) mapped to our points of interest and color scheme.Comparison of the warning forecast made after 3 min shows the high quality of GPS-based instant forecasting for the Tohoku event.
Predicted mareograms at 3 min origin time are compared to later buoy observations in the right panels of Fig. 3.For the nearest buoy (D21418), the model results in some amplitude and high-frequency deficit.This is probably because of inadequacy to fully resolve a high-slip area close to the trench (Romano et al., 2012), but additionally a submarine landslide may have been involved.Far-distant buoys (D21413 and D32412) show very good correlation.This is because for large distances and such a compact event, the tsunami initial condition acts almost like a point source and details of the slip distribution have only a minor effect.
In general, run-up and mareogram fit is good, total moment is well resolved, and accurate tsunami early warning is possible for nearby coastal areas as well as concerning transoceanic regions.

Discussion
GPS monitoring in regions close to subduction zones has been operational in several places worldwide for many years (e.g.Japan, Indonesia, Chile, Cascadia).Besides being used to derive slip and afterslip, obtained data is also used for establishing locking conditions at the plate interface (Loveless and Meade, 2010) or to derive rheological properties (Hoechner et al., 2011).However, to our knowledge there is not a single example for successful tsunami early warning based on GPS up to now.We stress the importance that existing stations be integrated into tsunami early warning systems, and that additional stations be set up at strategic positions specifically suitable for that purpose.
Proximity to the source, enabling GPS to deliver useful data within a short time, is also a requirement, since static displacement decays much faster with distance than seismic waves do (Fig. 1, upper right panel).Maximum useful distance is given by earthquake magnitude, which imposes a lower limit of earthquake size detectability on a specific GPS array.Stations close to the source, e.g. on islands between the trench and the mainland, like in the case of Sumatra, allow resolution of smaller earthquakes.However, large tsunamigenic subduction earthquakes can always be identified using mainland coastal GPS stations (Sobolev et al., 2007).In case of the Tohoku earthquake and tsunami, an accurate tsunami warning based on GPS could have been issued about three minutes after beginning of the earthquake.This lower limit is given by the time it takes to complete the rupturing process and the time it takes for the seismic waves to reach the receivers.As shown in Fig. 1, the full static displacement signal is then established at the nearest stations.Additional time is required to locally process GPS raw data using high precision satellite orbit estimates and to transmit displacement data to the warning center (10 s), to invert for the source model (30 s or less, depending on technique), and to compute tsunami propagation (10 s for nearby coastal regions), all in all less than 1 minute.Even though this approach is not yet implemented in an operational manner, we are confident that these timing conditions can also be obtained in an actual warning center.
The study by Ohta et al. (2012), which was published while the present paper was in preparation, reaches similar conclusions, even though there are significant methodological differences (long baseline real-time kinematic GPS (RTK-GPS) instead of PPP, inversion for a single fault plane including geometry instead of high-resolution discretization based on subduction interface geometry, homogeneous Earth model (Okada) instead of layered half-space, usage of tsunami wave form Green's functions instead of on-the-fly computation).Their proposed automatic algorithm applied to the Tohoku event produces a final fault model of magnitude 8.7 275 s after origin time.However, this magnitude is obtained also in their study already after 170 s.The reason why it is not reaching the full value of 9.0 may be the single fault approach, which does not allow for realistic slip distribution or application of the homogeneous half-space Earth model.
The Sumatra 2004 earthquake had a much larger spatial extent than the Tohoku earthquake of about 1600 km along the trench and a duration of about 10 min (Krüger and Ohrnberger, 2005).Of course, final magnitude is only recognizable after completion of the rupturing process.However, local rise time is much shorter than overall rupturing time, so once the rupture front has passed, nearby stations could provide local early warning information within a similar amount of time as for the Tohoku earthquake.
In Japan, the distance between coast and trench is about 200 km.This is actually too large to be able to fully resolve slip in the near-trench areas (Sobolev et al., 2007).Vertical displacement information from directly above the source could be provided by GPS buoys, but these are much more expensive in set up and maintenance than land stations.Geodetic data can also be obtained directly from the sea floor by a combination of acoustic and GNSS techniques (Sato et al., 2011), however this can only be done by extensive ship cruises after the event.Anyway, as noted above, the aim of this study is not to get the most accurate source model, but one which is adequate for tsunami early warning as quickly as possible, which is demonstrated in Fig. 2. Loveless and Meade (2011) compared various slip models of the Tohoku earthquake with their study of locking conditions at the Japan subduction zone.The correlation of the region with strongest coupling and maximum slip is striking.For a better understanding of the characteristics of the seismic cycle and governing processes however, much longer observation times are required.These can partially be obtained with paleogeodetic methods (Sieh et al., 2008).Nevertheless, the above mentioned study and others still hold several "locking hot spots".This should be motivation enough to implement permanent real-time GPS networks along all subduction zones and to provide data quality and access on a level as the well-established seismic networks have been providing for years.

Fig. 1 .
Fig. 1.Upper left panel: grey dots: GEONET stations; magenta: selected stations for real-time inversion; orange: epicentral rays used to select the subset of stations; circles: stations shown in the lower panels.Upper right panel: coseismic displacement as a function of epicentral distance for the selected stations.Lower panels: real-time displacement solutions at 30 s sampling interval.The yellow bars indicate the time points for which inversions are shown in Fig. 2.
et al.: Instant Tsunami early warning based on real-time GPS Fig. 2.Real-time GPS displacements (left panels) are inverted for slip at the fault (middle panels).These source models are used to compute estimated coastal tsunami wave heights (right panels).Lower right panel: observed wave heights.The table contains key characteristics of the three succeeding earthquake and tsunami models.M w : moment magnitude, H max : maximum observed horizontal displacement at the GPS stations, S max : maximum inverted slip at the fault, U max : maximum modeled uplift at the sea floor, V t : initial tsunami volume, E t : initial tsunami potential energy, W max : maximum estimated wave height at the coastal points of interest.

Fig. 3 .
Fig. 3. Left panel: Maximum wave heights in the Pacific predicted by the model inverted after 180 s.Thin black lines are isochrones in 1 h intervals.White dots denote positions of the buoys for which model and observation are shown in the right panels.
et al.: Instant Tsunami early warning based on real-time GPS