The North Sea Andrea storm and numerical simulations

A coupling of a spectral wave model with a nonlinear phase-resolving model is used to reconstruct the evolution of wave statistics during a storm crossing the North Sea on 8–9 November 2007. During this storm a rogue wave (named the Andrea wave) was recorded at the Ekofisk field. The wave has characteristics comparable to the well-known New Year wave measured by Statoil at the Draupner platform 1 January 1995. Hindcast data of the storm at the nearest grid point to the Ekofisk field are here applied as input to calculate the evolution of random realizations of the sea surface and its statistical properties. Numerical simulations are carried out using the Euler equations with a higher-order spectral method (HOSM). Results are compared with some characteristics of the Andrea wave record measured by the down-looking lasers at Ekofisk.


Introduction
A number of studies addressing rogue waves have been conducted theoretically, numerically, experimentally and based on the field data in the last decade.The occurrence of rogue waves, their generation mechanism, and detailed dynamic properties are now becoming clear.The state-of-the-art review on extreme and rogue waves can be found in recent review papers and books such as Socquet-Juglard et al., (2005), Dysthe et al., (2008), Kharif et al., (2009), Osborne (2010), Slunyaev (2010) and Onorato et al., (2013).
Unfortunately, there are few studies available based on field data, partly due to the limited number of rogue waves recorded in the ocean.Investigations of meteorological and oceanographic (met-ocean) conditions, in which extreme and rogue waves occur, together with field analyses of wave time series are of importance for a getting better insight of the mechanisms generating these abnormal waves.
It should be noted, however, that field data are usually recorded in 17-30 min periods every third hour and therefore are affected by sampling variability (uncertainty due to limited number of observations) making it more difficult to drawn firm conclusions from a field data analysis (see Bitner-Gregersen andHagen, 1990, 2003).Further, met-ocean conditions between each 3 h measurement are assumed to be stationary due to lack of information about their variability within the 3 h time period.There are some locations where continuous measurements of sea surface elevation are taken; however, they are spare.It is important to note that sea states recorded continuously every 17-30 min are usually not remaining stationary; therefore, not justifying the combination of 17-30 min wave records in one longer wave record.Thus, sampling variability will also be present in 17-30 min wave records extracted from the continuous measurements and in statistical properties derived from them.
A long wave time series is needed to obtain reliable estimates of extreme values of sea surface elevation characteristics and of their probability of occurrence, which are of importance for applications.The uncertainty due to sampling variability will particularly affect the higher-order statistical moments like skewness and kurtosis, which are more unstable than the estimates of significant wave height and spectral/zero-crossing wave period.To obtain reliable estimates of these higher-order statistical moments, 250-350 repetitions of a 17 min wave record will often be required, as demonstrated by, for example, Bitner-Gregersen and Hagen (2003).Therefore numerical wave models remain an important supporting tool in an analysis of field data.
As pointed out by Tomita (2009), both the numerical nonlinear wave models and the wave spectral models can be utilized in research of extreme and rogue waves and their use is encouraged.The complimentary nature of these models is clear; they give different information about a sea state.A spectral wave model (e.g., the WAM model) provides sea state description only in a form of the two-dimensional wave spectrum but does not give any information about the instantaneous position of the sea surface in a given sea state.Note also that it accounts for wind forcing and resonant wave interactions but not for quasi-resonance interactions, which are responsible for occurrence of modulational instability and hence rogue waves (Onorato et al., 2013).Phase-resolving wave models, however, provide the water surface elevation from which statistical properties of individual waves can be extracted and include quasi-resonance interactions.Further, these nonlinear wave models allow simulating a wave record for the required time duration and, by repeating the 17-30 min simulations, significantly reducing the uncertainty due to sampling variability in estimated sea surface characteristics and their probability of occurrence.Although the spectral wave model as well as the nonlinear numerical wave model are computationally intense, the great advance in enhancing computer power has made the coupling between these models feasible.
In the present study we demonstrate the complementary nature of the wave spectral model WAM and the numerical nonlinear wave model based on the Euler equations and solved with the HOSM proposed in West et al. (1987).The coupling is applied to investigate statistical properties of surface oscillations during the particularly severe Andrea storm, which crossed the central part of the North Sea on 8-9 November 2007.During this storm, on 9 November 2007 a rogue wave called Andrea was recorded at the Ekofisk field (Magnusson and Donelan, 2013).This wave is comparable in characteristics, both with respect to the wave height and wave crest criterion, to the well-known New Year wave (called also the Draupner wave) measured by Statoil at the Draupner platform on 1 January 1995 (Haver and Anderson, 2000).
The paper is organized as follows: Sect. 2 describes the Andrea storm and the Andrea wave.Section 3 addresses hindcast data used in the analysis while Sect. 4 is dedicated to characteristics of the Andrea storm and comparison of numerical results to some characteristics of the Andrea wave recorded at the Ekofisk field.Conclusions are summarized in Sect. 5.

The Andrea storm and the Andrea wave
A low pressure area entered the northern North Sea on 8 November 2007.It covered southern Norway and moved in the morning of 9 November towards southern Sweden.Strong westerly winds (50-55 knots) followed the low pressure area and a high wave field (significant wave height of 10-11 m) was built up in the north area of the Ekofisk field (see Fig. 1) in the afternoon of 8 November 2007.The wind slightly decreased at Ekofisk around 18:00-21:00 UTC (universal time coordinated); from 22 to 19 m s −1 .The strongest wind field passed the northeast and east of Ekofisk on 9 November around 06:00 UTC and generated waves of up to 11-12 m (for details see Magnusson and Donelan, 2013).The Andrea wave was recorded at Ekofisk by the downlooking lasers just past 00:00 UTC on 9 November 2007.This wave is comparable in characteristics to the well-known New Year wave (the Draupner wave) measured by Statoil at the Draupner platform on 1 January 1995 (Haver and Anderson, 2000).The characteristics of the Andrea wave, as reported by Magnusson and Donelan (2013), are compared to the New Year wave ones in Table 1.H s denotes the significant wave height, T p the spectral peak period, C max the maximum crest height in the wave record, H max the maximum zero-downcrossing wave height with the crest C max , CF is the maximum crest factor (crest criterion), HF is the maximum height factor (height criterion).
CF > 1.3 (or > 1.2) and HF> 2 within a 20 min wave record represent simplified definitions of a rogue wave (see e.g.Bitner-Gregersen and Toffoli, 2012).If both criteria are fulfilled a rogue wave can be classified as a double rogue wave (Krogstad et al., 2008).As seen in Table 1 both the New Year wave as well as the Andrea wave can be called a double rogue wave.Note that both waves are recorded in the North Sea from the platforms located over a water depth of ca.75 m.

Hindcast data
The hindcast data used in this study were retrieved from the European Centre for Medium-Range Weather Forecast (ECMWF) archive.Wave parameters were acquired at the nearest grid point to Ekofisk.The data cover the period of the Andrea storm history from 00:00 UTC 8 November 2007 to 00:00 UTC 11 November 2011 and are stored every 6 h.The selected grid point is at a water depth of 74 m and within a distance of ca.50 km from the Ekofisk field.The data include the wind speed as well as the significant wave height and spectral peak period for the total sea, wind sea and swell.The components of the wind sea are those that are still under the influence of the local wind forcing and are detected as the part of wave spectrum where the wind input source term is positive.The remaining part of the wave spectrum is considered as swell (see, e.g., Hauser et al., 2005, for details).
Figure 2 shows a history of the total significant wave height, spectral period and related sea state steepness during the Andrea storm.The significant wave height reaches its maximum at 06:00 UTC on 9 November 2007.This is consistent with the findings of Magnusson and Donelan (2013) based on the NORA10 (Norwegian 10 km Reanalysis Archive) hindcast data developed at the Norwegian Meteorological Institute with major support from a consortium of oil companies (see e.g., Aarnes et al., 2011).The maximum wave height of 9.8 m is associated with the largest spectral peak period and the highest sea state steepness of 0.14 during the storm.It should be noted that the same high steepness is observed before the significant wave height reaches its maximum.Because the hindcast data are sampled every 6 h it is not possible to detect the steepness at 00:40 UTC on 9 November 2007, when the Andrea wave was recorded at the Ekofisk field.
The probability of occurrence of rogue waves is related to mechanisms generating them.It is interesting to note that the output from a wave spectral model can be utilized when indicating a mechanism responsible for the occurrence of a rogue wave (e.g.Tamura et al., 2009), even though it may not always allow reaching the firm conclusions.We illustrate this below for the Andrea wave recorded at Ekofisk during the Andrea storm.
These mechanisms have also been considered in order to indicate a possible phenomenon responsible for generating the Andrea wave.The linear focusing is occurring very seldom and because the present study is addressing nonlinear waves the linear focusing has been eliminated from further considerations.Further, no strong current has been reported in the Ekofisk area representing the intermediate water depth ocean zone.Therefore wave-current interactions and shallow water effects seem not to be responsible for the occurrence of the Andrea wave.The two remaining rogue wave generation mechanisms, namely crossing seas and quasi-resonance nonlinear interactions (modulational instability), could be regarded as the only possible candidates that generated the Andrea wave.In order to select one of them the time history of wind sea and swell during the Andrea storm has been studied.
Figure 3 shows evolution of significant wave height for wind sea and swell during the Andrea storm at the nearest grid point to Ekofisk.Wind sea clearly dominates the total sea during the growth, peak and decay of the storm.Most of the time the significant wave height of swell is only slightly above 1 m, much lower than the wind sea significant wave height, which mostly reaches 10 m at the peak of the storm.Therefore the total sea of the Andrea storm is dominated by wind sea.The wind sea and swell have approximately the same energy (the significant wave height around 1 m) only at the beginning and end of the storm.
It is well established that two wave trains with similar energy and frequencies traveling at particular angles can trigger modulational instability and be responsible for the formation of rogue waves (Onorato et al., 2006a(Onorato et al., , 2010)).Such results have been confirmed through recent numerical simulations of the Euler equations and experimental work performed in the MARINTEK Laboratories (Toffoli et al., 2011).The investigations have showed that the kurtosis, a measure of the probability of occurrence of extreme waves, depends on an angle β between the crossing wave systems.The maximum value of kurtosis is achieved for 40 < β < 60.No such conditions where wind sea and swell have the same energy and spectral peak frequency and are crossing each other under the angle 40 < β < 60 have been identified in the hindcast data from the Andrea storm.The ECMWF hindcast data seems to point out that the Andrea wave might have occurred in the sea state more prone to extreme waves as a result of modulational instability (i.e., the sea state with relatively high steepness and not a particularly broad spectrum).
This conclusion is supported by Figs. 4 and 5 presenting evolution of the directional wave spectrum during the Andrea storm at the location considered.One wave system is seen in the period from 18:00 UTC 8 November 2007 to 06:00 UTC 9 November 2007 within which the Andrea wave was recorded.

Numerical simulations
Numerical simulations have been carried out to get further insight into the Andrea storm characteristics.Short-term wave records at the sampling interval of 6 h have been generated by solving the Euler equations with the HOSM as proposed by West et al. (1987).
In the case of constant water depth (h = 74 m in this study), the velocity potential (x,z,t) of an irrotational, inviscid, and incompressible liquid satisfies the Laplace's equation everywhere in the fluid.The boundary conditions are such that the vertical velocity at the bottom (z = −74) is zero, and the kinematic and dynamic boundary conditions are satisfied for the velocity potential (x,y,t)= (x,y, η(x,y,t),t) on the free surface; that is, z = η(x,y,t) (see Zakharov, 1968).The expressions of the kinematic and dynamic boundary conditions are as follows: where the subscripts denote the partial derivatives, and W (x, y, t) = z |η represents the vertical velocity evaluated at the free surface.
The time evolution of an initial surface elevation can be calculated from Eqs. (1) and (2).For this study, we have used the HOSM, which was independently proposed by West et al. (1987) and Dommermuth and Yue (1987).A comparison of these two approaches (Clamond et al., 2006) has shown that the formulation proposed by Dommermuth and Yue (1987) is less consistent than the one proposed by West et al. (1987) as it does not converge when the amplitude is very small; the latter, therefore, has been applied herein.The advantage of HOSM in comparison to other methods is that it allows simulating a large number of random realizations of the surface elevation, within a reasonable computational time, without limitations in terms of the spectral bandwidth.
HOSM uses a series expansion in the wave slope of the vertical velocity W (x, y, t) about the free surface.In the present study we have considered a third-order expansion so that the four-wave interaction is included (see Tanaka, 2001Tanaka, , 2007)).Under these circumstances, the solution presented herein is not fully nonlinear.The expansion is then used to evaluate the velocity potential (x, y, t) and the surface elevation η(x, y, t) from Eqs. ( 1) and (2) at each instant of time.The time integration is performed by means of a fourthorder Runge-Kutta method with a time step t = T p /200 (T p is the spectral peak period).All aliasing errors generated in the nonlinear terms are removed (see West et al., 1987, andTanaka, 2001, for details).A small time step is used to  minimize the energy leakage.Throughout the simulations the variation of total energy remains lower than 0.5 %.
The model works under the assumption that the water depth is uniform.At the Ekofisk area, including the location considered, the variation of bottom topography is negligible and hence such an assumption does not affect the end result of the simulations.It is worth mentioning, however, that where bottom topography is changing wave dynamics could be affected and thus a variable bathymetry should be considered (e.g., the numerical model of Fructus and Grue, 2007).
The HOSM requires as input an initial sea surface and velocity potential with periodic boundary conditions.For the purpose of the present study, the initial conditions are extracted from the hindcast wave spectrum.Initially, the spectrum E(ω, θ ) is converted into a wavenumber spectrum E(k x , k y ) with the linear dispersive relation.An input surface is then extracted by means of an inverse Fast Fourier Transform with random amplitude and random phase approximation.Under these circumstances, the initial surface is linear (and hence normally distributed).The nonlinearity, both in terms of bound wave and free wave modes, is automatically built by the Euler equations (cf.Dommermuth 2000;Ducrozet et al., 2007;Toffoli et al., 2010b).The velocity potential is then calculated from the surface elevation in accordance with linear wave theory.
Note that the converted wavenumber spectrum does not ensure constant k x and k y , which are needed to ensure period boundaries.Furthermore, standard output spectra only contain a very limited amount of frequencies and the resulting surface would be too course for statistical analysis.Therefore, in order to force periodic boundaries and refine the surface, wavenumbers are redefined by subdividing the domain from 0 to k max into a much larger number of modes spaced at constant intervals.The concurrent values of energy are obtained by linear interpolation.Note that this procedure may induce a shift in the peak period if the number of modes is not large enough.This is due to the fact that the new wavenumbers may not necessarily coincide with the original modes.To overcome this difficulty, 256 modes were applied in the present analysis in order to maintain the peak period unchanged.
It is interesting to mention that Dommermuth (2000) pointed out that initializing a nonlinear model with a linear surface poses problems due to the generation of spurious modes, which are nonphysical.Both the bound and free waves are affected by generations of standing waves making the free surface rougher than it should be.Dommermuth (2000) proposed an adjustment procedure that eliminates the spurious generation of standing waves.The analysis of regular wave packets indicates that the initial growth rate of the Benjamin-Feir instability is delayed because of the presence of standing waves but there is no apparent effect on the growth rate.The effect on irregular wave fields, however, is not completely clear yet.Nonetheless, comparisons of sta- tistical properties of the surface elevation from HOSM simulations initialized with linear surface (thus no adjustments) and laboratory experiments in directional wave basins (see, e.g., Toffoli et al., 2010aToffoli et al., , 2013, for infinite and finite water depth, respectively) showed a very good agreement both in terms of spatial/temporal evolution and maximum values of statistical moments.
In the definition of the initial surface, x and y were selected in order to ensure 10 wavelengths in the physical domain.The initial surface has been let to evolve according to Eqs. ( 1) and (2) for 70 peak periods.According to laboratory experiments (see e.g., Toffoli et al., 2010a), this is sufficient to ensure the full development of modulational instability (the mechanisms responsible for the formation of rogue waves), during typical stormy conditions.More generally, the timescale for the development of modulational instability and hence rogue waves is proportional to the (ka) −2 , where k is the dominant wave number and a is half the significant wave height.In order to have statistically significant results, 150 repetitions of the wave surface with different sets of random amplitudes and phases (see e.g., Toffoli et al., 2008Toffoli et al., , 2010a) have been simulated.
Figures 4 and 5 show evolution of the directional wave spectrum in the period 00:00 UTC 8 November 2007-18:00 UTC 11 November 2007 covering the time history of the Andrea storm.These spectra have been transformed to the wavenumber spectra E(k x , k y ) and used as an input to the numerical HOSM simulations.
Time histories of the crest ratio CF = C max /H s , skewness and kurtosis in the period 8-11 November 2007 (the Andrea storm time history) derived from the HOSM simulations are shown in Fig. 6 (A max denotes the maximum amplitude equal to C max ).The HOSM model predicts the maximum crest ratio C max / H s = 1.4 and the corresponding kurtosis of 3.35.Note that the crest ratio for the Andrea wave is higher, C max / H s = 1.63 (see Table 1), but the difference between the two ratios could be expected as the location considered is in the distance of 50 km from the Ekofisk field.
Temporal evolution of skewness, kurtosis and A max /A 0 (A 0 denotes the initial amplitude used in the simulations) within the 70T p time period obtained from the HOSM simulations of the peak time of the Andrea storm (18:00 UTC 8 November 2007, 00:00 UTC 9 November 2007, and 06:00 UTC 9 November 2007) are shown in Figs. 7, 8 and 9, respectively.In the location considered, the maximum kurtosis is observed already after 30 T p s simulations at 18:00 UTC on 8 November 2007.The results presented in Figs.6-9 point out that coupling of the wave spectral model and the nonlinear phase-resolving model could predict the occurrence of the Andrea wave.
It is interesting to note that a sea state capable of triggering rogue waves occurred several times during the Andrea storm history at the considered location.North Sea scatter diagrams of H s and T p indicate that sea states with similar steepness may occur frequently in some locations of the North Sea.Scatter diagrams, however, do not give any information about the frequency-directional wave spectrum and therefore firm conclusions regarding occurrence of rogue-prone sea states cannot be reached by studying scatter diagrams alone.

Conclusions
The study shows how the wave spectral WAM model and the HOSM model can be coupled to forecast/hindcast the occur- The considered location is approximately 50 km from the Ekofisk field, where the Andrea wave was recorded and therefore the wave characteristics derived from the numerical simulations cannot be uncritically compared with those of the Andrea wave.However, the results point out that coupling of a spectral wave model with a nonlinear phase-resolving model is able to predict the occurrence of an extreme wave event.
A spectral model coupled with the HOSM model provides statistical information about waves based on the actual hindcast/forecast spectrum, whether this is bimodal or unimodal.The coupled model's output allows deriving, among other wave properties, a distribution of sea surface elevation, the maximum wave crest, skewness and kurtosis.The proposed approach, although demonstrated for the Andrea storm, is of general character and can be applied to investigations of any other storm.Further, it represents a good supporting tool for an analysis of field data.
It needs to be noted that HOSM cannot be applied for prediction of rogue waves in very steep sea states because it does not include wave breaking.Although a wave spectral model like WAM provides valuable input to HOSM simulations the commonly used 6 h sampling interval may present a limitation when studying extreme and rogue waves; important information between sampling intervals may be missed.
The analysis shows that when the Andrea storm is passing the North Sea rogue waves can be expected in several locations, not only at Ekofisk where the Andrea wave was recorded.Rogue waves are produced during a storm's development; when a storm builds up and when it weakens, in the location considered by the study.They are not observed when the storm reaches the largest H s .
Uncertainties associated with the WAM and HOSM model, which will affect the presented results, have not been investigated in this study.Further investigations are still called for to evaluate impact of these uncertainties on extreme and rogue wave predictions.
The approach for coupling the wave spectral model with the nonlinear phase-resolving model presented herein can be considered for use with forecasting purposes.Although these models are computationally intense, the great advance in enhancing computer power has made the coupling between them feasible.

Figure 1 .
Figure 1.Location of the Ekofisk field.

Figure 2 .
Figure 2. History of significant wave height, spectral wave period and sea state steepness for the total sea during the Andrea storm.

Figure 3 .
Figure 3. History of significant wave height for wind sea and swell during the Andrea storm.

Figure 6 .
Figure 6.Time histories of the crest ratio CF= C max / H s , skewness and kurtosis in the period 8-11 November 2007 (A max is the maximum amplitude equal to C max ).

Figure 7 .
Figure 7. Temporal evolution of skewness, kurtosis and A max /A 0 (A 0 denotes the initial amplitude used in the simulations while A max is the maximum amplitude (crest)), 8 November 2007, 18:00 UTC.

Figure 8 .
Figure 8. Temporal evolution of skewness, kurtosis and A max /A 0 (A 0 denotes the initial amplitude used in the simulations while A max is the maximum amplitude (crest)), 9 November 2007, 00:00 UTC.

Figure 9 .
Figure 9. Temporal evolution of skewness, kurtosis and A max /A 0 (A 0 denotes the initial amplitude used in the simulations while A max is the maximum amplitude (crest)), 9 November 2007, 06:00 UTC.

Table 1 .
Characteristics of the Andrea and New Year waves.