Models for simulation of ocean surface waves, such as WW3, are formulated in terms of parameterizations for the nonlinear wave-wave interactions, S_nl (f, θ), with other source terms, such as wind input to waves, S_in (f, θ), and wave dissipation, S_ds (f, θ). WW3 is based on the well-known balance equation for spectral action density, expressed as n(f,θ), where (f, θ) denote frequency and direction; it is an open-source modern 3^rd generation wave model (WW3DG, 2016). Version 5.16 is used in this study. Detailed discussion of WW3 model physics and characteristics is given by WW3DG (2016) for both regional and global applications. As noted in the Introduction, the basic three source terms in deep water are the nonlinear wave-wave term S_nl (f, θ), wind input, S_in (f, θ), and wave dissipation, S_ds (f, θ). As described by Perrie et al. (2013), the implementation of TSA within WW3 follows the methodology used in implementing DIA or WRT in third generation wave models like WAM and WW3, respectively (WAMDI Group, 1988; Van Vledder, 2006; Tolman, 2009). We assume an explicit forward – time scheme for the difference equations, for the 2-dimensional ocean wave spectrum F(f, θ). In terms of the nonlinear wave-wave interactions S_nl , the integration is semi-implicit, requiring a diagonal term to estimate S_nl at succeeding time-steps. This term is the diagonal of the partial derivative of S_nl (f, θ) with respect to spectral energy F(f, θ), where f, θ are spectral frequency and direction. Thus, only array elements with equal f and θ in both the source and spectrum terms are used; the diagonal term may be written as Λ_i , where

which must now be determined for TSA.

In evolving wind and wave conditions, the TSA formulation works well in many sea state conditions, for example, simple fetch-limited or duration-limited wave growth. However, when the wind direction changes rapidly, the TSA formulation needs to be modified to allow the broad-scale term to take into account more complicated spectral situations such as multi-peaked spectra. If there is a significant misalignment of TSA’s broad-scale term with respect to waves generated by sudden wind direction changes, the original TSA formulation of Resio and Perrie (2008) may not able to provide a reliable representation of the nonlinear transfer, S_nl .

In its original formulation, TSA is based on the WRT formulation for S_nl (f, θ), which is due to Webb (1978); Tracy and Resio (1982) and Resio and Perrie (1991). In this approach, the nonlinear transfer of action density from one spectral wavenumber k₃ to another, k₁, is represented by a transfer function T(k₁,k₃),

which may be expressed as

where ϑ is the Heaviside function, k₄=k₁+ k ₂- k ₃ where k ₂= k ₂(s, k ₁, k ₃). Here, n_i is the action density at k_i , and the locus of possible wave-wave interactions is specified by the contour s satisfying the resonance condition,

and where ŋ is the local orthogonal to contour s.

In the original TSA of Resio and Perrie (2008), a given spectrum n_i is decomposed into a 1^st order, or broad-scale, term n ^ i and a 2^nd order, local-scale term n ′ i , where n ^ i is given a parametric JONSWAP-type form following Hasselmann et al. (1973), thereby depending on only a few parameters. The local-scale term n ^ i > is the residual,

with the same number of degrees of freedom as the input spectrum n_i . TSA becomes quite accurate if the parameterization for n ^ i can be selected so that n ^ i is small. However, to have optimal computational efficiency, the methodology will try to minimize the number of parameters used for n ^ i because application of large multi-dimensional sets of pre-computed matrices for n ^ i is time-consuming.

For complicated wave conditions, for example multi-peaked spectra n_i , the challenge of selecting n ^ i hat is dependent on a relatively small parameter set can be somewhat mitigated by application of multiple broad-scale terms, or for example in this paper, two broad-scale terms,

where n ^ 1 i and n ^ 2 i are both given JONSWAP-type parametric forms, corresponding to two peaks, f_{p 1} and f_{n 2}, in the given spectrum n_i . Therefore, following the original TSA methodology, the residual local-scale n ′ i term may be determined by subtracting,

Thus, in the usual manner of Resio and Perrie (2008), we partition the action density term, n_i , and write the transfer integral T in equation (3) in terms of the sum of interactions involving broad-scale terms, which we pose as n ^ 1 i + n ^ 2 i , denoted B, local-scale terms n ′ i , denoted L, and the cross interactions among n ^ 1 i + n ^ 2 i and n ′ i , denoted X. Thus, the nonlinear transfer interactions S_nl can still be represented as,

and B can be pre-computed and depends on JONSWAP-type parameters x_i ,

The objective of this slightly generalized TSA is to accurately, efficiently approximate L+X, by neglecting terms involving n ′ 2 and n ′ 4 thereby simplifying equation (8). This follows Resio and Perrie (2008) in assuming that the local-scale terms ( n ′ 2 and n ′ 4 ) are deviations around the associated broad-scale terms ( n ^ 2 and n ^ 4 ) which capture most of the spectral energy; and with their positive/negative differences and products, the former tend to cancel, as we move along the interaction loci. This approach is validated by mTSA’s ability to give results that compare well with those of WRT.

Thus, Resio and Perrie (2008) show that eliminating n ′ 2 and n ′ 4 gives

where N ∗ 3 is given by

and they use known scaling relations to obtain

with

where the so-called ‘pumping’ and ‘diffusion’ terms are Λ_p and Λ_d following Webb (1978)’s notation. Here, superscript p is the equilibrium-range power law, for example f^{- 4} or f^-5, and (ζ/ζ₀) is related to a linear scaling coefficient for the n ^ i terms; (k/k₀) is the ratio of the spectral peak wavenumber for the spectrum being integrated to that of the reference spectrum. Coordinates θ_* and k_* are

For a f^-5 JONSWAP spectrum, ς is the Phillips’ coefficient in equation (12), whereas for an f^-4-type spectrum like E(f) ≈ βf^-4, then ς is β, and generally any linear multiplicative term that scales the spectrum. The power of ς is the number of broad-scale densities ( n ^ i ) in the integrals used in matrices, Λ_d and Λ_p . The scaling factor for wavenumber k is from the wavenumber dimensions of the coupling coefficient (~k⁶), Jacobian (~k^1/2), and the phase space terms (dskdk~k³).

From equations (12) and (13), the diagonal terms for WRT are,

and neglecting of terms involving n ′ 2 and n ′ 4 , we find for TSA, or in this case, mTSA,

These terms are central to the mTSA semi-implicit implementation within WW3.

Operational wave models like WW3 are restricted in the sense that they have a finite discrete spectral grid. mTSA’s broad-scale terms n ^ 1 i + n ^ 2 i typically depend on a few parameters for each broad scale term, such as peak frequency f_p , peak direction, θ_p , Phillips coefficient β, peakedness γ, spectral width parameters σ_a and σ_b for the forward and rear faces of the spectral peak, and a spreading distribution, ~cos^m (θ-θ_p ) around the spectral peak direction θ_p . However, when the spectral peak f_p is too close to the highest frequency of the discrete spectral computational grid, it is not possible to define β in terms of the equilibrium range of the spectrum, in the usual manner,

where the equilibrium range is assumed as ~2 or 3 × f_p , and C_g is the group velocity (Donelan et al., 1985). In these cases, a simple practical approach is to define β in terms the highest discrete frequency above f_p , and below the equilibrium range, which may be nonexistent in this case (more on this below). This is an approximation in terms of the expected value for β; had the frequency grid extended to a higher limit with an equilibrium range, a more accurate estimate would be possible. In this way, a modified definition of β allows the WW3 forecast model to continue with the computation, providing an estimate for β. This approach is consistent with that of previous third generation wave models. However, the issue regarding the calculation of β can become critical if there are multiple spectral peaks, particularly regarding the region between two spectral peaks.

The mTSA approach allows generalization of the broad-scale term, allowing more than one broad-scale parametric term, corresponding to multi-peaked spectra. It is shown in test cases in the next section, that this modification allows the mTSA approach to rather accurately represent a fully integrated formulation of the nonlinear wave -wave transfer, in terms of the WRT estimates. For example, in test cases involving windsea-swell interactions, with two or more changing spectral peaks, each with differing peak frequency directions, we represent the double peak as two separate peaks in the mTSA formulation, and invoke the broad-scale parameterization twice to simulate each separately. The procedure for determining JONSWAP parameters for two broad-scale terms for a double-peaked spectrum is given in the Appendix.

An issue is the high-frequency equilibrium range. When there is only a single peak frequency, the equilibrium range is typically about 2~3 times the peak frequency f_p and is represented by spectral tail with f ^-4 variation, following Resio et al. (2004), and earlier studies. The f ^-4 spectral tail is matched to the upper limit of the discrete frequency spectrum, which allows an equilibrium range Phillips coefficient to be defined, denoted β in the notation of Resio et al. (2004). This is a key term in the broad-scale parameterization. Moreover, there are clearly instances where a given wind-wave spectrum may have a spectral tail which does not follow the f^-4 distribution. But that is not a problem for the mTSA methodology. In any case, whatever the distribution of the spectral tail, whether f^-4 or f^-5 or f^-4.5 etc., we construct the broad-scale spectrum using a JONSWAP-type parameterization as discussed in this section and in the Appendix. Whatever mismatch occurs is then reconciled by the residual local-scale term, as identified in equation (7).

In the mTSA formulation, while there may be two or more peaks, we represent the spectrum in terms of just two peaks in this paper, where the second peak f_p2 is required to be more than 2 frequency bins higher than the first peak f_p1 ; otherwise we represent the spectrum with a single broad-scale term, and mTSA becomes just the standard TSA formulation. For the second peak f_p2 , the upper part of the spectrum above f_p2 is used to define the equilibrium range and the Phillips coefficient β, in the usual manner. If there are enough frequency bins between the two peaks, we define an equilibrium range ~ 2-3×f_p2 , and calculate a Phillips-type coefficient, denoted here as “β“. Otherwise we just use the highest frequency bin above f_p2 , as a proxy to define the β term, to allow the simulation to proceed.

The region between the two peaks (f_p1 and f_p2 ) is a challenge. The peaks must be separated by at least one frequency bin, to allow definition of two broad-scale terms. In this case, if there are enough frequency bins between f_p1 and f_p2 , we define an equilibrium range ~ 2-3×f_p2 , and calculate a Phillips-type coefficient β for the broad-scale term associated with f_p1 , in the usual manner. However, if f_p1 and f_p2 are separated by only a very small region, then for the lower of the two spectral peaks, f_p1 , we represent the energy of its equilibrium range by the minimum spectral energy between f_p1 and f_p2 , thus defining β for the lower of the two spectral peaks. Thus, each of the two peaks, f_p1 and f_p2 , has its own β coefficient, peak direction, θ_p1 and θ_p2 , and its own broad-scale term.

For an assumed two-peaked spectrum, each of the two broad-scale terms requires a directional spreading distribution function. As with the original TSA formulation, the directional spreading distribution for each of the two broad-scale terms is assumed to be of the form ~cos^m (θ - θ_p ) where the integral exponent m is selected so that the broad-scale spreading at the respective spectral peak (f_p1 or f_p2 ) can approximate that of the given input spectrum F(f_p,θ), for whatever directional distribution this is.

The initial test case considers a sheared spectrum, with swell propagating to the west at 0°, and higher frequency wind-waves at higher frequency propagating to the north at 270° as shown in Figure 1C. Within WW3, the convention for winds is always the Meteorological Convention; direction from, clockwise from North. For waves, it is direction to, counterclockwise from East.

Figure 1 shows comparisons for the 1-D (1-dimensional) and 2-D action density and nonlinear transfer, S_nl , for the three formulations, namely DIA, multiple-TSA (denoted mTSA), and WRT. We do not show the single TSA (denoted sTSA) case, because as in the more basic JONSWAP-type cases considered by Perrie et al. (2013), the broad-scale term n ^ i is able to fit the swell spectral peak region rather well (Figure 1A), and therefore sTSA and mTSA are essentially the same in this case.

However, this case shows that there is a notable mismatch between n ^ i and the given test case action, n_i , in the high-frequency region, where the energy is almost zero. Thus, the role for the local-scale term n ′ i is relatively important in this latter region. Resultant 1-D estimates for the nonlinear transfer S_nl given by mTSA are able to match those of WRT well, compared to DIA, as shown in Figure 1B, as well for 2-D results shown in Figures 1D–F. By comparison, DIA results have magnitudes that are too large in the positive and negative lobe regions as shown in Figure 1B. While detailed patterns for mTSA and WRT are shown to compare well in Figures 1D, E, DIA results appear distorted in Figure 1F, particularly in the high frequency portion of the spectrum.

The second test case consists of simply letting the sheared spectrum in case 1 evolve, without wind input or wave dissipation, or other source terms, like the propagation of swell waves. As the nonlinear transfer is conservative, no change in total energy is expected. 2-D results are presented in Figure 2 after time evolution of 5 hours. Here, we see that, as in the initial conditions given in Figure 1, estimates for the nonlinear transfer S_nl given by both sTSA and mTSA are able to match those of WRT well, compared to DIA. As in the first test case, compared to the results from TSA and WRT, (which compare well with each other), DIA results have detailed patterns that appear distorted, particularly in the high frequency portion of the spectrum. The implications of these differences particularly appear to show up as energy growth in high frequency regions of the spectrum as indicated in Figure 2, e.g. the smaller waves that might be central to satellite backscattering from synthetic aperture radars.

The third test case considers the same initial wave spectrum as in the first case in Figure 1; however there is now a constant west-to-east wind blowing opposite the main westward swell direction, at 20 m/s, with an initial secondary wind-sea to the north, as shown in Figure 1, orthogonal to the assumed wind direction. See Figure 3. The main low-frequency swell direction is to the west. The formulations for wind input, S_in , and wave dissipation, S_ds , are given by the ST4 source terms of Ardhuin et al. (2010), as implemented within WW3.

Results are shown after the system has evolved for 10 hours. The wind speed of 20 m/s is relatively strong. After 10 hours of time evolution, the new wind-generated waves are the dominant feature in the spectrum of this system. However, in each of the simulations shown in Figure 3, we see that westward propagating swell remains mostly unchanged by the nonlinear transfer S_nl , regardless of which formulation is used, WRT, sTSA, mTSA, or DIA. Minor variations are obtained in the swell results, due to DIA formulation and more so in the sTSA results, compared to results from WRT or mTSA.

By comparison, results for the eastward propagating wind-sea driven by 20 m/s wind, imposed on the initial conditions of a northward wind-sea, and a westward propagating swell are a different story. We see that the results for the 2-D action density and the nonlinear transfer S_nl , given by mTSA, are largely able to compare to those of WRT relatively well, compared to those of DIA or sTSA. Differences between results of WRT and mTSA are comparatively minor.

By comparison, DIA results have variations in frequency regions around the spectra peak of the northward-propagating wind-sea region. In terms of direction, these effects are most notable in the northeastward direction, approximately diagonal between the northward propagating wind sea, and the new eastward-propagating wind-generated waves. Results from sTSA show notable biases throughout much of the spectrum, although the overall shape of the 2-D action density is rather similar to that of mTSA. For varying 2D distributions, it is not completely clear how to assign a single number to express error. If we take maximum Hs as a kind of qualitative expression of mismatch, then relative to WRT, mTSA has an error of about 1.9% in maximum Hs, compared to 8.7% for DIA.

The overall dominance of the new wind-generated waves propagating to the east is clear after 10 hours. The effects of nonlinear wave-wave interactions between the new wind-generated waves propagating to the east, and the initial conditions involving wind-sea propagating in the northern direction are relatively minor.

The fourth test case is similar to the third case, except now the 20 m/s wind is blowing south to north, orthogonal to the main east-to-west swell direction, and parallel to the initial secondary wind-sea, also to the north. See Figure 4. The formulations for wind input, S_in , and wave dissipation, S_ds , are given by the ST4 source terms of Ardhuin et al. (2010), as implemented within WW3. The simulation is for 5 hr.

As in test case 3, we see that westward-propagating swell remains largely unchanged by the nonlinear transfer S_nl , regardless of the formulation; WRT, sTSA, mTSA or DIA. Minor variations are obtained in the swell spectrum due to the DIA formulation, particularly in the southwestern direction, and more so for the results from sTSA, which shows apparent ‘smoothing’ of the westward swell spectrum. Results for the northward propagating wind-sea imposed on the initial conditions and the new generating north-propagating wind-sea waves are dominant features of the simulations.

As in the previous test cases (Figures 1–3), we see again that simulation results given by mTSA are able to match those of WRT rather well, compared to those from DIA or sTSA. By comparison, DIA results have more northerly-propagating wind-sea region and more directional spreading than results suggested by mTSA and WRT; these similar tendencies for more directional spreading and more smoothing are also notable in results from sTSA. The initial wind-sea in the northern direction is effectively ‘assimilated’ into the new wind-generated waves propagating in the northern direction. There is no apparent impact of the west-propagating swell on the newly-generated north-propagating wind-sea.

As in test case 3, the wind speed of 20 m/s is relatively strong, and after 5 hours of time evolution, the new wind-generated waves are the dominant feature in the 2-D spectrum of this system. The effect of the west-propagating swell on the new wind-sea is rather minor and only evident in directional components in the northwest direction, adjacent to the swell propagation direction. However, variations due to differing formulations for S_nl are evident. We see that after 5 hours evolution, the results for the simulation given by mTSA are able to match those of WRT relatively well, compared to DIA or sTSA.

By comparison, DIA results have magnitudes that are too large, in both the negative high frequency regions, and also in frequencies of the region around the spectra peak of the north-propagating wind-sea region. In terms of direction, the DIA effects result in more broadly distributed directional spreading of the new wind-sea, rather than to have dominantly north-propagating waves. This tendency of bias and excessive directional spreading is accentuated in the sTSA results. In terms of error, if we take maximum Hs as a qualitative expression for mismatch, then relative to WRT, mTSA has an error of about 2.0% in maximum Hs, compared to 4.3% for DIA.

The overall dominance of the new wind-generated waves propagating northward is clear after the 5-hour time evolution. The effects of nonlinear wave-wave interactions between the new wind-generated waves propagating to the north, and the initial conditions involving swell propagating to the west are relatively minor.

The computational domain for implementation of WW3 for the simulation of waves generated by Hurricane Teddy consists of the nested grid system shown below in Figure 6. This nested system has a relatively coarse-resolution (0.5°) large-scale grid which extends from 20°N to 65°N, and from 40°W to 75°W for the Northwest Atlantic. Within this domain, a relatively high-resolution subdomain is nested from 42°N to 52°N and from 55°W to 72°W, focused on the waters off northern New England and the Canadian Maritimes Provinces, as shown in Figure 6. These grid resolutions are selected to provide a relatively reliable degree of accuracy in simulating swell and wind-generated waves energy. The directional resolution is 10° and 29 frequency bins are used, spaced logarithmically using f_{n +1} = 1.10 f_n ranging from 0.04118 Hz to 0.5939 Hz. The model global time step is 600 s. As mentioned earlier, the ST4 source terms are used following Ardhuin et al. (2010), for wind input S_in and wave dissipation S_ds .

Wind fields to drive the wave models for Hurricane Teddy are obtained from Environment and Climate Change Canada (ECCC). For the coarse-grid Northwest Atlantic domain, 3-hourly ECCC global wind products are used with 0.24 resolution on a latitude-longitude grid. For the high-resolution Atlantic Canada domain, 1-hourly ECCC regional wind products are used with 10 km resolution based on a polar-stereographic projection. These are routine forecast products that are posted daily by ECCC (https://dd.weather.gc.ca/). The ECCC global wind data are already on a latitude-longitude grid, and so need no further processing. The ECCC regional model has wind components on the polar-stereographic grid of the weather model simulation, and need to be rotated to our latitude-longitude reference frame in order to be ingested into WW3, which then performs interpolations in space and time as needed.

Figure 6 compares significant wave height distributions from WW3 using the four formulations implemented for the nonlinear wave-wave transfer S_nl term. These are the three formulations used in this study; DIA, WRT, and mTSA. A fourth simulation is denoted mTSA4, which uses the tuning for ST4 for the Northwest Atlantic as determined by Perrie et al. (2018); the latter was a study of waves generated by three intense nor’easter storms, and different wave models implemented on coarse- and fine-resolution nested grid systems that are similar to those used in Figure 6. The tuning of ST4 reflects regional characteristics of the Northwest Atlantic and the southern Gulf of Maine, and was shown to give enhanced results.

In these simulations, WW3 is driven by the ECCC wind fields as Teddy propagated from the area around Bermuda, to Nova Scotia as described in the previous section. In terms of the Hs spatial distributions, the results from DIA and WRT appear to exhibit slightly larger Hs values that those of mTSA. However, as expected, this can be compensated for by applying the ST4 tuning suggested by Perrie et al. (2018), which results in higher values in the simulation results for mTSA4.

Overall, the qualitative features of the Hs area distributions and wave directions are similar in these four simulations, with respect to the propagation of the wave fields, the overall directional patterns of the waves, and the peak wave directions. Consistent with fetch-limited growth results reported in previous studies by Perrie et al. (2013), mTSA appears to give results that are biased low, whereas results from DIA are biased high compared to those of WRT, using the standard ST4 source terms. These differences in results from respective different S_nl formulations may be somewhat modulated by the inherent nonlinearity, numerical instability etc. present in the WW3 model system.

To estimate the reliability of the model simulations using the different nonlinear transfer S_nl formulations, we conducted comparisons with measured significant wave heights, Hs, with observations at four buoys in the model domain whose locations are shown in Figure 5B. Comparisons were made between time series of buoy measurements of significant wave heights, Hs, and the model simulations for Hurricane Teddy, using the implemented S_nl formulations; DIA, mTSA, mTSA4 and WRT. Buoy observations are generally reported at 30-minute intervals.

Figure 7 shows Hs time series and scatter plots at four selected buoys along Teddy’s storm track. These are buoy 41049 near Bermuda, 44008 off Cape Cod, 44137 off Nova Scotia and 44139 near the southern tip of the Grand Banks. These four were selected in order to have additional discussion regarding 1-D and 2-D spectra, in the sections that follow. A summary of statistics for specifically these 4 buoys, plus an additional 10 other buoys along, or near, the storm track of Teddy, indicated in Figures 5A, B is presented in Table 1, in terms of root mean square error (RMSE), bias, correlation coefficient (corr), and scatter index (SI). In this computation of statistics, only data within ± 3 days of the passage of the storm by a given buoy are used; these data are indicated in Figure 7. Note that the full observational time series are indicated by red stars (*), with the points used to calculated statistics (within ± three days of the peak) marked by black circles ʘ overlaying the red stars (*).

In terms of capturing the storm peak values and model biases, the Hs time series and scatterplots in Figure 7 appear to suggest similar behaviors for all formulations, DIA, mTSA, mTSA4 and WRT. Statistical indices in Table 1 also reflect this finding, with mTSA4 tending to overpredict the peak Hs values and mTSA tending to underpredict these values. This suggests the general approach, that in implementing WW3 for specific regional applications, it important to perform some careful tuning of the basic ST4 source terms of Ardhuin et al. (2010), to reflect the associated regional characteristics. This trend is also reported in the statistical indices in Tables 1 and 2 for the 14 buoys shown in Figure 5B along, or near, Teddy’s storm track. mTSA4 is able to outperform mTSA in terms of improved RMSE, reduced bias, improved correlation coefficient, and scatter index. However, although results from mTSA4 appear to compare somewhat favorably with measurements at buoys 41049 and 44008, the mTSA4 results at buoys 44137 and 44139 appear to have notable overestimates. Thus, the performance is not unequivocal.

Because ST4 was originally tuned for DIA by Ardhuin et al. (2010) for simulation of waves for global ocean studies it tends to perform well. In fact, statistical indices in Table 1 suggest that DIA can outperform the other simulations; those using mTSA, mTSA4 and WRT. Occasionally, for some of the statistical indices, mTSA4 can outperform the other simulations, for example at buoy 44066 at the edge of the Continental Shelf off the coast of Delaware USA.

As reported by Perrie et al. (2018) the tuning of ST4 consists of adjustments to parameters BETAMAX, the wind-wave growth parameter, and ZALP, the wave age shift of the long waves to account for gustiness, respectively, 1.75 and 0.008, to give optimal simulation skill for waves generated by three nor’easters for the Northwest Atlantic and Gulf of Maine region. It is anticipated that additional tuning could also produce more improvements to the performance of mTSA4. But that is not the objective of this study. In any case, it is interesting to compare the results from mTSA4 with those from mTSA and WRT. Table 2 suggests that for the four selected buoys, results from mTSA4 are generally better than those from WRT or mTSA.

Comparisons between 1-D wave spectra at the peak of the buoy measurements are shown in Figure 8 for Hurricane Teddy. The simulated 1-D spectra are all essentially dominated by single-peaked spectra. Some of the comparisons show large discrepancies between model simulations and observed data. It is evident that as the peak of the storm passes the buoy locations, each will experience changing wind directions, and therefore interactions between swell and wind-waves are present. These interactions are evident in the comparisons that are shown.

We consider the simulations from the different S_nl formulations. Buoy 41049 off Bermuda shows a dominant peak near 0.07 Hz and a secondary peak near 0.09 Hz. Only mTSA4 shows any (minor) indication of the secondary peak. All simulations underestimate the primary peak. mTSA provides a notable underestimate compared to the observations or compared to the other simulations like mTSA, DIA and WRT.

Results for buoy 44008 off Cape Cod are somewhat similar to those of 41049. The observed data indicate some indication of a secondary spectral peak at about 0.11 Hz, and a primary peak at about 0.07 Hz. All model simulations tend to also have indications of the occurrence of the secondary peak, e.g. a ‘wiggle’, but none simulate well the details suggested by the observed data. Also no simulation captures the primary peak, although mTSA4 and DIA appear to come close to so doing, whereas results from WRT and mTSA increasingly underestimate the observed peak data, respectively.

Results for buoy 44137 off Nova Scotia are similar to those of 41049, with a dominant peak around 0.06 Hz and a secondary peak at about 0.08 Hz. Different from the results at 41049, here the results from mTSA4 appear to capture the secondary peak, but notably overestimate results at the primary peak, as do results from WRT. By comparison, results from DIA and mTSA appear to provide a somewhat favorable simulation of the primary spectra peak.

Results for buoy 44139 on the Grand Banks are notable because of the double peak, a low frequency possible swell peak at about 0.06 Hz, and a higher wind-waves peak at about 0.09 Hz. In this case, all the simulations appear to provide some indication of the secondary wind-waves peak, although all present underestimates, with mTSA4 giving the best simulation. For the primary peak, mTSA4 provides an overestimate, whereas the other three manage to give somewhat reasonable simulations.

Comparisons between observed 2D wave spectra and model simulations for Hurricane Teddy at about its peak are shown in Figure 9 for buoys 41049 and 44008. Results for the buoy measurements are calculated following the Longuet-Higgins approximation for the Fourier expansion method as recommended by the NDBC website (www.ndbc.noaa.gov/measdes.shtml). The observed data in Figure 9 at both buoys show the response of the wave spectra to turning winds, as the hurricane passes by and as the primary peak modulates to the new direction of the developing wind-waves.

Model simulations suggest qualitatively similar results. Although the observed main directions of the low frequency primary peak and the developing secondary peak are approximately consistent with the simulations, the modelled maximum energy values are generally overestimated compared to the observed values, 92.6 m²Hz^-1rad^-1 at buoy 41049 and 38.3 m²Hz^-1rad^-1 at buoy 44008. The model simulations also appear to provide results with wider distributions of new wind-wave energy than is being generated in the new developing wind direction, compared to more restricted directional spreading suggested by the observed spectra. Often it is the other way around, with rather wide directional distributions reported for buoy wave data compared to narrow distributions estimated by wave models. An example of the latter can be found in the comparisons of 2-D modelled and measured spectra by Perrie et al. (2018), which may be attributed to the Longuet-Higgins approximation for the Fourier expansion method.

Overall, the directional distributions resulting for the three S_nl formulations do not differ significantly, except in terms of the magnitude of the simulated 2D spectral peaks compared to the observed data. At buoy 41049, magnitudes of peak 2D spectral values are approximately the same as observed, and mTSA4 is too high. See Table 3. At buoy 44008, magnitudes of peak 2D spectral values are approximately the same as observed for results from DIA and WRT, whereas mTSA is too low and mTSA4 is too high. Therefore, as mentioned before, although the ST4 tuning of the two parameters, BETAMAX, the wind-wave growth parameter, and ZALP, the wave age shift of the long waves to account for gustiness, may improve simulations of Hs in Tables 1 and 2, and the 1D spectra in Figure 7, that is not always the case for the 2D spectra.

Another possible constraint on the models is numerics, in terms of the shifting of direction of the spectral peak, and spectral direction distributions. WW3 uses third-order upwind propagation. This is the mechanism that can contribute to the model’s ability to shift the dominant wave directions in response to changing wind directions.

The presentation of results is not complete without mention of computational efficiency. However, the focus of this study has been computational accuracy, rather that computational efficiency. The new mTSA code has not been optimized with MPI (Message Passing Interface) or other methodologies, whereas WRT has had such optimizations. Therefore, in its present formulation, mTSA does not run efficiently. For example, whereas mTSA allows a very large reduction in the number of computations needed to approximate the full integration for the Boltzmann integral for the wave-wave interactions, the separation within the spectrum is presently quite demanding and has not been optimized.

A summary of computational efficiency of mTSA relative to DIA and other formulations for S_nl is given in Table 4. In this comparison, FBI is the full Boltzmann integration representation of these quadruplet interactions, which is similar to WRT, and has been used extensively in earlier comparison studies of TSA, such as in Resio and Perrie (2008), with similar run times, ~110 × DIA. By comparison, the present mTSA methodology is about ~100 × DIA, whereas previous older parameterizations of these formulations, FBI-4 and TSA-4, which incorporate alternating frequency and angle computational loops to accelerate the efficiency, and additional parameterizations to attempt improved accuracy, have a computational efficiency in the range of about ~ 26 to 30 × DIA. Future work will focus on optimizing mTSA and the need to enhance computational efficiency.