3.1 Exact solution to linear problem

For weak pressure distributions, $\unicode[STIX]{x1D716}\ll 1$ , the problem (2.1)–(2.5) can be linearised. The linearised version has the exact solution (Wehausen & Laitone Reference Wehausen and Laitone1960)

where $\tilde{p}(k,\unicode[STIX]{x1D713})=\exp (-k^{2}/(4\unicode[STIX]{x03C0}^{2}F^{4}))/(\unicode[STIX]{x03C0}F^{4})$ is the Fourier transform of the pressure distribution (2.5), $H(x)$ is the Heaviside function and the path of integration over $k$ is taken below the pole $k=k_{0}$ , where $k_{0}=\sec ^{2}\unicode[STIX]{x1D713}$ (thus this integral can be interpreted as a Cauchy principal value integral plus half of the residue at $k=k_{0}$ ).

3.2 Linear spectrogram

A spectrogram of a ship wake is generated by first taking a cross-section of the wave surface at a constant value of $y$ to create a wave signal, $s(t)$ , where $t=x$ by changing the reference frame to move with the uniform flow (recall the dimensionless speed is unity). The spectrogram data are then given by the square magnitude of a short-time Fourier transform:

where the window function, $h(t)$ , is an even function with compact support. In this paper we will use the Blackman–Harris 92 dB window function (Harris Reference Harris1978). The results are placed in a time–frequency heat map of angular frequency $\unicode[STIX]{x1D714}$ against scaled time $t/y$ and colour intensity on a log scale, $\log _{10}(S(t,\unicode[STIX]{x1D714}))$ .

Spectrograms computed from the linear solution (3.1) for the Froude numbers $F=0.3$ , $0.7$ , $1$ and $1.5$ are presented in figure 2 together with a solid curve which we refer to as the linear dispersion curve. More details are provided in the following subsection, but for now we note the linear dispersion curve has two branches, the lower branch corresponding to the transverse wave component of the spectrogram and the upper branch corresponding to the divergent wave component. In this figure we have chosen to fix $y=100$ , but note that for sufficiently large $y$ (say $y>25$ ), the spectrograms appear the same on this scale. On the other hand, if $y$ is chosen to be too small (that is, the sample is too close to the pressure disturbance), there will be some unwanted blurring between the two branches in the spectrogram.

There are two key features of the spectrograms in figure 2. First, the high intensity portion of the spectrograms (the lighter coloured part) in parts (a–d) appears to be centred on the linear dispersion curve. Thus, for this problem of linearised flow past a pressure distribution, the linear dispersion curve provides an excellent prediction for the dominant wave signals propagating past a representative point in space. The second key feature is that the high intensity portion is confined to the lower branch of the dispersion curve for low Froude numbers (approximately $F<0.4$ ), while for larger Froude numbers it follows the fold and the upper branch. This result is consistent with observations of ship wave patterns that suggest slowly moving ships give rise to wakes dominated by transverse waves, while fast vessels produce a wave train dominated by divergent waves. We note that this second feature is not evident if we treat the idealised problem of flow past a point pressure (as in Torsvik et al. Reference Torsvik, Soomere, Didenkulova and Sheremet2015a ), since the solution to the point-pressure problem does not depend on any parameter values apart from the pressure strength $\unicode[STIX]{x1D716}$ (there is no length scale and so no Froude number).

Figure 2.

Spectrograms of linear ship waves for the Froude numbers $F=0.3$ , $0.7$ , $1$ and $1.5$ . The solid curve is the dispersion curve given by (3.7). The colour intensity is given by $\log _{10}(S(t,\unicode[STIX]{x1D714}))$ , where $S(t,\unicode[STIX]{x1D714})$ is given by (3.2). In each case, the wave signal $\unicode[STIX]{x1D701}$ , plotted against $t/y$ , is shown above the related spectrogram.

Figure 3.

This schematic of ship wave dispersion illustrates a ship moving from the point M to the point O with a non-dimensional speed 1 over the time period $t_{0}$ . The waves generated by the ship at the point M propagate at all angles $\unicode[STIX]{x1D703}$ from the sailing line and wavenumbers $k$ with phase velocity, $c_{p}$ , and group velocity, $c_{g}$ . The large black semicircle is traced out using the group velocity for all $\unicode[STIX]{x1D703}$ . The small grey semicircle follows the same rules but is generated at a different point. The grey line, tangent to the two semicircles, represents the geometric edge of the wake system. Finally, the point I is the location of the sensor.

3.3 Linear dispersion curve

The linear dispersion curve in figure 2 is independent of the speed of the disturbance (that is, independent of $F$ ) and can be determined via the following geometric arguments. First we consider a ship moving along the path MO a distance of $t_{0}$ as shown in figure 3. At the point M the ship will generate waves in all directions $\unicode[STIX]{x1D703}$ with phase velocity

and group velocity, $c_{g}$ . To obtain the frequency at a given point $\text{I}=(t,y)$ , we must determine the angle $\unicode[STIX]{x1D703}$ and the associated wavenumber $k$ . We require the following properties: an implicit function $k(\unicode[STIX]{x1D703})$

where $\unicode[STIX]{x1D6FA}(k)=\sqrt{k}$ is the dispersion function for waves on an infinite-depth fluid; and a frequency relation given by projecting the wave vector $(-k\cos \unicode[STIX]{x1D703},k\sin \unicode[STIX]{x1D703})$ onto the direction the sensor is moving relative to the ship, $(1,0)$ , and then taking the absolute value,

Finally, we use geometry to determine $\unicode[STIX]{x1D703}$ , which can be found as a root of the equation

where $\unicode[STIX]{x1D6FC}(k)=c_{g}/c_{p}=1/2$ is the ratio of group velocity to phase velocity in a fluid of infinite depth. Solving (3.3)–(3.6) for $\unicode[STIX]{x1D714}$ gives

The above argument is essentially equivalent to that provided by Wyatt & Hall (Reference Wyatt and Hall1988), although they do not provide the formula (3.7).

The linear dispersion curve (3.7) is the solid curve in figure 2. As mention above, it has two branches, the upper branch $\unicode[STIX]{x1D714}_{1}$ and lower branch $\unicode[STIX]{x1D714}_{2}$ , and a fold where the two branches meet at the point ( $t/y$ , $\unicode[STIX]{x1D714}$ )=( $\sqrt{8}$ , $\sqrt{3/2}$ ). The upper branch corresponds to the divergent waves; it approaches the line $\unicode[STIX]{x1D714}=t/2y$ for large $t/y$ and is represented by the smaller grey circle in figure 3. The lower branch, represented by the large black circle in figure 3, corresponds to the transverse waves and approaches $\unicode[STIX]{x1D714}=1$ for large $t/y$ . The fold represents the wedge boundary of the wave train in this geometric representation, providing the well-known Kelvin’s ship wake angle, $\text{arctan}\,(1/\sqrt{8})$ (Thomson Reference Thomson1887).

4.1 Numerical scheme

There are no known exact solutions to the fully nonlinear problem (2.1)–(2.5). To compute numerical solutions we reformulate the equations using a boundary integral method, and construct a system of nonlinear equations $\boldsymbol{F}(\boldsymbol{u})=\mathbf{0}$ using collocation (Forbes Reference Forbes1989; Părău & Vanden-Broeck Reference Părău and Vanden-Broeck2002; Părău et al. Reference Părău, Vanden-Broeck and Cooker2007). The key feature of this system is that the Jacobian is fully dense, which provides challenges in terms of accuracy and run time. We have revisited this scheme recently with a preconditioned Jacobian-free Newton–Krylov (JFNK) method (Pethiyagoda et al. Reference Pethiyagoda, McCue, Moroney and Back2014a ). The preconditioner used in the JFNK method was the block-banded linear preconditioner. The solution was computed over a mesh of $N$ points in the $x$ -direction and $M$ points in the $y$ -direction, hereafter written as $N\times M$ , with mesh spacing in the $x$ and $y$ -directions given by $\unicode[STIX]{x0394}x$ and $\unicode[STIX]{x0394}y$ , respectively. Unfortunately, the largest computed domain provided by this method does not provide sufficient resolution in the frequency domain of the spectrogram to produce a reasonably accurate computed spectrogram. To overcome this limitation we have developed a subdomain stitching method described here.

The subdomain stitching method is an iterative procedure that computes subdomains one at a time and connects them together, as shown in figure 4. For clarity, subscripts are added to the parameters $N$ , $M$ , $\unicode[STIX]{x0394}x$ and $\unicode[STIX]{x0394}y$ denoting which subdomain they relate to. The first subdomain is computed on a $N_{1}\times M_{1}$ mesh that begins upstream of the pressure distribution as in the unaltered method, shown as the first white rectangle in figure 4. This subdomain is then truncated from the downstream end onto a smaller $N_{1}^{\prime }\times M_{1}$ mesh, where $N_{1}^{\prime }<N_{1}$ is the number of points in the $x$ -direction of the truncated domain. The truncation step is performed under the assumption that there exists a domain truncation error near the downstream boundary. The next subdomain is chosen such that $\unicode[STIX]{x0394}y_{1}=\unicode[STIX]{x0394}y_{2}$ , $\unicode[STIX]{x0394}x_{1}=\unicode[STIX]{x0394}x_{2}$ , $M_{1}\leqslant M_{2}$ and is placed such that the upstream border of the second subdomain matches up with mesh points from the previous subdomain. The second subdomain is placed such that the upstream points coincide with the downstream points of the truncated first subdomain (first solid red rectangle in figure 4). After the subdomain has been placed, the radiation conditions at the $M_{1}$ upstream points are changed from the algebraic decay of the original method (Pethiyagoda et al. Reference Pethiyagoda, McCue, Moroney and Back2014a ) to make sure there is continuity between subdomains (i.e. $\unicode[STIX]{x1D701}_{N_{1}^{\prime },j}^{1}=\unicode[STIX]{x1D701}_{1,j}^{2}$ , $\unicode[STIX]{x1D719}_{N_{1}^{\prime },j}^{1}=\unicode[STIX]{x1D719}_{1,j}^{2}$ , etc. for $j=1\ldots M_{1}$ where $\unicode[STIX]{x1D701}_{k,j}^{i}$ and $\unicode[STIX]{x1D719}_{k,j}^{i}$ are the surface height and velocity potential for the $k$ th point in the $x$ -direction and the $j$ th point in the $y$ -direction on the $i$ th subdomain, respectively). The integral in the boundary integral equation can then be split into an integral over a known surface, the previously computed subdomain that does not overlap with the intended solution, and an integral over the unknown surface to be computed. This procedure is repeated with $\unicode[STIX]{x0394}y_{1}=\unicode[STIX]{x0394}y_{i}$ , $\unicode[STIX]{x0394}x_{1}=\unicode[STIX]{x0394}x_{i}$ and $M_{1}\leqslant M_{2}\leqslant \cdots \leqslant M_{i}$ until a satisfactory domain size is reached. For example, figure 4 shows a free-surface profile corresponding to a solution computed with six subdomains. The preconditioner used in the JFNK method only needs to be updated for subdomain $i$ if $N_{i-2}=N_{i-1}=N_{i}$ and $M_{i-2}=M_{i-1}=M_{i}$ is not satisfied.

Figure 4.

Plan view of a nonlinear free-surface solution computed with $F=0.7$ and $\unicode[STIX]{x1D716}=0.01$ using the subdomain stitching method. The solution is on a $1625\times 551$ mesh constructed from six $301\times 551$ subdomains indicated by the dashed and solid rectangles.

When computing highly nonlinear solutions it is prudent to use a bootstrapping method which takes the solution for a smaller value of $\unicode[STIX]{x1D716}$ as the initial guess to compute the solution for a larger value of $\unicode[STIX]{x1D716}$ . This presents two methods of calculating solutions: compute the solution over all subdomains for a constant value of $\unicode[STIX]{x1D716}$ before increasing $\unicode[STIX]{x1D716}$ (subdomains first); or compute the solution for all desired $\unicode[STIX]{x1D716}$ values one subdomain at a time ( $\unicode[STIX]{x1D716}$ first). Computing subdomains first has the advantage that it enables the full solution to be calculated for a given value of $\unicode[STIX]{x1D716}$ before continuing, as opposed to the $\unicode[STIX]{x1D716}$ first method which outputs the full solutions for all $\unicode[STIX]{x1D716}$ only at the end of the procedure. The $\unicode[STIX]{x1D716}$ first method will fail to return a solution over the full domain if one of the values of $\unicode[STIX]{x1D716}$ chosen does not have a solution (or takes unreasonably long to converge to the solution); however, computing $\unicode[STIX]{x1D716}$ first is faster as the preconditioner only needs to be updated when the subdomain is changed. Therefore, the choice of method depends on whether or not a solution is known to exist for all chosen values of $\unicode[STIX]{x1D716}$ .

The main disadvantage of the stitching method described above is that spurious numerical waves are introduced at every boundary between subdomains. These numerical waves are roughly two-dimensional, small in amplitude, sinusoidal in nature and more prominent for lower Froude numbers. Due to these numerical waves, a solution constructed via the stitching method will not be in complete numerical agreement with a single domain solution. Spectrograms, however, are robust against such numerical waves. In fact, due to the sinusoidal nature of the numerical waves, they appear as a band of colour intensity at the constant value $\unicode[STIX]{x1D714}=1$ (that is, at the same frequency as the transverse waves along the centreline $y=0$ ) and so are readily identifiable and can be easily ignored.

4.2 Second-order dispersion curves

At this stage we have discussed the linear dispersion curve (3.7) which does not account for nonlinear effects. In order to include nonlinearity, we now derive an expression for further dispersion curves using a weakly nonlinear analysis. We follow the approach of Hogben (Reference Hogben1972), who proposed a second-order solution to a ship wave problem provided by a finite number of $N$ monochromatic waves and their interaction terms,

where the $a_{n}$ (for $n=1,\ldots ,N$ ) are the primary monochromatic wave amplitudes, $b_{rs+}$ and $b_{rs-}$ are non-zero amplitudes of the interacting waves and $R_{n}=k_{n}(x\cos \unicode[STIX]{x1D703}_{n}+y\sin \unicode[STIX]{x1D703}_{n})$ defines the primary wavenumber $k_{n}$ and direction $\unicode[STIX]{x1D703}_{n}$ . By considering a solution which consists of only the two monochromatic waves ( $N=2$ ) from linear theory, the amplitudes $a_{1,2}$ can easily be identified by comparing with the stationary phase approximation to the linear solution (3.1) which depends on the non-dimensional parameters $\unicode[STIX]{x1D716}$ and $F$ (the details of this calculation are not included here). The wavenumber $k_{n}$ and direction $\unicode[STIX]{x1D703}_{n}$ are given by solving (3.4) and (3.6) together to give

The amplitudes of the interacting waves $b_{rs\pm }$ are proportional to $a_{r}a_{s}$ . For further details, see Hogben (Reference Hogben1972).

To determine the location of the additional dispersion curves that arise from the second-order solution (4.1), the actual wave amplitudes $a_{n}$ and $b_{rs\pm }$ are not required (provided they are non-zero). Instead, by examining the phases of the interacting waves (that is, $R_{r}+R_{s}$ and $R_{r}-R_{s}$ ), we can easily calculate the additional dispersion curves by either doubling the linear frequencies, adding them together, or taking the difference between them:

We refer to the curves described by (4.4)–(4.5) as being the second-order dispersion curves as they come from analysing the second-order solution (4.1).

The linear and second-order dispersion curves are shown in figure 5. We emphasise that, like their linear counterparts, these second-order dispersion curves do not depend on $\unicode[STIX]{x1D716}$ or $F$ . Instead, they indicate which parts of the time–frequency domain could be highlighted in a spectrogram for a nonlinear wave pattern. Whether one or more of the branches are actually highlighted depends on the strength of the nonlinearity $\unicode[STIX]{x1D716}$ and the speed of the ship $F$ , as we see in the following subsection.

Figure 5.

The linear and second-order dispersion curves. The linear dispersion curve (with branches $\unicode[STIX]{x1D714}_{1,2}$ ) is given by (3.7). The second-order dispersion curves (with branches $\unicode[STIX]{x1D714}_{3,4}$ and $\unicode[STIX]{x1D714}_{5,6}$ ) are given by (4.4) and (4.5), respectively.

4.3 Spectrograms computed using our nonlinear simulations

We now present in figures 6–8 spectrograms computed using our fully nonlinear numerical simulations described in § 4.1 together with our linear and weakly nonlinear predictions (3.7) and (4.4)–(4.5). For all of the nonlinear spectrograms in these figures, we have fixed our $y$ -coordinate to be $y=35$ , which we have found to be sufficiently large to prevent unwanted blurring between branches of the dispersion curves.

Figure 6 shows the spectrograms computed for Froude number $F=0.2$ , representative of a slow moving ‘ship’, and the four pressure strengths $\unicode[STIX]{x1D716}=0.15$ , $0.75$ , $1.65$ and $2.25$ . Before proceeding to analyse these spectrograms, we note two features that are not of particular concern in this study. First, the low-frequency intensity present in all spectrograms in figure 6 (the horizontal band near $\unicode[STIX]{x1D714}=0$ ) is caused by the local change in surface height due to the presence of the pressure distribution, and is not part of the far-field wave train. If we had the capacity to employ much larger computational domains for our numerical solution, then we could have fixed the $y$ -coordinate to be larger than $y=35$ , which would begin to eliminate this effect. Either way, we ignore this low-frequency band. Second, a signal is visible in the spectrogram as a band of intensity at $\unicode[STIX]{x1D714}=1$ ahead of the dispersion curve (for $t/y<\sqrt{8}$ ). This effect is due to spurious numerical waves that are caused by numerical truncation and the stitching method, as discussed at the end of § 4.1.

Returning to the important trends in figure 6, we see that for low nonlinearity ( $\unicode[STIX]{x1D716}=0.15$ ), the intensity in the spectrogram (figure 6 a) occurs on the lower branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{2}$ . As the nonlinearity increases, the spectrogram intensity moves up the lower branch towards the fold of the linear dispersion curve (figure 6 b, $\unicode[STIX]{x1D716}=0.75$ ). These observations are consistent with our linear spectrograms in figure 2, which show the dominant part of the time–frequency signal for slowly moving ships lies on the transverse wave component (the lower branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{2}$ ). For even higher nonlinearity ( $\unicode[STIX]{x1D716}=1.65$ in figure 6 c and $\unicode[STIX]{x1D716}=2.25$ in figure 6 d), additional high intensity portions appear along the second-order mode $\unicode[STIX]{x1D714}_{4}$ , which comes from transverse waves interacting with themselves (recall $\unicode[STIX]{x1D714}_{4}=2\unicode[STIX]{x1D714}_{2}$ ). There is even a hint of a further signal at $\unicode[STIX]{x1D714}=3\unicode[STIX]{x1D714}_{2}$ , which would represent an even higher-order mode. Finally, the extremely nonlinear solution in figure 6(d) exhibits multiple modes appearing in a vertical section in line with the fold ( $t/y=\sqrt{8}$ ). These multiple modes bear a striking resemblance to the so-called leading waves shown in figure 1. Thus we can conclude that the leading wave component identified by Torsvik et al. (Reference Torsvik, Soomere, Didenkulova and Sheremet2015a ) is possibly due to nonlinearity (steep nonlinear waves).

Figure 6.

Spectrograms of nonlinear ship waves for non-dimensional pressure strength number $\unicode[STIX]{x1D716}=0.15$ , $0.75$ , $1.65$ , $2.25$ and Froude number $F=0.2$ . The solid curves are the linear and second-order dispersion curves given by (3.7), (4.4) and (4.5). The colour intensity is given by $\log _{10}(S(t,\unicode[STIX]{x1D714}))$ where $S(t,\unicode[STIX]{x1D714})$ is given by (3.2). In each case, the wave signal $\unicode[STIX]{x1D701}$ , plotted against $t/y$ , is shown above the related spectrogram.

For $F=0.7$ , representative of a faster moving ‘ship’, the maximum intensity of the spectrogram for the linear solution is around the fold of the linear dispersion curve (see figure 2 b). This behaviour is replicated in the nonlinear spectrograms in figure 7(a), which is for the small value $\unicode[STIX]{x1D716}=0.01$ . For moderate nonlinearity, shown in figure 7(b), an additional mode of intensity appears along $\unicode[STIX]{x1D714}_{3}$ and $\unicode[STIX]{x1D714}_{5}$ . Further increasing nonlinearity allows for a clearer realisation of the additional modes and the appearance of colour intensity along $\unicode[STIX]{x1D714}_{6}$ (figure 7 c). For a highly nonlinear solution, $\unicode[STIX]{x1D716}=0.15$ (figure 7 d), the region of intensity distorts away from the linear dispersion curve, leading to high intensity regions of the time–frequency map appearing to the left of the fold (that is, for $t/y<\sqrt{8}$ ). Thus, for high nonlinear flows past a pressure distribution, there exist a visible part of the wave train that appears outside of Kelvin’s wedge, as discussed by Pethiyagoda, McCue & Moroney (Reference Pethiyagoda, McCue and Moroney2014b ) (that is, some highly nonlinear solutions have apparent wake angles that are greater than Kelvin’s angle).

Figure 7.

Spectrograms of nonlinear ship waves for non-dimensional pressure strength number $\unicode[STIX]{x1D716}=0.01$ , $0.05$ , $0.11$ , $0.15$ and Froude number $F=0.7$ . The solid curves and colour intensity are the same as in figure 6.

Figure 8 is for $F=1$ , representing an even faster ‘ship’. Here the spectrogram follows the same trend as $F=0.7$ with a clearer colour intensity along the second-order dispersion curves $\unicode[STIX]{x1D714}_{3,5}$ and an absence of colour intensity along $\unicode[STIX]{x1D714}_{6}$ . Thus, for highly nonlinear flows due to faster ships, we expect additional high intensity portions of the spectrogram along the second-order modes that come from divergent waves interacting with themselves ( $\unicode[STIX]{x1D714}_{3}=2\unicode[STIX]{x1D714}_{1}$ ) and divergent waves interacting with transverse waves ( $\unicode[STIX]{x1D714}_{5}=\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{2}$ ). Further, we see that the distortion observed in figure 8(d) occurs later than that observed for $F=0.7$ (figure 7 d). The different locations for the distortions is consistent with the location of the greatest intensity in the spectrograms for the exact linear solution.

Figure 8.

Spectrograms of nonlinear ship waves for non-dimensional pressure strength number $\unicode[STIX]{x1D716}=0.01$ , $0.05$ , $0.11$ , $0.15$ and Froude number $F=1$ . The solid curves and colour intensity are the same as in figure 6.

We close this subsection by commenting on the ultimate fate of nonlinear solutions as our parameter $\unicode[STIX]{x1D716}$ increases. From Bernoulli’s equation (2.2), we see there is an absolute upper bound for the free-surface height, $\unicode[STIX]{x1D701}_{\mathit{upper}}=1/2$ . Thus, for a given solution, the maximum surface height $\unicode[STIX]{x1D701}_{\mathit{max}}<\unicode[STIX]{x1D701}_{\mathit{upper}}$ . This maximum surface height typically occurs at the crest of one of the divergent waves, although it can also occur on the centreline, depending on the parameter values. As $\unicode[STIX]{x1D716}$ increases, the problem becomes more nonlinear, and we find the maximum surface height $\unicode[STIX]{x1D701}_{\mathit{max}}$ also increases. Qualitatively, we observe that for increasing $\unicode[STIX]{x1D716}$ , the waves themselves have sharper crests and broader troughs. We speculate that ultimately there exists a critical value, $\unicode[STIX]{x1D716}_{c}(F)$ , such that $\unicode[STIX]{x1D701}_{\mathit{max}}\rightarrow \unicode[STIX]{x1D701}_{\mathit{upper}}$ as $\unicode[STIX]{x1D716}\rightarrow \unicode[STIX]{x1D716}_{c}$ , although our present numerical scheme is not capable of computing solutions that are sufficiently nonlinear to explore this issue in any detail.

4.4 Comparison with experimental data

We shall now make some comparisons with the experimental spectrogram shown in figure 1. The high-speed ferry that produced the wake in question is named the Star, and sails from Tallinn, Estonia to Helsinki, Finland. The reported operating speed of the Star at the time the data were measured was $14.2~\text{m}~\text{s}^{-1}$ (Parnell et al. Reference Parnell, Delpeche, Didenkulova, Dolphin, Erm, Kask, Kelpšaite, Kurennoy, Quak and Räämet2008). Its length and width are 186 m and 27.7 m, respectively. This corresponds to a half-length based Froude number of $F_{L}\approx 0.47$ and half-width based Froude number of $F_{W}\approx 1.21$ . Obviously, the ferry has an aspect ratio which is much larger than unity. However, it is still instructive to compare with our mathematical model which assumes an axisymmetric disturbance.

Figure 9.

A scaled version of figure 1, with the dispersion curves given by (3.7), (4.4) and (4.5) overlaid. The intensity of the spectrogram closely follows the linear and second-order dispersion curves $\unicode[STIX]{x1D714}_{1-6}$ within the leading wave component of the spectrogram and follows the divergent wave portion of the linear dispersion curve $\unicode[STIX]{x1D714}_{1}$ outside the leading wave. The transverse wave component of the spectrogram does not appear to follow the lower branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{2}$ , for large $t/y$ .

Considering the spectrogram in figure 1, we now overlay the linear and second-order dispersion curves in order to compare our theoretical results with experimental data. To do so we scale the axes based on the speed of the ship, $U$ , and the minimum distance to the sensor, $y$ , and align the time the ship is closest to the sensor with $t=0$ . As the exact speed, distance and passing time are not known, we match the divergent wave intensity and fold location to the upper branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{1}$ . The resulting scaled spectrogram is shown in figure 9. Note that in this example we do not consider the transverse wave component $\unicode[STIX]{x1D714}_{2}$ when matching between our theoretical and experimental results, because the lower branch of the linear mode in the experimental spectrogram in figure 1 does not appear to be horizontal for large $t/y$ , as the theory predicts. We return to this point in § 5.

By matching the experimental spectrogram as just described, we find the speed of the ship to be $15.75~\text{m}~\text{s}^{-1}$ , the distance from the sensor to be 2.5 km and the time the ship is closest to the sensor to be at 21:18. Comparing the properties of the ship’s voyage with the reported operating speed and the closest distance to the sensor given by Torsvik et al. (Reference Torsvik, Soomere, Didenkulova and Sheremet2015a ), we see that our calculated speed is a slight over estimation of the reference value, $14.2~\text{m}~\text{s}^{-1}$ . This could be due to a number of issues, for example the effects of a steady underlying current or finite-depth effects that subtly change the shape of the dispersion curve. More encouragingly, having a clear divergent wave component means our estimated distance to the sensor falls in the reference bounds of 2.5–3 km, where 2.5 km corresponds to the outgoing shipping lane.

Visually, figure 9 shows very good agreement between the dispersion curves and the experimental spectrogram in the time frame of the leading wave. In particular, for the region near $t/y=\sqrt{8}$ (roughly $2<t/y<4$ ), the high intensity part of the spectrogram follows the linear dispersion curve (made up of branches $\unicode[STIX]{x1D714}_{1,2}$ ) and the second-order curves (with branches $\unicode[STIX]{x1D714}_{3-6}$ ). This exercise shows how important it is to understand the consequences of steep nonlinear waves when interpreting spectrograms. It is worth noting that the experimental spectrograms presented by Wyatt & Hall (Reference Wyatt and Hall1988) (calculated from a tug boat named Quapaw) show very similar features, including the prominent high intensity regions in the leading wave component. Indeed, Wyatt & Hall (Reference Wyatt and Hall1988) comment that these properties are likely due to nonlinearity.

Finally, we note that to the right of the leading wave region of the spectrogram in figure 9 (for roughly $t/y>4$ ), there is an obvious disagreement between the experimental results and the lower (transverse) branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{2}$ . We consider this issue further now.

5.1 Linear dispersion curve with acceleration

One feature of the experimental spectrogram in figure 9 that is not yet explained is the obvious difference between the transverse branch of the linear dispersion curve, $\unicode[STIX]{x1D714}_{2}$ , and the transverse component of the experimental spectrogram data. In an attempt to explain this discrepancy, we consider the case that the ship is accelerating from rest up to its cruising speed. We are motivated to explore this approach for two reasons: a ship travelling at a slower speed will generate transverse waves of higher frequency (from the dispersion relation); and the transverse waves are caused by waves generated earlier than those that generate divergent waves.

We illustrate the second property by referring to figure 10(a,b). In figure 10(a), the solid curve represents a hypothetical ship speed versus time for a $y$ value that is fixed to be the closest distance to the sensor. The actual curve is found by fitting to data, as explained below. The points A, B and C in the figure are associated with waves that are generated at a time $t_{gen}<0$ (figure 10 a) and detected by the sensor at a later time $t_{sen}>0$ (figure 10 b). With this particular velocity profile, the wave represented by the point A in figure 10(a,b), generated at $t_{gen}\approx -6.05$ and detected by the sensor at $t_{sen}\approx 18.34$ , falls on the transverse branch of the dispersion curve. The wave represented by the point B is generated later than A, at $t_{gen}\approx -4.4$ , but still lies on the transverse branch of the dispersion curve. It is detected earlier by the sensor, at $t_{sen}\approx 5.61$ . Waves generated at later times will eventually lie on the divergent wave branch, such as the wave represented by C. Thus, from this argument we see that if the ship is travelling at a slower speed when wave A is generated compared to when wave B is generated, the frequency of the wave A along the dispersion curve will be greater than the frequency of wave B if both waves fall on the transverse branch of the dispersion curve. Therefore, the transverse branch of the linear dispersion curve will grow in frequency as $t/y$ increases if the ship is accelerating.

Figure 10.

(a) A plot of example data (squares) for non-dimensional ship speed and the approximately fitted velocity function (5.3). (b) A plot of the comparison between the linear dispersion curve for an accelerating ship (5.2, solid) and constant-speed linear dispersion curve ( $\unicode[STIX]{x1D714}_{1,2}$ , dashed). (c) The dispersion curves computed for the accelerating ship overlaid on the scaled version of figure 1. The labels A, B and C in (a) and (b) each identify a wave generated at the ship at a time in (a) and received by the sensor at the time in (b).

To determine the location of the linear dispersion curve for an accelerating ship, we first define the non-dimensional displacement of the ship, $X(t)$ , and its velocity, $U(t)$ , with the following properties:

Here $t=0$ corresponds to the time at which the ship is closest to the sensor. Thus, under these properties, the ship is accelerating during an interval of time before $t=0$ but is moving with constant speed by the time $t=0$ . The dispersion curve can be defined parametrically in terms of $\unicode[STIX]{x1D703}$ by first determining where the ship generated the wave, $X_{gen}=-y\cot \unicode[STIX]{x1D703}$ , and thus the time the wave was generated, $t_{gen}\leqslant 0$ such that $X(t_{gen})-X_{gen}=0$ . The time taken for the wave to reach the sensor can then be calculated by $t_{wave}=y\,\text{cosec}\,\unicode[STIX]{x1D703}/c_{g}$ where the group velocity $c_{g}=c_{p}/2=U(t_{gen})\cos \unicode[STIX]{x1D703}/2$ . Finally, the variable-speed version of (3.4) and (3.5) given by $k(\unicode[STIX]{x1D703})=\sec ^{2}\unicode[STIX]{x1D703}/U^{2}$ and $\unicode[STIX]{x1D714}=Uk(\unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D703}$ , respectively, are used to derive the dispersion curve

5.2 Comparison with experimental data

The specific velocity function $U(t)$ will affect the shape of the new dispersion curve (5.2), thus in order to match with the experimental spectrogram (figure 1) we require a rough approximation for the velocity of the ship observed in the experiment. We used the shipping traffic website run by the organisation MarineTraffic (MarineTraffic.com) to track the velocity of the Star on a particular day as it left port at Tallinn and travelled to the point closest to the sensor (the location of the sensor is given by Parnell et al. Reference Parnell, Delpeche, Didenkulova, Dolphin, Erm, Kask, Kelpšaite, Kurennoy, Quak and Räämet2008). The data we obtained are non-dimensionalised by scaling the velocity by the cruising speed, $U_{cruise}$ , and by scaling time by $U_{cruise}/g$ . Time is shifted so that $t=0$ corresponds to when the ship was closest to the sensor. An example of the data (for a day in February 2016) is given as the squares in figure 10(a). It appears the ferry initially accelerates as it leaves port then slows down as it turns onto its sailing line before finally accelerating up to its cruising speed. For this specific dataset $U_{cruise}\approx 12.86\,\text{m}~\text{s}^{-1}$ . Considering only the final acceleration phase, the velocity data for this particular ship look approximately like an error function, thus we represent it by

where $t_{shift}=-535/74$ and $\unicode[STIX]{x1D6FD}=100/37$ are parameters that are determined by roughly fitting to the data. The velocity profile (5.3) with these parameters is the solid curve in figure 10(a).

Figure 10(b) shows a comparison between the new linear dispersion curve (5.2) with $U(t)$ given by (5.3) for an accelerating ship (solid line) and the constant-speed linear dispersion curve ( $\unicode[STIX]{x1D714}_{1,2}$ , dashed). The transverse branch of the dispersion curve for the accelerating ship clearly increases in frequency as $t/y$ increases; conversely, the divergent branch is very close to the constant-speed dispersion branch. The second-order dispersion curves for the accelerating ship have been derived using the method in § 4, and the full set of dispersion curves has been overlaid onto the experimental spectrogram in figure 10(c), with the same axis scaling as figure 9. We see that the new linear dispersion curve does a much better job of predicting the transverse wave in the spectrogram. Of course the fitting process used here is very rough, but our goal is only to suggest a reasonable approximation of the velocity $U(t)$ on a given day (the data for the ferry observed in 2008 are no longer freely available). Given the clear improvement in the comparison between the measured spectrogram and the theoretical prediction, we are confident that the increase in frequency in the measured transverse wave is due to the observed ship accelerating before passing by the sensor (and not due to nonlinearity, say). We acknowledge that other explanations are possible, for example shoaling or other finite-depth effects such as wave refraction. Or perhaps another explanation entirely.