For a PDF version
of this article >> Click Here (approx. 4.4M)
|Journal of Oceanography, Vol. 62, pp. 705 to 716, 2006
HF Radar Detection of Tsunamis
BELINDA J. LIPA*, DONALD E. BARRICK, JOHN BOURG and BRUCE B. NYDEN
CODAR Ocean Sensors, 125 La Sandra Way, Portola Valley, CA 94028, U.S.A.
(Received 4 January 2006; in revised form 4 May 2006; accepted 2 June 2006)
This paper demonstrates that HF radar systems can be used to detect tsunamis well before their arrival at a coastline. We solve the equations of motion and continuity on the ocean surface using models to simulate the signals produced by a tsunami approaching the east U.S. coast. Height and velocity profiles are derived along with expressions for the radar-observed current velocities in terms of bathymetry and tsunami height and period. Simulated tsunami-generated radial current velocities are superimposed on typical maps of radial velocity generated by a Rutgers University HF radar system. A detection parameter is defined and plotted to quantify the progress of the tsunami, which is shown to be detectable well before its arrival at the coast. We describe observations/warnings that would have been provided by HF radar systems at locations in the path of the 2004 Indian Ocean tsunami.
Keywords: Remote sensing, modeling, HF radar oceanography
Existing tsunami watch systems are based on computer modeling programs that warn against the possibility of earthquake-generated tsunami impacts, and attempt to predict their strength and location. These computer models include ocean scale bathymetry and could run continually for possible epicenter points, with the input parameters updated in the event of a potentially dangerous earthquake. A protocol for the rapid dissemination of seismic data and tsunami model predictions between foreign governments is being established, but no system as yet exists for local detection of an actual incoming wave with a significant warning capability.
Barrick (1979) originally proposed the use of shorebased HF radar systems for tsunami warning. HF radar systems presently operate continuously from many coastal locations around the globe, monitoring ocean surface currents and waves to distances of up to 200 km. Tsunami watch software could run in the background, activating a warning should a tsunami be detected. This information would be available to local authorities and would be invaluable if international communications fail or are too general in their predictions. Global models may be inadequate for localized areas for which the bathymetry is inaccurate; for example, many developing nations do not conduct regular bathymetric/hydrological surveys. In addition, when a quake epicenter is close to shore there may be insufficient time for the international communication
* Corresponding author. E-mail: firstname.lastname@example.org
Copyright© The Oceanographic Society of Japan/TERRAPUB/Springer
||chain to be activated; in such cases local systems would
provide the only advance warning. It would alleviate the
false-alarm problems that plague existing tsunami watch
systems. Finally, computer rediction models and early
warning schemes apply only to tsunamis generated by
earthquakes; HF radars would also be able to detect tsunamis
due to underwater rockslides and tidal bores.
This paper demonstrates that modern HF radar systems
already in use for observations of ocean currents
and waves are capable of detecting tsunamis well before
their arrival at many coastlines. The warning time increases
with the width of the continental shelf. This study
was inspired by the disastrous tsunami generated by an
earthquake off Sumatra on December 26, 2004.
2. Theoretical Analysis
In the open ocean, a tsunami has a low height, an
extremely long wavelength and travels at great speed. It
can be detected in real time only by bottom-mounted pressure
sensors. However, a tsunami slows down as it moves
into shallower water, allowing it to be detected by shorebased
HF radar systems well before its arrival. Although
the wave period remains invariant, the height and the orbital
velocity increase and the wavelength decreases. The
Navier-Stokes equation and the equation of continuity,
see for example Lamb (1932), Kinsman (1965), form the
basis for tsunami modeling. In this section we present
two approaches to solving these equations. Subsection 2.1
reproduces the closed-form expressions for tsunami parameters
originally presented by Barrick (1979), which
are based on linear wave theory, assuming piece-wise
constant depths and a sinusoidal profile for the tsunami.
|In Subsection 2.2 the depth is allowed to vary and arbitrary
height profiles may be used; equations must then be
solved numerically. In Subsection 2.3 we give an expression
for water depth below which linear wave theory can
be expected to break down and nonlineartities will prevail.
2.1 Linear approximation
As the wavelength of a typical tsunami is far greater
than the water depth d, the defining equations can be simplified
by expanding them as Taylor series in d/λ(d),
where λ(d) is the wavelength, and retaining only the first
terms in the series. Barrick (1979) assumed the water
depth to be piecewise constant and, using a sinusoidal
waveform for the tsunami, derived the following equations
for tsunami properties. For shallow-water waves like
tsunamis, phase velocity and group velocity of the wave
are equal and given by:
where g is the acceleration due to gravity. The height of
the tsunami is expressed in terms of its value in deep
water, which was taken to have a depth of 4000 m:
where h4000, the height of the tsunami in deep water and
d are expressed in meters. The maximum surface orbital
or particle velocity for water of depth d is given by:
The relationship between the wave period P and wavelength
λ(d) is given by the shallow-water dispersion equation,
which defines how the wavelength changes with
The time for the tsunami to cover a distance L terminating
at the radar site is given in terms of the phase velocity
For a tsunami, the phase/group velocity of the wave
is much faster than orbital velocities, but from (1) it slows
down in shallow water as d1/2, increasing the warning
available from the time of detection. From (2) and (3),
||the height increases as d–1/4, while the maximum particle
velocity increases more quickly as d–3/4, increasing the
signal seen by the radar.
2.2 A more general model
The linear shallow-water wave methods of Barrick
(1979) and Subsection 2.1 are based on the assumption
of a piece-wise constant depth and a sinusoidal profile
representation of a tsunami. To validate these analytic results,
we now consider a more realistic general hydrodynamic
model, for which the depth is variable, using a
Gaussian-modulated sinusoid for the waveform.
The water depth is now taken to be a function of
horizontal position (x, y), and the height η above a mean
plane and velocity υ˜ are assumed to vary with position
and time. The following two equations form the basis for
tsunami modeling, e.g., Lamb (1932); Kinsman (1965):
Equation (6) represents Newton’s second law, involving
the expression of the horizontal pressure gradient and
acceleration in terms of the Navier-Stokes equation. Equation
(7) embodies the incompressibility of water, i.e., the
horizontal velocity must vary with position (x, y) as water
is confined between the two interfaces: bottom and
surface. It is assumed that the water velocity υ˜ is independent
of depth, which is a good approximation for the
time and space scales of tsunamis and the basis of most
numerical tsunami models (Mofjeld et al., 2000). Remaining
terms in the Navier-Stokes equations (e.g., Coriolis
and frictional effects) and differences between Lagrangian
and Eulerian time differentiation, are ignored as they are
unimportant in the tsunami scenario involving the propagation
of massive amounts of energy across large distances
on hourly time scales.
Because the tsunami height is small compared to the
depth except during the final run-up onto shore, the second
term on the left side of (7) can be ignored. Then (6)
and (7) can be combined to give two variations of the
Fig. 1. (a) Typical radial vector map measured by a long-range
SeaSonde located at Tuckerton, New Jersey, 3:30 pm January
27, 2005. Constant depth contours are labeled in meters.
(b) Schematic map with the wave fronts of an approaching
tsunami that are approximately parallel to the coast and
the depth contours. Velocity vectors are perpendicular to
the wave fronts, and the radial component observed by the
radar is defined by (15).
For constant depth d, (8) and (9) reduce to the classic
constant-depth plane-wave representations of tsunami
velocity and height.
As the tsunami moves into shallow water, refraction
||causes the water motion to be essentially normal to the mean large-scale isobaths. It is then valid to express the position-dependence of variables as a function of a single parameter x, which is defined as the distance from shore along a direction normal to the mean large-scale isobaths.
Reducing (9) to its one-dimensional form and using the method of separation of variables, it is shown in Appendix that solutions for arbitrary time dependence can be represented as:
is the solution to the following ordinary differential equation:
and p(ω) is the Fourier transform of the time-dependent waveform representing the tsunami in deep water as it approaches the continental slope; c is the group velocity in deep water (of depth D o) at the edge of the slope, which from (1) is given by:
To obtain the tsunami parameters, these equations must be solved numerically. The functional dependences of the tsunami parameters were found to be consistent with those for the simpler model, as given by Eqs. (1)– (5). It is therefore valid to use these analytic expressions to estimate the tsunami characteristics and available warning time. It is shown in Appendix that energy conservation principles also confirm these functional forms.
2.3 Limits on linear wave theory
When linear wave theory breaks down in shallow water, currents become stronger and hence they are more easily detected. In the nonlinear regime, wave energy is converted to massive water transport or surge. There are no simple equations to describe this final progress of the tsunami.
Barrick and Lipa (1986) expressed the surface height normalized by the wavelength as a perturbation expansion in h(d)/λ(d). For a solitary periodic wave traveling in the x-direction, when λ(d) >> d, the first- and secondorder terms for the surface height can be written as:
Fig. 2. Parameters plotted versus distance offshore from the Rutgers SeaSonde:
(a) Water depth, (b) tsunami wavelength,
(c) maximum orbital velocity, (d) wave height produced by the tsunami.
|Linear wave theory breaks down when the magnitude of the nonlinear term approaches that of the linear term. We assume that linear wave theory becomes invalid when the nonlinear term exceeds 50% of the linear term, which occurs at a depth given by:
For depths less than this threshold, higher-order wave nonlinearities become significant based on this definition of the linear domain. This depth is approximately four meters for typical deep-water tsunami amplitudes of 0.5 m. 3. Simulations A single HF radar system produces maps of the radial component of the ocean surface current velocity using radar echoes from short wind-driven waves. The tsunami produces surface velocities that superimpose on the slowly varying ambient current velocity background. They have a characteristic signature due to their coherence over large distances allowing them to be detected when they arrive in the radar coverage area.
We based our simulations on Eqs. (1)–(5) and on satellite measurements of the 2004 tsunami made by US and French satellites (Gower, 2005), taking values of 0.5 m for the deep-water wave amplitude and 600 km for the wavelength. From (4) the corresponding period of the wave is about 50 minutes, independent of depth.
||3.1 Tsunami signals in SeaSonde radial velocity maps Simulated tsunami current velocities were superimposed on current velocities measured by the Rutgers University SeaSonde system located at Tuckerton, New Jersey, which operates with a transmit frequency of 4.54 MHz. Voltage time series from the three SeaSonde antennas are Fourier-transformed to give complex frequency spectra. System parameters are set to provide spatial/temporal resolutions of 6 km/0.5 hours respectively and a radial current velocity resolution of 13 cm/s. The system measures the component of the current velocity that is radial to the radar to typical ranges of 150 to 200 km. Figure 1(a) shows the location of the radar, the bathymetry, and a typical map of measured radial current velocities.
As discussed by Barrick (1979), a sinusoidal tsunami wave appears as a periodic surface current. Its wave orbital velocity at the surface transports the much shorter waves seen by the radar, adding to the ambient current field and producing a signature detectable by the radar. We assume that the water depth over the radar coverage area can be adequately represented by depth contours parallel to shore. A tsunami will move perpendicular to these contours, as illustrated schematically in Fig. 1(b). The approximate depth profile is shown in Fig. 2(a). At time t and range x measured perpendicular to the coast from the radar, the surface current velocity produced by the tsunami can be written as:
|where the maximum surface orbital velocity and wavelength are given by (3) and (4) respectively. The component of this velocity radial to the radar is given by:
where θ is the angle between the depth contours and the line drawn to the radar, as shown in Fig. 1(b).
Figures 2(b)–(d) show the variations in wavelength, maximum orbital velocity and waveheight plotted as a function of distance perpendicular to the coast from the radar as the tsunami advances over the continental shelf, which is approximately 125 km wide in this region. It follows from (5) that approximately an hour and a half will elapse from the time the tsunami reaches the edge of the continental shelf before it strikes the coast.
Figure 3 shows the simulated radial velocity maps produced by the tsunami wave coming directly onshore, at half-hour intervals starting when the tsunami arrives at the outer edge of the radar coverage area and ending when it reaches the coast. These radial velocities would superimpose on the ambient current velocities, with radar map samples shown every half hour for the combined radial velocity field. The tsunami radial velocities shown in
Fig. 3 were added to the four measured radial velocity maps from 4 pm–5:30 pm, January 27, 2005. Including the measured radial map for 3:30 pm shown in Fig. 1, this gives a set of five radial maps starting before the tsunami enters the radar coverage area and ending when it reaches the coast. Section 4 describes a method to detect the tsunami, based on these radial maps.
As noted in Section 2 for a deep water amplitude of 0.5 m, the water depth at which linear wave theory breaks down is only about four meters. The equations given above will therefore apply over most of the radar coverage area, until the first surge is only moments away from the shore. As the wave moves into shore and the water depth decreases from this threshold value, the water velocity seen by the radar will increase as the group/phase velocity decreases.
3.2 Height and horizontal particle velocity profiles Tsunami waveforms typically have a weak trough followed by a large crest, trailed by a weaker trough. We used a Gaussian-modulated sinusoid to model this waveform and calculated the height and horizontal particle velocity profiles as this wave propagates towards the New Jersey coast, using the bathymetry profile shown in Fig. 2(a). Figure 4 shows results obtained by numerically solving (10) and (11) at three half-hourly intervals before the wave reaches shore. In the first, the wave is beyond the slope, more than 1.7 hours from impact. Both height and velocity curves are normalized to unity in deep water, allowing observation of their growth with time.
||As the Gaussian wave feature moves from deep to shallow water off New Jersey, Fig. 4 shows that it decreases in speed, compresses in horizontal scale, and increases in amplitude, taking nearly two hours to cross the continental shelf. Horizontal particle velocity seen by the radar increases more rapidly, as d–3/4; height increases at a slower rate, as d–1/4, consistent with Subsection 2.1.
For a tsunami amplitude of 0.5 m in 4000-m deep water, the maximum horizontal particle velocity is approximately 2.5 cm/s. As water depth decreases to 40 m, 40 km from shore, the particle velocity increases by a factor of 32, to 80 cm/s. We demonstrate in the next section that this signal is detectable at this point by a shorebased radar, 84 minutes before impact on the shore.
4. Detection Method
We applied the following procedure to demonstrate that the approaching tsunami can be detected from the simulated radial data.
(i) In the 2.5-hour period, most variation in the composite velocity field is produced by the changing tsunami currents. To reduce the effect of the varying ambient current field over this time period, we subtract the radial velocities for January 27, 2005 3:30 pm (with no tsunami signal) from those at 4:00, 4:30, 5:00, 5:30 pm. In practice, this would be done on a continuous basis.
(ii) Each composite radial velocity Vr is normalized as follows:
where the angle θ is defined as in (15).
(iii) The radar coverage area is then divided into 6 km bands parallel to the depth contours, the normalized radial velocities in a given band are averaged and the standard deviation calculated.
(iv) We then define a parameter Q(x) for each band as the absolute value of the average radial velocity normalized by the standard deviation:
When the current velocity field is dominated by tsunami velocities which are typically coherent over bands parallel to the depth contours, the value of Q(x) will increase sharply. We use the criterion Q(x) > 1 to indicate the presence of a tsunami in a given band. Q(x) is plotted in Fig. 5 both with and without the tsunami signal for the four radial velocity maps shown in Fig. 3. With no tsunami present, the value of Q(x) is
Fig. 3. Radial velocity vectors simulated for the tsunami at half-hour intervals January 27, 2005: (a) As the tsunami arrives at the
outer edge of the continental shelf, (b) 30 minutes after arrival, (c) 60 minutes after arrival, (d) 90 minutes after arrival as the
tsunami reaches the coast.
Fig. 3. (continued).
Fig. 4. Height and velocity profiles obtained by numerically solving (6) and (7) for the New Jersey depth profile shown in Fig. 2(a), normalized to their values in deep water (4000 m depth): (a) Height profile as tsunami arrives at the edge of the continental shelf, (b), (c) height profiles at successive half-hour intervals, (d) horizontal water particle velocity: (e), (f) velocity profiles at successive half-hour intervals.
|always less than unity. With the tsunami approaching the shore, it can be seen that Q(x) begins to exceed unity when the tsunami is about 125 km from the radar, where the particle velocity due to the tsunami is 4 cm/s. Q(x) increases as the tsunami advances, demonstrating that the tsunami can be detected using this simple method. It follows from (5) that 84 minutes will elapse from detection at 125 km before it strikes the coast.
Clearly, more sophisticated detection schemes could be developed than the simple method described here.
5. Warning Times at Different Locations
In this section, we further consider the time delay after detection before the tsunami strikes the shore. This depends to a large degree on the water depth over the radar coverage area. We illustrate this by calculating the time delay from detection to arrival at the coast for the Indian Ocean tsunami approaching hypothetical HF radar systems at two locations where the adjoining continental shelves are respectively broad and narrow: Penang, Malaysia, and Chennai, India. Tsunami parameters are given in Section 3. Bathymetric charts for the two locations
||are shown in Fig. 6. Figure 7 shows the water depth plotted as a function of range from the radar toward deep water. We assume that the tsunami signal is detectable in the radial velocity map when the particle velocity due to the tsunami exceeds 4 cm/s, as for Tuckerton, New Jersey (see Section 4). The maximum orbital velocities and arrival times follow from (3) and (5) using the depth profiles for the two locations shown in Fig. 7.
Off Penang in the Straits of Malacca, within the radar range (taken to be 150 km), the water depth is always less than 100 m and the corresponding particle velocity always exceeds 40 cm/s, far above the detection threshold of 4 cm/s. The tsunami will therefore be detected at the outer range of 150 km. From (5) it will be 130 minutes before it reaches the shore.
In contrast, off Chennai, the continental shelf is narrow, and the particle velocity does not reach the detection threshold of 4 cm/s until the water depth decreases to 2000 m, 60 km from the shore. With the detection method described in Section 4, the tsunami will be detected only 18 minutes before it reaches the shore, but more detailed modeling could improve this margin.
Fig. 5. The tsunami-detection parameter Q(x) plotted versus distance from shore, calculated for January 27, 2005: (a) 4:00 pm
(b) 4:30 pm (c) 5:00 pm (d) 5:30 pm. Circles: Q(x) calculated from measured SeaSonde radial velocities. Crosses: Q(x)
calculated from composite maps: tsunami currents superimposed on measured SeaSonde radial velocities.
|In a field application a feasibility study would be required for each location, based on the type of radar system planned and taking into consideration typical current regimes for the location, in addition to the bathymetry.
In conclusion, it has been demonstrated that HF radar systems in operation today are capable of detecting tsunamis and providing vital information well before impact when the adjacent continental shelf is wide. Tsunami detection software could run in parallel with software producing current velocities and local wave information at the many HF radars in operation around the coastlines of the world.
Additional methods need to be developed to predict the height and direction of the tsunami from simulated HF radar data. Simulations are also required to determine the amount of information on the tsunami that can be obtained for different radar system operational parameters. Radar data used in this study were obtained with the operational system parameters presently in use. Simulations can determine the optimal time/spatial resolution to use for tsunami detection. It is to be expected that a tsunami will be more difficult to detect if it is small or if the background current velocities are high and rapidly varying. In a field application a feasibility study would be required for each location, based on radar transmit frequency and taking into consideration typical current regimes for the location, in addition to the bathymetry.
||Figure 8 shows that shallow-water bathymetry applies over much of the area of South-East Asia that was devastated by the recent tsunami: large portions of Indonesia, Malaysia, Thailand, all of Cambodia and a good portion of southern Vietnam all are surrounded by shallow waters. This would suggest that vast areas of the Indo- Malaysian shores form an ideal region for the use of HF radar to monitor tsunamis. This study indicates that HF radar installations in such locations could provide vital capability for the detection and measurement of deadly approaching tsunamis.
We thank Ms. L. Pederson and Mr. Jimmy Isaacson of CODAR Ocean Sensors for many valuable suggestions. We are grateful to Rutgers University Coastal Ocean Observation Laboratory for the SeaSonde data used in the simulations. This was provided by through grants from ONR, NSF, NOAA and DOD.
Fig. 6. Bathymetric maps for (a) Penang, Malaysia. (b) Chennai,
India. The simulated radar position is shown by black circles.
If one reduces Eq. (8) to its 1D scalar version, one obtains:
Separation of variables is a classic method that assumes a partial differential equation (PDE) such as that above can be written in terms of factors that are each functions of only one of the two independent variables: x and t. The justification for this assumption is that it works, and does indeed produce complete solutions for the partial differential equation. Thus we assume: E(x)H(t) ≡ d(x)v(x, t). When we substitute this into (A1) above, we can group terms so that those involving x are separated from those involving t. I.e.,
Fig. 7. Water depth plotted as a function of range from the radar toward deep water. (a) Penang, Malaysia (b) Chennai, India. Note that the vertical scales for (a) and (b) differ by a factor of 10.
Fig. 8. Bathymetric map of area in South East Asia, showing extensive shallow water. Constant depth contours are labeled in meters.
|where D o is the depth at x = L, i.e., d(L) = D o, with L defining the beginning of the continental slope, i.e., the break from deep water. The group speed in deep water is given by
Each term must be set to a constant, because the variables are “independent”, and hence each can take on a range of values that do not influence the other. So set the second term equal to –ω2; hence the first term is equal to +ω2 so that they sum to zero. The solution class for the second term is the familiar trigonometric or complex exponential functions, i.e.,
which has solutions given by:
The first term then satisfies the ordinary differential equation given by (11), which has no convenient closed form solution unless the depth is constant, in which case its solution has the same complex exponential forms in space as H(t) does vs. time, i.e., traveling sinusoids.
In general, the frequency ω is an arbitrary constant, unless we want to specify time dependence to be sinusoidal. For a particular time waveform for the tsunami velocity (or height) that we might select, we can represent this as a Fourier distribution (or series if one prefers) at a variety of frequencies. Thus we represent depth times velocity (i.e., water transport) as a distribution:
where we will specify the time dependence at the distance x = L representing the beginning of the continental slope, where deep water ends and shallow begins.
where f(t) is our chosen time form to represent velocity (times depth) at the edge of the shelf. In general for arbitrary time dependence and non-constant depth profile, one uses an ODE solver to obtain E(x, ω)—and thence E(L, ω)—at an array of frequencies ω, and then Fourier transforms both sides of (A4) in order to obtain p(ω). Thence we solve (A3) via another Fourier transform (or FFT) to get the velocity profile. For the New Jersey shelf
||profile and Gaussian wavetrain, we used the MATLAB “ode45” routine, obtaining both the velocity and the height profiles vs. time in a single call, as one is the space derivative of the other, as seen by the 1D versions of (6) and (7) earlier.
As a simple example, let us assume constant depth Do and a Gaussian-modulated sinusoid such as that demonstrated in Fig. 4. Thus,
or a wave traveling in the –x direction toward shore. The form of the Gaussian wavetrain is:
where in this case we use F = 100 to produce one dominant crest surrounded by two weaker troughs. All of the integrations called for in (A4) and (A3) can be found in standard tables, and the result for the velocity profile in water of constant depth D o is:
where αo = 1/cT is the spatial wavenumber and c is wave group speed. We have normalized (A7) so that the velocity amplitude is unity.
Employing the concept of a plane wave, it is not difficult to demonstrate the strange d–3/4 depth dependence used in Section 2 for orbital/particle velocity. We base this on conservation of energy principles. As we have shown, for water of varying as well as constant depth (e.g., Fig. 4), a wave feature is compressed horizontally as the water depth decreases. Total energy must also be compressed horizontally to match the extent of the wave feature. So the standard expression for kinetic energy is KE = mv2/2, where m is mass. We can write this as:
where ρ is water density; di is the total vertical extent of the water column; λi is the length or wavelength of the feature over which energy is conserved; and vi is the horizontal velocity of the water mass. (The transverse dimension for the one-dimensional case is the same and so is ignored.) Wavelength is related to depth and period P through the shallow-water dispersion relation (4), i.e.,
|When these definitions are used and common factors are dropped, Eq. (A8) above gives the ratio of velocities for water of two different depths:
From this and the relations between height and velocity— (6) and (7)—(or based on energy conservation principles formulated around height instead of velocity), one can likewise establish the law for height increase from deeper to shallower water as:
Barrick, D. E. (1979): A coastal radar system for tsunami warning. Remote Sensing of the Environment, 8, 353–358.
Barrick, D. E. and B. J. Lipa (1986): The second-order shallow- water hydrodynamic coupling coefficient in interpretation of HF radar sea echo. IEEE J. Oceanic Eng., OE-11, 310–316.
Gower, J. (2005): Jason 1 Detects the 26 December 2004 Tsunami. EOS, 86, 37–38.
Kinsman, B. (1965): Water Waves. Prentice Hall, Inc., Englewood Cliffs, New Jersey, 676 pp.
Lamb, H. (1932): Hydrodynamics. Dover Publications, New York, 738 pp.
Mofjeld, H., V. Titov, F. Gonzalez and J. Newman (2000): Analytic theory of tsunami wave scattering in the open ocean with application to the North Pacific. NOAA Tech. Memo. OAR PMEL 116.