Tsunamis are long waves that evolve substantially, through spatial and temporal spreading from their source region. Here, we introduce a new analytical solution to study the propagation of a finite strip source over constant depth using linear shallow-water wave theory. This solution is not only exact, but also general and allows the use of realistic initial waveforms such as N-waves. We show the existence of focusing points for N-wave-type initial displacements, i.e. points where unexpectedly large wave heights may be observed. We explain the effect of focusing from a strip source analytically, and explore it numerically. We observe focusing points using linear non-dispersive and linear dispersive theories, analytically; and nonlinear non-dispersive and weakly nonlinear weakly dispersive theories, numerically. We discuss geophysical implications of our solutions using the 17 July 1998 Papua New Guinea and the 17 July 2006 Java tsunamis as examples. Our results may also help to explain high run-up values observed during the 11 March 2011 Japan tsunami, which are otherwise not consistent with existing scaling relationships. We conclude that N-waves generated by tectonic displacements feature focusing points, which may significantly amplify run-up beyond what is often assumed from widely used scaling relationships.
The devastating effects of tsunamis, near- and far-field, became widely recognized following the 2004 Boxing Day tsunami. Run-up measurements over the periphery of the Indian Ocean showed considerable variation in near- and far-field impact. This variation, while it can be inferred from direct numerical simulations of the evolution of the initial wave, remains largely counterintuitive .
Once generated, tsunamis evolve substantially through spreading as seen in figure 1, for the 2011 Japan event similar to the 2004 tsunami. In their map of global propagation patterns of the 26 December 2004 tsunami, Titov et al.  observed two main factors affecting directionality: the configuration of the source region  and the waveguide structure of mid-ocean ridges . Continental shelves also act as waveguides  and likely caused the persistent ringing for the Pacific coasts of South and North America in the 2004 tsunami. Similar effects were observed during the 2010 Chilean and the 2011 Japan tsunamis; in both cases, strong currents persisted in California ports for days [8,9].
Most studies of past tsunamis have concentrated on modelling specific historic or scenario events and have invariably focused on numerically estimating coastal effects and devastation. The standard practice is for initial conditions to be established directly from estimates of the seismic parameters, then their subsequent evolution is deterministic. Coastal inundation involves often supercritical overland flow depths [10,11], as observed in the numerous videos of the 11 March 2011 Japan tsunami, and remains the most temperamental aspect of computations, as small coastal features affect flooding patterns to first order .
It is thus not surprizing why analyses of directivity and focusing in the open ocean remain few. During the 2011 Japan tsunami, Guam, 2750 km from the source, experienced maximum run-up of less than 0.60 cm, while in Irian Jaya, 4500 km away, the reported run-up reached 2.6 m. It is our objective here to supplement the few existing substantive studies of the physics of deep-sea evolution of tsunamis and suggest that observations that are counterintuitive may be explainable through the classic field theory.
Tsunamis triggered by submarine earthquakes have a finite crest (strip) length, which is believed to be calculable adequately from estimates of the seismic parameters and scaling relationships. The initial profile is dipole-shaped, directly reflecting regions of uplift and subsidence of the seafloor, yet the generated wavefield exhibits finger-like radiation patterns, a process often referred to as directivity. Ben-Menahem  defined a directivity function using the source length and the rupture velocity. Later, Ben-Menahem & Rosenman  used linear theory to calculate the radiation pattern from an underwater moving source and showed that tsunami energy radiates primarily at a direction normal to a rupturing fault. Okal et al.  reported field observations of the 1946 Aleutian tsunami in the far-field, and concluded that a large slow earthquake and a landslide must have occurred concurrently to have caused the observed far-field distribution and near-field run-up. Okal  then identified differences in directivity patterns between tsunamis from landslides and dislocations.
Directivity arguments alone, however, cannot explain the complexity of the radiated patterns in oceans with trenches and seamounts. Berry  discovered how such underwater features may concentrate tsunamis into cusped caustics, causing large local amplifications at specific focal points. He used linear dispersive theory that describes fairly accurately the evolution of tsunamis in the open ocean. Nonlinear effects become important only as the tsunami evolves over near shore topography, through shoaling and refraction.
In our analysis, we will not develop exact results for onland inundation; when needed to model few interesting past events, we will use numerical methods to calculate the run-up. We note that, in terms of 1+1 dimensional run-up, Synolakis  developed an exact solution to both linear and nonlinear shallow-water theory for the canonical problem of a non-periodic long wave climbing up a beach. Synolakis & Skjelbreia  then discussed the evolution of the maximum during shoaling. Carrier & Yeh  developed an analytical solution based on the methodology defined by Carrier  to evaluate the propagation of finite crest length sources of Gaussians over flat bathymetry and discussed the directivity. Their solution involves computation of complete elliptic integrals of the first kind, with singularities, as also in Carrier et al. , a difficulty subsequently resolved by Kânoğlu  and Kânoğlu & Synolakis . Moreover, in its current form, the Carrier & Yeh  model can only be used for Gaussians, and cannot even be applied for more standard long wave models such as solitary waves. Their solution can be extended for N-waves by superposing Gaussians, but the calculation of the singular integrals is challenging and involves approximations.
We note that Tadepalli & Synolakis  proposed a paradigm change for analytical studies of the impact of long waves and introduced N-waves as more realistic initial waveforms for tsunamis. Indeed, the initial waveform of real tsunamis is dipole-shaped, with leading-elevation or -depression waves depending upon the polarity of the seafloor deformation and the observation location. Marchuk & Titov  described the process of tsunami wave generation by rectangular positive and negative initial ocean surface displacements, and their results suggested unusual amplification.
In summary, we will first present a new general analytical solution for the linear shallow water-wave equation for propagation of a finite crest length source without any restriction on the cross section of the initial profile, i.e. Gaussian, solitary or N-waves, and without singular elliptic integrals . We will then show that focusing points exist for N-wave-shaped sources and persist in the corresponding dispersive solutions. Last, we will apply our analytical solution to the 17 July 1998 Papua New Guinea (PNG), the 17 July 2006 Java and the 11 March 2011 Japan events to explain some extreme run-up observations.
We will then examine focusing and large local amplification, not by considering the effects of underwater diffractive lenses, which anyway are calculable since Berry , but by considering the dipole nature of the initial profile. We do not of course purport to explain specific patterns of real transoceanic tsunamis, but only suggest that, in addition to the Berry focusing from bathymetric lenses, large amplification might result from focusing dependent on the shape and orientation of the initial wave.
We note that our analytical solution is not intended to replace numerical models that are inadvertent for identifying the impact of scenario events or historic tsunamis. Just as analytical results for idealized problems help establish the scaling of natural phenomena, far easier than repeated numerical computations over large parameter ranges, our basic wave theory analysis helps interpret puzzling field observations.
2. Analytical solution
We use the linear shallow-water wave equation to describe a propagation problem over a constant water depth d as a governing equation. In terms of the free surface elevation η*(x*,y*,t*), the dimensional governing equation is 2.1where g is the gravitational acceleration. Dimensionless variables are introduced as 2.2Here, l0(d0)=d and are the characteristic length (the depth), and the time scales, respectively. The dimensionless form of the governing equation then takes the form 2.3with an initial surface profile η(x,y,t=0)=η0(x,y) and zero initial velocity ηt(x,y,t=0)=0.
The Fourier transform pair over the space variables (x,y) is 2.4a,bwhere , η=η(x,y,t), and k and l are the wave numbers in the x and y directions, respectively. We transform the governing equation (2.3) into 2.5and the initial conditions to and , using (2.4a). The solution of (2.5) under these conditions is now straightforward; where Back-transformation through (2.4b) gives 2.6
We express the finite-crested initial waveform as the product of two independent functions, 2.7as in Carrier & Yeh . Here, f(x) describes the transverse extent of the initial wave profile, g(y) represents the streamwise cross section of the source, such as a solitary wave or an N-wave (figure 2). Representation of the initial wave in the form of (2.7) is advantageous, because it allows evaluation of the Fourier transforms of f(x) and g(y) independently; 2.8
Unlike Carrier & Yeh , who used Gaussians, we prefer hyperbolic functions to define lateral cross-sectional profiles g(y) such as solitary or N-waves. We also use a hyperbolic function in the transverse direction to define finite crest length f(x), i.e. 2.9In (2.9), x0 is the starting point of the source and L is its crest length, as shown in figure 2. The parameter γ in (2.9) is determined by the lateral cross section of the initial wave; it is either γ=γs for solitary waves or γ=γn for N-waves. The factor is included so that the amplitude of f(x) is equal to unity, in the limit and , for a given γ. In that case, the problem reduces to a single propagation direction, and (2.7) represents an infinitely long source. Given that 2.10(see appendix A for details), the transform of (2.9) takes the following form: 2.11
Spectra for cross-sectional profiles g(y) have been given by Synolakis  for solitary waves, and in Tadepalli & Synolakis  for generalized N-waves. A solitary wave with amplitude H can be described by gs(y)=Hsech2γs(y−y0) with . Its transform is given by Synolakis  as , with αs=π/(2γs). Consequently, an initial finite crest wave with solitary wave cross section can be described with 2.12and its transform by 2.13The generalized N-wave profile is defined by Tadepalli & Synolakis  as , where ε is a scaling parameter which ensures that the initial wave amplitude is H. The steepness of the wave is controlled by the parameter p0 in . The locations of depression and elevation parts of an N-wave are controlled by y1 and y2. The transform of the generalized N-wave is given by Tadepalli & Synolakis  as 1 with αn=π/(2γn). As a result, an initial finite crest wave with an N-wave cross section is 2.14with the corresponding transform, 2.15In appendix B, we revisit the so-called axisymmetric water wave problem analysed by Carrier & Yeh , then their extension to the finite crest length solution, and we address shortcomings of their solution.
In addition, we follow the linear dispersive analytical solution of Kervella et al.  for the potential flow equation.2 The solution is similar to (2.6) except the dispersion relation and it is given by 2.16where , and, again with
3. Results and discussion
In what follows, we will use realistic initial wave profiles to show the existence of focusing points. Then, we will discuss possible consequences of focusing points for the 17 July 1998 PNG, the 17 July 2006 Java and the 11 March 2011 Japan tsunamis. We will then compare our results obtained with linear non-dispersive theory with results from linear dispersive, nonlinear non-dispersive and weakly nonlinear weakly dispersive theories.
When a typical N-wave initial source (figure 3a) propagates over a constant depth, the initial wave profile splits into two outgoing waves as presented in figure 3b,c, i.e. leading-elevation N-wave (LEN) and leading-depression N-wave (LDN). This is consistent with the inferences of Tadepalli & Synolakis [25,28] and also with field observations, for example, those after the 26 December 2004 Great Sumatran tsunami. In Male, Maldives the tsunami manifested itself as an LEN, as elsewhere to the west of the Sumatran subduction zone. In Phuket, Thailand, it manifested itself as an LDN, as elsewhere to the east of the subduction zone . More interestingly, however, as the LEN and LDN travel in opposite directions, in the path of the LDN, a positive wave from the centre of elevation part and two positive waves from the sides of depression arrive simultaneously at a point along the bisector line, as shown in figure 3b,c. This is the focusing point, and in its vicinity abnormal tsunami wave height is observed (figure 3d).
We now study specific N-wave initial source with H=0.001, L=30, p0=15, y1=0, y2=2.3 and ε=0.04 and compare two (1 spatial+1 temporal)- and three (2 spatial+1 temporal)-dimensional propagation results in figure 4. In both cases, LDN and LEN propagate in opposite directions. However, while two-dimensional propagation results show that the initial wave splits into two waves with identical elevation and depression heights propagating in opposite directions, as expected from classic linear wave theory, three-dimensional propagation produces waves propagating with different elevation and depression heights in each direction, along the bisector. Moreover, because of focusing, the wave height increases at first on the leading-depression side, and then decreases monotonically. On the leading-elevation side, the decay is monotonous.
In addition, we compare analytical solutions of linear non-dispersive (2.6) with linear dispersive (2.16) theories, and with numerical solutions of nonlinear non-dispersive (MOST) and with weakly nonlinear weakly dispersive3 (Zhou et al. ) theories (figure 5). The focusing points persist in predictions using all four approximations of the governing equations of hydrodynamics, and the differences among the four are almost indiscernible, in these parameter ranges of geophysical interest.
(a) The 17 July 1998 Papua New Guinea tsunami
We now consider the 17 July 1998 PNG tsunami, an iconic catastrophe that brought into worldwide attention the impact of submarine landslides . A crucial feature for the landslide hypotheses was the unusually high run-up values observed over a fairly small coastal area, an observation which led to the development of source discriminants . We will now suggest that the extreme run-up values might have also been due to focusing of the LDN.
We used the landslide source suggested by Synolakis et al. , with an initial wave of approximately −18 m depression followed by the +16 m elevation, and source length of L=1, as shown in figure 6a, and present the maximum wave height distribution over the entire flow field in figure 6b. Focusing is apparent, manifesting itself as a second local maximum, the global being the maximum of the initial waveform.
In figure 6c, we present the maximum wave height envelopes along different directions for the wave propagating into the leading-depression side, where we observe the focusing points in each direction. Sissano Lagoon was the area with the maximum impact where the maximum number of casualties was reported. This region is approximately r=25 km away from the source, between 30° and 45° from the source orientation close to the 30° radial line (figure 6d). As seen from figure 6c, along 30°, the shoreline would face a wave 1.4 times larger than if the region was r=50 km away. Using the run-up variation of Kânoğlu , also implied by Tadepalli & Synolakis , where R∼H3/4, the equivalent run-up on the target shoreline could be 1.3 times greater.
We then investigate the effect of the source crest length (L) has on the location and amplitude of the focusing point (figure 6e). As L increases, the focusing point moves further away, as expected, since focusing is an effect of the finite crest length. Overall, the maximum wave height along the leading-depression wave side is higher than along the direction of the leading-elevation wave side, at any point past the focusing point. Also, increasing L up to a certain value increases the maximum wave height at the focusing points, then it is constant. In figure 6f, we present the effect of the steepness parameter over the location of the focusing point and the maximum wave height value at the focusing point. Decreasing the wave steepness translates the focusing point further away. However, it does not change the maximum wave heights at the focusing points.
In addition, we investigate the effects nonlinearity and dispersion have on focusing points, using the PNG initial condition as a test case. When the initial wave is steep (p0=15), the wave is more dispersive, and we observe a slight increase of the maximum at the focusing point in dispersive solutions (figure 7). However, all four approximations of shallow-water wave theory—linear non-dispersive, linear dispersive, nonlinear non-dispersive, weakly nonlinear weakly dispersive—produce almost identical results when p0=2, as seen in figure 8.
(b) The 17 July 2006 Java, Indonesia tsunami
The magnitude Mw 7.8, 17 July 2006 earthquake off the south coast of Western Java, Indonesia generated an unexpectedly large tsunami. The tsunami affected over the 300 km of the coastline with over 600 casualties. A pronounced run-up peak exceeding 20 m was measured at Permisan (figure 9a), although overall run-up values ranged from 5 to 7 m at surrounding areas . They inferred that the focused run-up heights were suggestive of a local submarine slump or mass movement, following the invariants presented by Okal & Synolakis . Here, we suggest another explanation for such a pronounced run-up.
As discussed in the previous section, a wave from the centre of the elevation side and two waves from the sides of the depression (shown by the arrows in figure 9b) arrive simultaneously at the focusing point. Therefore, we consider Permisan as the approximate focal point, and we present reverse tsunami travel time (RTTT) contours, i.e. the time it will take for a source over a contour to reach the specific shoreline point, Permisan in this case (figure 9b). Over the RTTT plot, we also included the tsunami source functions from the propagation database of the United States National Oceanic and Atmospheric Administration’s Center for Tsunami Research (NCTR). Detailed explanation of the NCTR’s tsunami propagation database and its application to real-time forecasting can be found in Tang et al.  and Wei et al. .
We thus identified a scenario where positive waves from the centre of the elevation part of the source and boundaries of the depression part arrive simultaneously and fit the RTTT contours well. The arrival time is approximately 40 min, consistent with accounts of eyewitnesses , as shown in figure 9b.
We then used MOST with its Community Interface for Tsunamis (ComMIT)  to simulate the Javan tsunami with the proposed source over the real bathymetry/topography and found the maximum wave amplitudes as shown in figure 9c. Because of limited bathymetric and topographic data for the region, it was not possible to make satisfactory computations to evaluate the run-up distribution along the shoreline. The model results suggest approximately 5 m run-up at Permisan. However, our computations over a constant depth bathymetry clearly shows that Permisan is close to a focusing point, and possibly the unusual run-up partly explainable with classic field theory rationalizations.
(c) The 11 March 2011 Japan tsunami
Finally, we examine whether focusing amplified the coastal effects of the 2011 Japan tsunami using the source distribution identified in real time by Tang et al. . We note that, as to the time of this writing, the unusually high run-up in Northern Japan remains unexplained, with some having already suggested a co-seismic landslide. In figure 10a, we show the maximum wave height distribution and also included maximum wave height contours from propagation of the wave over a constant depth of 5000 m. This is done to exclude bathymetric focusing effects. The results presented in figure 10a reveal that focusings occur just in front of the shorelines where maximum run-up heights were measured (figure 10b).
We considered three-dimensional long wave propagation over a constant depth basin. We solved the linear shallow-water wave equation as an initial value problem subject to realistic initial conditions, such as N-waves. We showed the existence of focusing point for evolution of dipole sources.
We then discussed unusually pronounced run-up observations from recent events which might have been exasperated by focusing, using the 1998 PNG and the 2006 Java tsunamis as particular examples. Our results strongly imply that focusing increases the shoreline amplification of the tsunami. Moreover, as shown in the Javan case, our methodology was useful in identifying the tsunami generation mechanisms considering the focusing location and tsunami travel times. Also, preliminary examination of the field survey data and analysis of the 25 October 2010 Mentawai Islands, Sumatra, Indonesia tsunami  suggests a similar focusing effect. We additionally presented initial modelling results for the 11 March 2011 Japan tsunami as indicative of possible focusing.
We note that the three-dimensional focusing identified here supplements the focusing by bathymetric lenses proposed by Berry  and points to a very rich and complex evolutionary process for long waves in the ocean. While for most practical applications, numerical solutions will be used, event-specific computations can seldom help identify basic wave phenomena. Both focusing mechanisms suggest extreme care when attempting to further regularize ill-posed seismic inversions from run-up measurements, using simplistic multiplicative factors to scale sources. Furthermore, our analysis suggests that the occurrence of extreme run-up at a given locale should not always automatically lead to the inference of offshore landslides, at least not before focusing is eliminated.
Ending, we dedicate this paper to the memory of the late Prof. D. H. Peregrine. In the first meeting of the Tsunami Risk and Strategies for the European Region (TRANSFER) project, which was supported by the European Union, he suggested that it would be useful to compare numerical solution results with an analytical solution for three-dimensional spreading, so that dispersive characteristics of numerical schemes can be investigated. We thus suggest that the analytical solution developed here could also be used as a benchmark, in addition to the analytical solutions presented in Synolakis et al.  for establishing the veracity of the rapidly multiplying newly minted tsunami numerical models.
U.K. is indebted to NOAA’s Pacific Marine Environmental Laboratory for support during the every summer visits from 2004 to 2012. U.K. and B.A. express gratitude for the support of BAP-08-11-DPT2002K120510 in Middle East Technical University. U.K., B.A. and C.E.S. acknowledge the partial support from the Scientific and Technological Research Council of Turkey project no. 109Y387 and from the General Secretariat for Research and Technology, Greece project no. 10TUR/1-50-1 in the framework of the joint research program between Turkey and Greece. T.S. acknowledges the support of Frederick Dias through grants from the ESDP of ENS-Cachan, the cultural service of the French Embassy in Dublin and the Strategic and Major Initiatives scheme of University College, Dublin. This work was partially supported by the TRANSFER project, EU contract no. 037058, FP6-2005 Global-4. This study was partially funded by the Pacific Marine Environmental Laboratory (PMEL) and the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under NOAA cooperative agreement no. NA10OAR4320148, contribution no. 3806 (PMEL), no. 1857 (JISAO). C.E.S. gratefully acknowledges the support of the National Science Foundation of the US under grant nos CMMI-0928905, CMMI-1135768 and CMMI-1313839.
We split into real and imaginary parts, i.e. A1
to evaluate the Fourier transform of f(x), equation (2.9). The first integral on the right-hand side (r.h.s.) vanishes as the integrand is an odd function evaluated over a symmetric interval. The second integral has an even integrand, hence, (A1) reduces to A2The integral on the r.h.s. is the Fourier sine integral of and it is readily available in integral tables, or can easily be evaluated by using computer algebra systems such as Mathematica as A3provided Re γ>0, which is satisfied for the present problem. Hence, A4
Appendix B. Comparison with the solution of Carrier & Yeh 
Carrier & Yeh  first developed an analytical solution for axisymmetric waves—single Gaussian humps—in order to analyse the directivity of finite crest sources. They chose as characteristic length and depth scales the source breadth l0=W—breadth of the Gaussian hump or elongated hump—and the source amplitude d0=H (figure 2) and the time scale . The dimensionless form of the governing equation (2.1) is the same as (2.3), with this scaling. Carrier & Yeh  introduced the change of variable and obtained the axisymmetric equation B1from (2.3). They developed an axisymmetric solution to (A1) using a Gaussian hump η(r,t=0)=2 e−r2=2 e−(x2+y2) with zero initial velocity ηt(r,t=0)=0, as B2However, they observed that the resultant integral becomes inconvenient to compute for large distances (r) and times (t), although it is well behaved for small values of r and t. Hence, following Carrier et al. , they expressed the solution as B3with B4where is the complete elliptic integral of the first kind. They showed that the solution of (A1) exhibits a self-similar behaviour (t∼r). They further replaced the complete elliptic integral with modified Bessel functions, through trial-and-error, using the self-similarity of the solution. Then they extended their axisymmetric solution to an elongated source defined by, B5
and found the radiation pattern from a strip when t>20, as B6in which B7where I1/4 and K1/4 are the modified Bessel functions of the first- and second-kind of order , respectively.
With the simple methodology we presented in §2, we find an exact analytical solution which does not resort to singular elliptic integrals for the axisymmetric problem. Given a Gaussian hump η0(x,y)=2 e−r2=2 e−(x2+y2), as in figure 11a, its Fourier transform is , where again , and thus the exact axisymmetric wave solution becomes B8through (2.6). The time evolution of Gaussian hump is obtained from (A8) and given in figure 11b, which compares its results with numerical solution predictions. Our results not only show an exact comparison with Carrier & Yeh  and MOST computed results, but also we had no difficulty in direct numerical integration of (A8). In addition, we could observe the self-similarity (t∼r) of the solution, exactly as used by Carrier & Yeh  to extend the axisymmetric wave solution into a solution for a strip source.
The extension of our axisymmetric solution to a strip source is also straightforward. The Fourier transform of the elongated source (A5) used by Carrier & Yeh  is . Thus, our solution integral is given by B9In figure 12, we consider an elongated source with fixed length L=20. Time series of water surface elevation are evaluated along different angular directions (θ=0°, 45° and 90°), through direct numerical integration of (A9). The effect of the source length L on wave height distribution is demonstrated in figure 13. We let L vary, and we compute time series along directions normal (θ=0°) and parallel (θ=90°) to x-axis. Our results are identical to those of Carrier & Yeh  and are also verified through numerical modelling. Moreover, our analytical solution is simple, yet versatile in terms of application of different initial conditions such as single humps, solitary waves and N-waves expressed in the form of Gaussian and hyperbolic functions. Furthermore, our solution involves no approximation and no curve fitting.
↵2 Instead, a variant of shallow-water wave theory which captures the effect of dispersion to the lowest order of the Boussinesq approximation  can also be used, i.e. It produces almost identical results as the potential flow solution.
↵3 By weakly nonlinear, we mean the model which retains the lowest-order nonlinear terms. Similarly, weakly dispersive refers to the model that considers the lowest-order dispersive terms.
- Received January 9, 2013.
- Accepted January 30, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.