Accurate and Fast Generation of Irregular Short Crested Waves by Using Periodic Boundaries in a Mild-Slope Wave Model

: In this work, periodic lateral boundaries are developed in a time dependent mild-slope equation model, MILDwave, for the accurate generation of regular waves and irregular long and short crested waves in any direction. A single wave generation line inside the computational domain is combined with periodic lateral boundaries. This generation layout yields a homogeneous and thus accurate wave ﬁeld in the whole domain in contrast to an L-shaped and an arc-shaped wave generation layout where wave diffraction patterns appear inside the computational domain as a result of the intersection of the two wave generation lines and the interaction with the lateral sponge layers. In addition, the performance of the periodic boundaries was evaluated for two different wave synthesis methods for short crested waves generation, a method proposed by Miles and a method proposed by Sand and Mynett. The results show that the MILDwave model with the addition of periodic boundaries and the Sand and Mynett method is capable of reproducing a homogeneous wave ﬁeld as well as the target frequency spectrum and the target directional spectrum with a low computational cost. The overall performance of the developed model is validated with experimental results for the case of wave transformation over an elliptic shoal (Vincent and Briggs shoal experiment). The numerical results show very good agreement with the experimental data. The proposed generation layout using periodic lateral boundaries makes the mild-slope wave model, MILDwave, an essential tool to study coastal areas and wave energy converter (WEC) farms under realistic 3D wave conditions, due to its signiﬁcantly small computational cost and its high numerical stability and robustness. In Figure 8, a comparison is made between the target frequency spectrum (S t ) and the simulated frequency spectra (S 0 and S 45 ) for the cases of θ = 0 ◦ and θ = 45 ◦ . The surface elevations, η , at the electronic wave gauges, which are positioned at the centre of the computational domain, are recorded from t = 40 T p to t = 600 T p with a sampling interval of 0.2 s. The recorded data are processed in segments of 2048 points per segment. A taper window and an overlap of 20% are applied for smoother and statistically more signiﬁcant spectral estimates. The resulting frequency spectra agree well with S t apart from those that correspond to wave frequencies higher than 2 f p , since there is no energy for these frequencies in the MILDwave model. Furthermore, it is observed that the periodic lateral boundaries do not affect the waves since for both wave propagation angles, the frequency spectra are similar.


Introduction
Numerical wave propagation models are commonly used as engineering tools for the study of wave transformation in coastal areas. Berkhoff [1] derived the first elliptic mild-slope wave equation and based on this, the parabolic model [2] and the hyperbolic model [3] have been developed to predict the transformation processes of regular waves, such as wave refraction, diffraction, shoaling, and reflection. In the parabolic model, wave reflection and diffraction in the direction of wave propagation are neglected and hence it suffers from low accuracy in cases where these phenomena are significant. On the other hand, the hyperbolic model, which takes into consideration wave reflection and diffraction in the direction of wave propagation, provides higher accuracy, but more computational time compared to the parabolic model. Moreover, to study the transformation of random waves, time dependent mild-slope equations have been developed. Radder and Dingemans [4] suggested a set of canonical equations, which are based on the time dependent mild slope equations and are derived using the Hamiltonian theory of surface waves. Booij [5] proved that these equations are valid for sea bottom slopes up to 1/3. However, Suh et al. [6] extended the latter equations by including higher order bottom effect terms proportional to the square of the bottom slope and to the bottom curvature to study wave propagation on rapidly varying topography.
For numerical prediction of the wave field in a coastal area, waves are generated along the offshore boundary of the computational domain and propagate towards the coastline. However, to be able to apply a sponge layer to absorb waves reflected towards the wave generation boundary, waves should be generated inside the numerical domain and not along the boundary. This internal wave generation technique in combination with numerical wave absorbing sponge layers was firstly proposed by [7] for Peregrine's [8] classical Boussinesq equations. Later, Lee and Suh [9] achieved wave generation for the mild slope equations of [3] and [4] by applying the source term addition method. The main observation was that the energy of the incident waves can be properly obtained from the viewpoint of energy transport. To generate multidirectional waves, they applied an L-shaped wave generator, which is generally composed of two wave generation lines; one parallel to the x-axis and one parallel to the y-axis of the numerical domain as well as wave absorbing sponge layers behind the two wave generation lines. However, wave diffraction patterns appear inside the computational domain as a result of the intersection of the two wave generation lines and due to the interaction with the lateral sponge layers. To deal with this problem, Lee and Yoon [10] proposed an arc-shaped wave generation line; two parallel lines connected to a semicircle to avoid wave diffraction caused at the intersection of the previous wave generation lines. Further, Kim and Lee [11] used an arc-shaped source band that gives smaller errors than the Lee and Yoon [10] method, especially for a coarse grid size. Recently, Lin and Yu [12] proposed a promising method for non-reflective boundaries in a mild-slope wave model to avoid the use of sponge layers. However, in non-reflective boundaries, the level of re-reflection strongly depends on the initial approximations since the characteristics of the reflected waves (i.e., wave angle, wave celerity) inside the numerical domain cannot be estimated a priori.
In the present paper, a wave generation layout using periodic lateral boundaries is developed where a single internal wave generation line parallel to the y-axis is combined with periodic lateral boundaries at the top and bottom of the domain. With this technique, the information leaving one end of the numerical domain enters the opposite end and thus no lateral sponges are required. In this way, the wave diffraction patterns that appear inside the computational domain as a result of the intersection of the two wave generation lines and due to the interaction with the lateral sponge layers (see Figure 1) are avoided. In Figure 1, the water surface elevation, η, is presented for a regular wave field generated by an L-shaped wave generator (red dashed lines) at an angle, θ = 45 • , from the x-axis. The absence of uniformity of the surface elevation along the crests and troughs, which can be observed in the same figure, is a sign of the existence of diffraction patterns inside the numerical domain.
Periodic boundaries are implemented in a time dependent mild-slope equation model, MILDwave, developed at Ghent University [13] in order to accurately generate regular and irregular waves in any direction. In addition, the performance of the periodic boundaries is evaluated for two different wave synthesis methods to generate short crested waves, a method proposed by Miles [14] and a method proposed by Sand and Mynett [15]. The mild-slope equations of Radder and Dingemans [4] without the extension of Suh et al. [6] are the basic equations employed in the phase-resolving model, MILDwave, used to simulate wave transformation processes, such as refraction, shoaling, reflection, transmission, and diffraction, intrinsically [16]. MILDwave has previously been used to predict wave diffraction and wave penetration inside harbours [17] and to study wave power conversion applications [18][19][20]. One of the challenges in the field of renewable energies is to determine the optimal geometrical layout for wave energy converter (WEC) farms, targeting the maximum possible energy production and the correct assessment of the impact of WEC farms on the wave field. To do so, accurate and detailed numerical modelling of WEC farms under realistic 3D wave conditions is considered crucial. This kind of application requires a homogeneous wave field in the whole numerical domain. To the present authors' knowledge, periodic boundaries, which are commonly used in non-hydrostatic models [21] and Boussinesq models [22], have not been used before in a mild-slope wave model to study short crested waves. So, the novelty of the present work concerns the capability of accurate and fast generation of such homogeneous wave fields with periodic boundaries in MILDwave. In engineering applications, where non-linearities are not significant, mild-slope wave models are preferred instead of Boussinesq models or non-hydrostatic models due to their significantly smaller computational cost and their high numerical stability and robustness, and thus further development of these models should be encouraged. The implementation of periodic boundaries is important as it introduces noteworthy improvements in mild-slope models, which can then make full use of their benefits for the study of WEC farms under oblique long crested regular and irregular waves or short crested waves (real sea waves). Periodic boundaries are implemented in a time dependent mild-slope equation model, MILDwave, developed at Ghent University [13] in order to accurately generate regular and irregular waves in any direction. In addition, the performance of the periodic boundaries is evaluated for two different wave synthesis methods to generate short crested waves, a method proposed by Miles [14] and a method proposed by Sand and Mynett [15]. The mild-slope equations of Radder and Dingemans [4] without the extension of Suh et al. [6] are the basic equations employed in the phaseresolving model, MILDwave, used to simulate wave transformation processes, such as refraction, shoaling, reflection, transmission, and diffraction, intrinsically [16]. MILDwave has previously been used to predict wave diffraction and wave penetration inside harbours [17] and to study wave power conversion applications [18][19][20]. One of the challenges in the field of renewable energies is to determine the optimal geometrical layout for wave energy converter (WEC) farms, targeting the maximum possible energy production and the correct assessment of the impact of WEC farms on the wave field. To do so, accurate and detailed numerical modelling of WEC farms under realistic 3D wave conditions is considered crucial. This kind of application requires a homogeneous wave field in the whole numerical domain. To the present authors' knowledge, periodic boundaries, which are commonly used in non-hydrostatic models [21] and Boussinesq models [22], have not been used before in a mild-slope wave model to study short crested waves. So, the novelty of the present work concerns the capability of accurate and fast generation of such homogeneous wave fields with periodic boundaries in MILDwave. In engineering applications, where non-linearities are not significant, mild-slope wave models are preferred instead of Boussinesq models or non-hydrostatic In the next section, a short description of the mild-slope wave propagation model, MILDwave, and the methodology of the implemented periodic boundaries are presented. Section 3 provides a detailed overview of the model results for regular and irregular long and short crested waves. In Section 4, the validation results are presented where the accuracy of the developed model is compared with experimental data. The last section provides conclusions and a summary discussion of the present study.

The Mild-Slope Wave Propagation Model, MILDwave
MILDwave is a mild-slope wave propagation model based on the depth integrated time dependent mild-slope equations of Radder and Dingemans [4]. These are given in Equations (1) and (2) and where η(x, y, t) and φ(x, y, t) are, respectively, the surface elevation and the velocity potential at the still water level, ∇ is the horizontal gradient operator, t is the time, and g is the gravitational acceleration. The coefficients, A and B, are calculated using Equations (3) and (4), respectively: with C the phase velocity and C g is the group velocity for a wave with a wave number, k, angular frequency, ω, wavelength, L w , and frequency, f. For irregular waves, C, C g , k, and ω are replaced by the wave characteristics for the carrier frequency, f. A finite difference scheme, as described in [23], is used to discretize and solve Equations (1) and (2). The domain is divided in grid cells with dimensions, ∆x and ∆y, and central differences are used both for spatial and time derivatives. Both η and φ are calculated at the center of each grid cell at different time levels, (n + 1/2)∆t and (n + 1)∆t as: where A and B are given by Equations (3) and (4) and the superscripts and subscripts stand for the relevant time step and the relevant cell of the grid, respectively. In the current configuration of MILDwave, the boundary conditions are formed in such a way that one layer of ghost cells is considered at each boundary. The values at the cells closest to the boundary are copied to the ghost cells and thus the layer of the ghost cells acts as a fully reflective boundary (solid wall). However, the effect of the fully reflective domain boundaries is negligible because absorbing sponge layers are applied at the outside boundaries to dissipate the incoming wave energy.
The grid cell size, ∆x = ∆y, is chosen so that L w /20 ≤ ∆x = ∆y ≤ L w /10 (for irregular waves, L w = shortest wave length (maximum wave frequency)) while the time step meets the Courant-Friedrichs-Lewy criterion to guarantee a stable and consistent result.
In MILDwave, waves are generated along a wave generation line near the offshore boundary by applying the source term addition method proposed by [9] where the source term propagates with the energy velocity. According to this method, additional surface elevation, η * , is added with the desired energy to the calculated surface elevation, η, at the wave generation line for each time step and is given by: Energies 2019, 12, 785

of 21
where ∆x is the grid size in the x-axis, ∆t is the time step, θ is the wave propagation angle with respect to the x-axis, η I the water surface elevation of incident waves, and C e is the energy velocity given by Equation (8): where the bar indicates that the variable is associated with the carrier angular frequency, ω.
It has been proven that the model of [4] can be used to simulate the transformation of long and short crested waves. To generate long crested irregular waves, a parameterized JONSWAP (Joint North Sea Wave Observation Project) spectrum, S(f) (Equation (9)), or a TMA (Texel, Marsen, Arsloe) spectrum (Equation (10)), which can be applied in shallow water, has been used as an input spectrum: where H s is the significant wave height, f p is the peak wave frequency, γ is the peak enhancement factor, α is the Phillips constant, and σ is the spectral width parameter, which depends on the value of the wave frequency: The frequency dependent factor, φ(f, d), which takes into account the effect of a finite water depth, d, is given by: To generate short crested waves, two different wave synthesis methods are employed in the present study, a single summation method proposed by [14] and a second single summation method proposed by [15]. In single summation models, each wave component must have a unique frequency and each frequency component can only travel in one direction, while several wave components are travelling in the same direction. According to the Miles method, the surface elevation is defined as follows: with the wave amplitude, A nm = 2S(f nm )D(f nm , θ m )M∆f∆θ, the wave angular frequency, ω nm = 2πf nm = 2π(M(n − 1) + m)∆f + 2πf min , the frequency interval, ∆f = f max −f min NM−1 , the wave propagation angle, θ m = (m − 1)∆θ + θ 0 − θ max , the wave propagation angle interval, ∆θ = 2θ max M−1 , the random phase, ε nm , and the maximum discrete wave direction, θ max . In this way, the directional spreading function is discretized into M equally spaced wave directions ranging from θ min to θ max .
Sand and Mynett [15] proposed a method in which the directional spectrum is decomposed as follows: η(x, y, t) = N ∑ n=1 2S(f n )∆f cos(ω n t − k n (x cos θ n + y sin θ n ) + ε n ) (14) In this method, the wave propagation angles are selected at random according to the cumulative distribution function of the directional spreading function, D(f, θ), and are assigned to each frequency component. The Miles method [14] provides an accurate representation of the targeted spreading function shape, but introduces different localized distortions in the frequency spectrum to obtain a Energies 2019, 12, 785 6 of 21 close fit to the spreading function. On the other hand, the Sand and Mynett method [15] yields an exact match to the frequency spectrum.
Several semi-empirical proposed formulations of the directional spreading function, D(f, θ), have been reported, most of which consider the spreading function to be independent of the wave frequency.
Here, an alternative of the well-known spreading function of [24] is employed [25]: where s 1 is the directional spreading parameter, Γ is the Gamma function, and θ 0 is the wave propagation angle.
To derive the relation between the directional spreading parameter, s 1 , and the spreading standard deviation, σ θ , which is called the directional width, the methodology of [26] is followed: Substitution of Equation (15) into Equation (16) yields: The relation between s 1 and σ θ is indicated in Figure 2. Hence, fixed values of s 1 and σ θ can be determined for wind and swell waves: Energies 2019, 12, 785 7 of 23 The relation between s and σ is indicated in Figure 2. Hence, fixed values of s and σ can be determined for wind and swell waves: Figure 2. Relation between the directional spreading parameter, s , and the spreading standard deviation, σ .

Periodic Boundaries
To create a homogeneous wave field of oblique long crested regular and irregular waves or short crested waves, periodic lateral boundaries have been implemented. In the present study, a single wave generation line parallel to the y-axis is combined with periodic lateral boundaries at the top and bottom of the domain. In this way, the information leaving one end of the numerical domain enters the opposite end and thus the required model length in this direction is reduced. This technique can lead to a more homogeneous wave field than an L-shaped and an arc-shaped wave generation since no wave diffraction problems are caused by the presence of lateral sponge layers and the intersection of the two wave generation lines (Figure 1). Figure 3 is a schematic representation of the way that the periodic boundaries are working. Firstly, the left part of the figure presents the numerical domain under investigation, with the squares Figure 2. Relation between the directional spreading parameter, s 1 , and the spreading standard deviation, σ θ .

Periodic Boundaries
To create a homogeneous wave field of oblique long crested regular and irregular waves or short crested waves, periodic lateral boundaries have been implemented. In the present study, a single wave generation line parallel to the y-axis is combined with periodic lateral boundaries at the top and bottom of the domain. In this way, the information leaving one end of the numerical domain enters the opposite end and thus the required model length in this direction is reduced. This technique can lead to a more homogeneous wave field than an L-shaped and an arc-shaped wave generation since no wave diffraction problems are caused by the presence of lateral sponge layers and the intersection of the two wave generation lines ( Figure 1). Figure 3 is a schematic representation of the way that the periodic boundaries are working. Firstly, the left part of the figure presents the numerical domain under investigation, with the squares representing a random number of cells where Nx and Ny is the number of cells in the x and y direction, respectively. A layer of ghost cells is present next to each vertical boundary, acting as a fully reflective wall as described in the previous section. On the other hand, the dashed lines, which are parallel to the x-axis, represent the periodic boundaries.

Results by the Model, MILDwave
In this section, the applicability and accuracy of periodic boundaries to propagate regular and irregular long and short crested waves over a flat bottom is examined. For the evaluation of the numerical simulation results, a common method is followed. The results for a series of test cases are presented in terms of a disturbance coefficient, k , for the entire numerical domain. For regular waves, k is given by the ratio, H H ⁄ , where H is the local wave height and H is the wave height at the wave generation boundary. In the case of irregular waves, the significant wave heights are used to calculate the k value (H H ⁄ ). For all test cases, the numerical basin is 20L long and 20L wide (where L is the wave length and is equal to L = 100 m) while a uniform water depth, d = 7.5 m, is used. Thus, in the following numerical tests, the target wave field is a homogeneous one.
A homogeneous wave field is characterized by the lack of significant fluctuations from the target, k , value, which is equal to k = 1.

Generation of Regular Waves
Firstly, oblique regular waves are examined to compare three different wave generation layouts. The first wave generation layout ( Figure 4a) is an L-shaped wave generator, which is composed by two orthogonal lines as proposed by [9]. The two wave generation lines intersect at a point outside the sponge layers and they extend to the computational domain boundary. The hatched areas indicate the areas of the numerical domain covered by wave absorbing sponge layers. These sponge layers are necessary to absorb waves reflected by the numerical domain boundaries and by structures. The second wave generation layout ( Figure 4b) is based on the method suggested by [10], who proposed an arc-shaped wave generation layout where two parallel lines are connected to a semicircle. Figure  4c shows the third wave generation layout that is the main research object of the present paper, in which a single wave generation line parallel to the y-axis is combined with periodic lateral boundaries at the top and bottom of the domain. This means that the computational grid is repeated in the ydirection and hence waves are passing freely from the one periodic boundary to the other. At the position of the periodic boundaries, Equation (5) is solved considering that the two boundaries are adjacent, as it can be observed from the right part of Figure 3, yielding Equations (19) and (20) instead of Equation (5) for the bottom (j = 1) and top (j = Ny) layer of cells, respectively: This implementation is valid under the assumption that the bathymetry close to the lateral boundaries is uniform for the waves to pass from one boundary to the other without any distortion and thus the continuity to be satisfied. In addition, the distance between the two boundaries should be an integer multiple of the wavelength to avoid phase differences. To ensure that this requirement is valid, the wave propagation angle, θ, of each wave component is slightly adjusted such that Equation (21) is satisfied [27]: where k is the wave number of each wave component, α is an integer, and L per (L per = ∆y * Ny) is the periodicity length. This means that the simulated wave propagation angle slightly deviates from the input wave angle. To minimize this deviation, the periodicity length, L per , of the computational domain should be carefully chosen to be larger than one periodic wave length, L y (L y = L w / sin θ) [28]. Hence, in the case of irregular long crested waves, each wave component propagates in a slightly different direction than the mean value and a directional spread of the wave energy occurs. However, it has been verified that the spreading standard deviation, σ θ , is very small (σ θ < 2 • ) and the waves can be considered as uni-directional.

Results by the Model, MILDwave
In this section, the applicability and accuracy of periodic boundaries to propagate regular and irregular long and short crested waves over a flat bottom is examined. For the evaluation of the numerical simulation results, a common method is followed. The results for a series of test cases are presented in terms of a disturbance coefficient, k d , for the entire numerical domain. For regular waves, k d is given by the ratio, H/H GB , where H is the local wave height and H GB is the wave height at the wave generation boundary. In the case of irregular waves, the significant wave heights are used to calculate the k d value (H s /H sGB ). For all test cases, the numerical basin is 20 L w long and 20 L w wide (where L w is the wave length and is equal to L w = 100 m) while a uniform water depth, d = 7.5 m, is used. Thus, in the following numerical tests, the target wave field is a homogeneous one. A homogeneous wave field is characterized by the lack of significant fluctuations from the target, k d , value, which is equal to k d t = 1.

Generation of Regular Waves
Firstly, oblique regular waves are examined to compare three different wave generation layouts. The first wave generation layout (Figure 4a) is an L-shaped wave generator, which is composed by two orthogonal lines as proposed by [9]. The two wave generation lines intersect at a point outside the sponge layers and they extend to the computational domain boundary. The hatched areas indicate the areas of the numerical domain covered by wave absorbing sponge layers. These sponge layers are necessary to absorb waves reflected by the numerical domain boundaries and by structures. The second wave generation layout (Figure 4b) is based on the method suggested by [10], who proposed an arc-shaped wave generation layout where two parallel lines are connected to a semicircle. Figure 4c shows the third wave generation layout that is the main research object of the present paper, in which a single wave generation line parallel to the y-axis is combined with periodic lateral boundaries at the top and bottom of the domain. This means that the computational grid is repeated in the y-direction and hence waves are passing freely from the one periodic boundary to the other.
Due to the nature of the governing equations of the numerical model, MILDwave (linear mild-slope equations), and the fact that a flat empty numerical basin is examined, the performance of the model for different wave heights and periods is similar. Thus, regular waves with a wave period of T = 12 s, wave height of H = 1 m, and wave propagation angles from θ = 0 • to θ = 45 • with respect to the x-direction are simulated. The numerical basin consists of an inner computational domain without the sponge layers with dimensions of 10 L w × 10 L w for the first two layouts and 10 L w × 20 L w for Energies 2019, 12, 785 9 of 21 the third one, as well as sponge layers with a width of 5 L w at the edges of the wave basin. The sponge layer formula as described in [18] is used to minimize reflections, where the free surface elevation is multiplied on each new time step with an absorption function that has the value of 1 at the start of the sponge layer and smoothly decreases until it reaches the value of 0 at the end. As a result, the reflection is zero. The grid cell size is chosen so that ∆x = ∆y = L w /20 = 5 m. To obtain a steady state wave field, waves are generated for a duration of 1000 s with a time step of ∆t = 0.25 s.   Figure 5 shows the plan views of the resulting disturbance coefficient, k , and the water surface elevation, η, in the inner computational domain at a time of t = 80T, where T is the wave period, for each wave generation layout for a wave propagation angle of θ = 30°. The black solid contour lines indicate the region where the surface elevation is out of the target range (η > 0.5 m, η < −0.5 m). It is observed that the third wave generation layout with the periodic lateral boundaries yields a more homogeneous wave field across the whole domain, which is not disturbed by unwanted wave diffraction patterns in contrast to the other two wave generation layouts. These wave diffraction disturbances are caused by the lateral sponge layers and the intersection of the two wave generation lines, and thus the k values fluctuate around the target value (k = 1).  Figure 5 shows the plan views of the resulting disturbance coefficient, k d , and the water surface elevation, η, in the inner computational domain at a time of t = 80 T, where T is the wave period, for each wave generation layout for a wave propagation angle of θ = 30 • . The black solid contour lines indicate the region where the surface elevation is out of the target range (η > 0.5 m, η < −0.5 m). It is observed that the third wave generation layout with the periodic lateral boundaries yields a more homogeneous wave field across the whole domain, which is not disturbed by unwanted wave diffraction patterns in contrast to the other two wave generation layouts. These wave diffraction disturbances are caused by the lateral sponge layers and the intersection of the two wave generation lines, and thus the k d values fluctuate around the target value (k d t = 1). To obtain a quantified difference between the three different wave generation layouts, the absolute error is calculated in a specific test area far from the sponge layers as in [10]. The error is determined as follows: To obtain a quantified difference between the three different wave generation layouts, the absolute error is calculated in a specific test area far from the sponge layers as in [10]. The error is determined as follows: where H c is the local resulting wave height at a point of the numerical domain and H t is the target wave height. The test area in which the absolute error is calculated is 6 L w long and 6 L w wide and is indicated in Figure 4 (area of crosses). In Figure 6, the mean (E mean ) and the maximum (E max ) value of the calculated absolute errors for each wave generation layout for a range of different wave propagation angles are plotted. It is observed, as seen from Figure 5 for θ = 30 • , that the periodic boundaries' wave generation layout provides the smallest value of E mean and E max for all θ values. For the L-shaped wave generation layout and the periodic boundaries' wave generation layout, E mean is the lowest for θ = 45 • with values of 0.042 and 0.007, respectively, while for the arc-shaped wave generation layout, E mean is the lowest for θ = 30 • with a value of 0.021. The trend of E max (Figure 6b) is similar to that of the values of E mean from Figure 6a. The E max values of the oblique regular wave field generated using periodic boundaries are significantly smaller than the resulting E max values for the L-shaped and the arc-shaped layouts for all wave directions. In addition, the comparison between the L-shaped wave generation layout and the arc-shaped wave generation layout reveals that in the case of the latter, the wave diffraction due to the intersection of the two generation lines is avoided [10]. However, for wave directions of θ larger than 30 • , E max is similar for both layouts with values around 0.105 while for the periodic boundaries' wave generation layout, it is only 0.012. Hence, this proves that a single wave generation line combined with periodic lateral boundaries provides results of a higher accuracy compared to the other two wave generation layouts and is capable of providing a homogeneous wave field of short crested waves with minimal error. where H is the local resulting wave height at a point of the numerical domain and H is the target wave height. The test area in which the absolute error is calculated is 6 L long and 6 L wide and is indicated in Figure 4 (area of crosses). In Figure 6, the mean (E ) and the maximum (E ) value of the calculated absolute errors for each wave generation layout for a range of different wave propagation angles are plotted. It is observed, as seen from Figure 5 for θ = 30°, that the periodic boundaries' wave generation layout provides the smallest value of E and E for all θ values. For the L-shaped wave generation layout and the periodic boundaries' wave generation layout, E is the lowest for θ = 45° with values of 0.042 and 0.007, respectively, while for the arcshaped wave generation layout, E is the lowest for θ = 30° with a value of 0.021. The trend of E (Figure 6b) is similar to that of the values of E from Figure 6a. The E values of the oblique regular wave field generated using periodic boundaries are significantly smaller than the resulting E values for the L-shaped and the arc-shaped layouts for all wave directions. In addition, the comparison between the L-shaped wave generation layout and the arc-shaped wave generation layout reveals that in the case of the latter, the wave diffraction due to the intersection of the two generation lines is avoided [10]. However, for wave directions of θ larger than 30°, E is similar for both layouts with values around 0.105 while for the periodic boundaries' wave generation layout, it is only 0.012. Hence, this proves that a single wave generation line combined with periodic lateral boundaries provides results of a higher accuracy compared to the other two wave generation layouts and is capable of providing a homogeneous wave field of short crested waves with minimal error.

Generation of Irregular Long Crested Waves
We consider a test case of uni-directional irregular waves fitting a JONSWAP spectrum, with a peak wave period, T p = 12 s, a significant wave height, H s = 1 m, and a peak enhancement factor, γ = 3.3. The frequency range is confined between 0.75 f p and 2 f p , which covers 94% of the total energy, to ensure numerical stability. Periodic lateral boundaries are applied at the top and bottom of the domain similar to Section 3.1, while the length of the sponge layers at both ends of the numerical basin is determined by the longest wave length in the incident wave spectrum. The grid cell size and the time step, ∆t, are defined by considering the highest frequency (2 f p ) and the lowest frequency (0.75 f p ), respectively.
A single wave generation line is placed at the left boundary after the sponge layer. To obtain a steady state wave field, time series of the surface elevation of irregular waves are generated for a duration of 7200 s with a time step of ∆t = 0.2 s. Figure 7 shows the plan views of the k d values and the water surface elevations in an inner domain (part of the domain without the left and right sponge layers) of 1 km in both directions, at a time instant of t = 7000 s, for wave propagation angles of θ = 0 • and θ = 45 • . It is clear that for both wave propagation angles, the k d is uniform all over the inner domain with values close to k d t = 1. In addition, it is observed by comparing the surface elevation of the two different wave directions that for the case of θ = 45 • , the surface is not as uniform along the crest as it is for the case of θ = 0 • . This is happening due to the fact that in order to ensure that waves are periodic, the wave propagation angle, θ, of each wave component is slightly adjusted and a directional spread of the wave energy occurs as described in detail in Section 2.2. Hence, the wave propagation angle varies from 43.2 • to 47.7 • for the present wave conditions. However, the deviation from the mean wave direction is small and the waves can be considered as uni-directional.

Generation of Irregular Long Crested Waves
We consider a test case of uni-directional irregular waves fitting a JONSWAP spectrum, with a peak wave period, T = 12 s, a significant wave height, H = 1 m, and a peak enhancement factor, γ = 3.3. The frequency range is confined between 0.75 f and 2 f , which covers 94% of the total energy, to ensure numerical stability. Periodic lateral boundaries are applied at the top and bottom of the domain similar to Section 3.1, while the length of the sponge layers at both ends of the numerical basin is determined by the longest wave length in the incident wave spectrum. The grid cell size and the time step, Δt, are defined by considering the highest frequency (2f ) and the lowest frequency (0.75f ), respectively.
A single wave generation line is placed at the left boundary after the sponge layer. To obtain a steady state wave field, time series of the surface elevation of irregular waves are generated for a duration of 7200 s with a time step of Δt = 0.2 s. Figure 7 shows the plan views of the k values and the water surface elevations in an inner domain (part of the domain without the left and right sponge layers) of 1 km in both directions, at a time instant of t = 7000 s, for wave propagation angles of θ = 0° and θ = 45°. It is clear that for both wave propagation angles, the k is uniform all over the inner domain with values close to k = 1. In addition, it is observed by comparing the surface elevation of the two different wave directions that for the case of θ = 45°, the surface is not as uniform along the crest as it is for the case of θ = 0°. This is happening due to the fact that in order to ensure that waves are periodic, the wave propagation angle, θ, of each wave component is slightly adjusted and a directional spread of the wave energy occurs as described in detail in Section 2.2. Hence, the wave propagation angle varies from 43.2° to 47.7° for the present wave conditions. However, the deviation from the mean wave direction is small and the waves can be considered as uni-directional.  In Figure 8, a comparison is made between the target frequency spectrum (S t ) and the simulated frequency spectra (S 0 and S 45 ) for the cases of θ = 0 • and θ = 45 • . The surface elevations, η, at the electronic wave gauges, which are positioned at the centre of the computational domain, are recorded from t = 40 T p to t = 600 T p with a sampling interval of 0.2 s. The recorded data are processed in segments of 2048 points per segment. A taper window and an overlap of 20% are applied for smoother and statistically more significant spectral estimates. The resulting frequency spectra agree well with S t apart from those that correspond to wave frequencies higher than 2 f p , since there is no energy for these frequencies in the MILDwave model. Furthermore, it is observed that the periodic lateral boundaries do not affect the waves since for both wave propagation angles, the frequency spectra are similar.
In Figure 8, a comparison is made between the target frequency spectrum (S ) and the simulated frequency spectra (S and S ) for the cases of θ = 0° and θ = 45°. The surface elevations, η, at the electronic wave gauges, which are positioned at the centre of the computational domain, are recorded from t = 40 T to t = 600 T with a sampling interval of 0.2 s. The recorded data are processed in segments of 2048 points per segment. A taper window and an overlap of 20% are applied fo smoother and statistically more significant spectral estimates. The resulting frequency spectra agree well with S apart from those that correspond to wave frequencies higher than 2f , since there is no energy for these frequencies in the MILDwave model. Furthermore, it is observed that the periodic lateral boundaries do not affect the waves since for both wave propagation angles, the frequency spectra are similar.

Generation of Irregular Short Crested Waves
Finally, we consider a test case of irregular short crested waves. The parameters of the inpu frequency spectrum are identical as for the long-crested irregular waves (H = 1 m, T = 12 s, γ = 3.3) described in the previous section. Waves are generated for a duration of 7200 s with Δt = 0.2 s while to generate short crested waves, two different wave synthesis methods have been employed A method proposed by [14] and a method proposed by [15] as described in detail in Section 2.1.
The directional wave spectrum was measured by a group of five wave gauges placed in the centre of the computational domain, using a "CERC 5" configuration as introduced by [29]. The recorded normalised spreading function distributions, D , are compared with the target distribution

Generation of Irregular Short Crested Waves
Finally, we consider a test case of irregular short crested waves. The parameters of the input frequency spectrum are identical as for the long-crested irregular waves (H s = 1 m, T p = 12 s, γ = 3.3) described in the previous section. Waves are generated for a duration of 7200 s with ∆t = 0.2 s while to generate short crested waves, two different wave synthesis methods have been employed: A method proposed by [14] and a method proposed by [15] as described in detail in Section 2.1.
The directional wave spectrum was measured by a group of five wave gauges placed in the centre of the computational domain, using a "CERC 5" configuration as introduced by [29]. The recorded normalised spreading function distributions, D n , are compared with the target distribution, D n,t , in Figure 9 for both wave synthesis methods and for the spreading standard deviation, σ θ = 10 • (swell waves) and σ θ = 30 • (wind waves). The normalised spreading function is calculated by integrating the calculated 3D spectrum over all wave frequencies and normalizing it. The difference between the two measured normalised spreading function distributions of the two methods is small and the agreement with the target distribution is very good. The Sand and Mynett method agrees slightly better with the target distribution, D n,t , than the Miles method especially for the case of σ θ = 30 • .  Figures 10 and 11 show the plan views of the k values throughout the inner computational domain and the surface elevations at a time instant of t = 7000 s for each wave synthesis method for a mean wave propagation angle of θ = 0° and spreading standard deviation of σ = 10° and σ = 30°. The results are presented in an inner domain with dimensions of 10L × 10L , where L is the wave length that corresponds to the peak wave frequency, f . It is clear that the Sand and Mynett method yields a more homogenous wave field than the Miles method. Specifically, the maximum (E ) value of the calculated absolute errors for the Miles method is 0.051 for σ = 10° and 0.085 for σ = 30°, respectively, while for the Sand and Mynett method, E is 0.013 for σ = 10° and 0.029 for σ = 30°, respectively. It is worth mentioning that in the Sand and Mynett method, the number of wave components, N , is equal to the number of wave frequencies, N. This is happening due to the fact that the wave propagation angles are randomly selected according to the cumulative distribution function of the directional spreading function, D(f, θ), and are assigned to each wave frequency component. In contrast to that, in the Miles method, the number of the wave components, N , is equal to the product of the number of wave frequencies, N, multiplied with the number of wave angles, M. Thus, in the Sand and Mynett method, (i) the computational time, (ii) the number of corrections of the wave propagation angles in order to ensure the periodicity of the waves, and (iii) the non-homogeneity of the generated wave field are much smaller than those for the Miles method, making the Sand and Mynett method preferable when periodic lateral boundaries are applied in a mild-slope wave model. Specifically, the developed model simulation time for a test case of short crested waves of a duration of 7200 s in a numerical basin, which is 2 km long and 2 km wide, is approximately 3 min on a PC with an Intel Core CPU (Central Processing Unit) 2.90 GHz. The results are presented in an inner domain with dimensions of 10 L p × 10 L p , where L p is the wave length that corresponds to the peak wave frequency, f p . It is clear that the Sand and Mynett method yields a more homogenous wave field than the Miles method. Specifically, the maximum (E max ) value of the calculated absolute errors for the Miles method is 0.051 for σ θ = 10 • and 0.085 for σ θ = 30 • , respectively, while for the Sand and Mynett method, E max is 0.013 for σ θ = 10 • and 0.029 for σ θ = 30 • , respectively. It is worth mentioning that in the Sand and Mynett method, the number of wave components, N tot , is equal to the number of wave frequencies, N. This is happening due to the fact that the wave propagation angles are randomly selected according to the cumulative distribution function of the directional spreading function, D(f, θ), and are assigned to each wave frequency component. In contrast to that, in the Miles method, the number of the wave components, N tot , is equal to the product of the number of wave frequencies, N, multiplied with the number of wave angles, M. Thus, in the Sand and Mynett method, (i) the computational time, (ii) the number of corrections of the wave propagation angles in order to ensure the periodicity of the waves, and (iii) the non-homogeneity of the generated wave field are much smaller than those for the Miles method, making the Sand and Mynett method preferable when periodic lateral boundaries are applied in a mild-slope wave model. Specifically, the developed model simulation time for a test case of short crested waves of a duration of 7200 s in a numerical basin, which is 2 km long and 2 km wide, is approximately 3 min on a PC with an Intel Core CPU (Central Processing Unit) 2.90 GHz.
The results show that the MILDwave model with the addition of the periodic boundaries is capable of reproducing a homogeneous short-crested wave field in the whole computational domain as well as the target frequency spectrum and the target directional spectrum. Thus, in future studies, a WEC farm and its effect on the near and far field can be examined under real sea waves by coupling the developed model with a wave-structure interaction solver as described in [19,20]. More precisely, during this coupling methodology, firstly, the short-crested wave field is calculated by the developed model at the position of the WEC farm. Subsequently, this wave field (time series of surface elevation) is used as an input for a wave-structure interaction solver to simulate the diffracted and radiated wave field generated by the presence of the WEC farm. Finally, the resulting wave field is coupled back to MILDwave and propagates throughout the whole numerical domain. As a result, with the coupling methodology in MILDwave, the impact of WEC farms, especially the far field, can be studied using accurate real sea waves. is used as an input for a wave-structure interaction solver to simulate the diffracted and radiated wave field generated by the presence of the WEC farm. Finally, the resulting wave field is coupled back to MILDwave and propagates throughout the whole numerical domain. As a result, with the coupling methodology in MILDwave, the impact of WEC farms, especially the far field, can be studied using accurate real sea waves.

Numerical Validation Using the Vincent and Briggs Shoal Experiment
To demonstrate the strategic importance of periodic boundaries in the generation of a short crested wave field over varying water depths, simulations are conducted for waves propagating over a shoal. A study of regular and irregular wave propagation over an elliptic shoal, which has been performed by [30], is used here for validation of the present numerical model.
The bathymetry of the experimental setup of Vincent and Briggs is implemented in MILDwave, as is illustrated in Figure 12, is defined as: where x and y are the coordinates centered at the center of the shoal and d is the bottom level at any point inside the shoal area. The shoal is similar to that used in the experiments of [31], with a minimum water depth of d = 15.24 cm at the center of the shoal while the area around the shoal has a constant water depth of d = 45.72 cm. In addition, the simulated shoal is located at the center of the numerical domain to take advantage of the lateral periodic boundaries and create a fully

Numerical Validation Using the Vincent and Briggs Shoal Experiment
To demonstrate the strategic importance of periodic boundaries in the generation of a short crested wave field over varying water depths, simulations are conducted for waves propagating over a shoal. A study of regular and irregular wave propagation over an elliptic shoal, which has been performed by [30], is used here for validation of the present numerical model.
The bathymetry of the experimental setup of Vincent and Briggs is implemented in MILDwave, as is illustrated in Figure 12, is defined as:  (24) where x and y are the coordinates centered at the center of the shoal and d e is the bottom level at any point inside the shoal area. The shoal is similar to that used in the experiments of [31], with a minimum water depth of d min = 15.24 cm at the center of the shoal while the area around the shoal has a constant water depth of d = 45.72 cm. In addition, the simulated shoal is located at the center of the numerical domain to take advantage of the lateral periodic boundaries and create a fully developed wave field at the region of the shoal. The grid cell size is ∆x = ∆y = 0.05 m. To obtain a steady state wave field, waves are generated with a duration of 400 s with a time step of ∆t = 0.01 s.  The experiments of Vincent and Briggs are conducted for three kinds of incident breaking and non-breaking waves, i.e., regular, irregular long crested, and irregular short crested waves. In the present study, only irregular long and short crested waves, where large scale breaking is not involved, are examined. The experimental and numerical input wave parameters for non-breaking series are listed in Table 1. Following [30], we use the TMA spectrum [32] as the target wave frequency spectrum in which the spectral energy density, S(f), depends on the parameters, α (Phillip's constant), f (peak wave frequency), γ (peak enhancement factor), and σ (spectral width parameter), as described in Section 2.1. The parameter, γ, is assigned values of 2 for the broad frequency spectrum and 20 for the narrow frequency spectrum. Similarly, these frequency spectra are combined with narrow (σ = 10°) and broad (σ = 30°) directional spreading, while the α value is selected to correspond to the target wave height.
The directional spreading function, D(f, θ), given in Equation (15) is used instead of the wrapped normal spreading function used by [30]. However, the directional spreading parameter, s , is calculated according to Equation (17) in order to achieve the same distribution as the experimental one. Finally, the mean wave propagation angle is θ = 0° and the number of wave components, N, varies from 50 to 400, with the highest value for the case of broad-banded directional spreading. The The experiments of Vincent and Briggs are conducted for three kinds of incident breaking and non-breaking waves, i.e., regular, irregular long crested, and irregular short crested waves. In the present study, only irregular long and short crested waves, where large scale breaking is not involved, are examined. The experimental and numerical input wave parameters for non-breaking series are listed in Table 1. Following [30], we use the TMA spectrum [32] as the target wave frequency spectrum in which the spectral energy density, S(f), depends on the parameters, α (Phillip's constant), f p (peak wave frequency), γ (peak enhancement factor), and σ (spectral width parameter), as described in Section 2.1. The parameter, γ, is assigned values of 2 for the broad frequency spectrum and 20 for the narrow frequency spectrum. Similarly, these frequency spectra are combined with narrow (σ θ = 10 • ) and broad (σ θ = 30 • ) directional spreading, while the α value is selected to correspond to the target wave height.
The directional spreading function, D(f, θ), given in Equation (15) is used instead of the wrapped normal spreading function used by [30]. However, the directional spreading parameter, s 1 , is calculated according to Equation (17) in order to achieve the same distribution as the experimental one. Finally, the mean wave propagation angle is θ = 0 • and the number of wave components, N, varies from 50 to 400, with the highest value for the case of broad-banded directional spreading. The measured frequency spectra are compared with the target spectrum (S t ) in Figure 13 for the two values of the peak enhancement factors used here (see Table 1). The agreement is very good especially for the narrow banded spectrum where the measured one replicates the desired target spectrum. measured frequency spectra are compared with the target spectrum (S ) in Figure 13 for the two values of the peak enhancement factors used here (see Table 1). The agreement is very good especially for the narrow banded spectrum where the measured one replicates the desired target spectrum.  Figure 14 shows the comparison of normalised wave heights, H/H , (where H is the local significant wave height and H is the significant wave height at the wave generation boundary) between numerical model results and experimental data along the measurement transect shown in Figure 12 for the test cases of non-breaking waves listed in Table 1. Additionally, to evaluate the model, the root mean square error (RMSE) and the skill factor for the normalised wave heights are calculated as: where O and P indicate the observed and predicted values, respectively. Very good agreement between the numerical model and the experimental model is observed (Figure 14 and Table 2). MILDwave with the addition of the periodic boundaries correctly predicts the wave focusing behind the shoal for the case of irregular long crested waves where the wave height is strongly affected by the change of the water depth. The model gives a maximum normalised wave height of around 2.05 at y = 0 for both the broad frequency spectrum (Test Case U3) and the narrow frequency spectrum (Test Case U4). On the other hand, a broad directional spreading distribution (Test Case B3, Test Case B4) yields much less spatial wave height variation and the effect of wave refraction is significantly diminished.  Figure 14 shows the comparison of normalised wave heights, H/H 0 , (where H is the local significant wave height and H 0 is the significant wave height at the wave generation boundary) between numerical model results and experimental data along the measurement transect shown in Figure 12 for the test cases of non-breaking waves listed in Table 1. Additionally, to evaluate the model, the root mean square error (RMSE) and the skill factor for the normalised wave heights are calculated as: where O and P indicate the observed and predicted values, respectively. Very good agreement between the numerical model and the experimental model is observed ( Figure 14 and Table 2). MILDwave with the addition of the periodic boundaries correctly predicts the wave focusing behind the shoal for the case of irregular long crested waves where the wave height is strongly affected by the change of the water depth. The model gives a maximum normalised wave height of around 2.05 at y = 0 for both the broad frequency spectrum (Test Case U3) and the narrow frequency spectrum (Test Case U4). On the other hand, a broad directional spreading distribution (Test Case B3, Test Case B4) yields much less spatial wave height variation and the effect of wave refraction is significantly diminished.  Figure 14. Comparison of normalised wave heights, H/H , between numerical model results (blue solid lines) and experimental data (red circles) along the measurement transect of Figure 12. for the test cases of non-breaking waves of Table 1.   Figure 12. for the test cases of non-breaking waves of Table 1.

Conclusions
In the present paper, periodic lateral boundaries were developed in a time dependent mild-slope equation model, MILDwave, for the accurate and fast generation of regular and irregular waves in any direction.
Initially, three different wave generation layouts were examined and compared. The first wave generation layout is an L-shaped wave generator, which is composed of two orthogonal lines. The second wave generation layout is an arc-shaped wave generator where two parallel lines are connected to a semicircle. The third wave generation layout consists of a single wave generation line parallel to the y-axis and periodic lateral boundaries. Numerical experiments were conducted for regular waves with wave propagation angles from θ = 0 • to θ = 45 • . The results were presented by means of contour plots, in terms of a disturbance coefficient (k d ) for the calculation domain, allowing easy comparisons between the different layouts. The wave generation layout with the periodic boundaries showed the best results since this layout leads to a homogeneous wave field, which is not disturbed by unwanted wave diffraction patterns in contrast to the other two wave generation layouts.
Then, this wave generation layout was used to simulate irregular long and short crested waves by using two different wave synthesis methods: A method proposed by Miles [14] and a method proposed by Sand and Mynett [15]. The results indicate that the MILDwave model with the addition of the periodic boundaries is capable of reproducing a homogeneous wave field (especially when the Sand and Mynett method is used) as well as the target frequency spectrum and the target directional spectrum. Finally, the developed model was also used to study wave transformation over an elliptic shoal (Vincent and Briggs shoal experiment), where very good agreement was observed between the numerical model and the experimental results. The aforementioned observations indicate that periodic boundaries make the mild-slope wave model, MILDwave, an essential tool to generate multi-directional waves and study their transformation due to its significantly small computational cost and its high numerical stability and robustness.
The next step is to combine the developed model with the coupling methodology described in [19,20] to study a WEC farm under real sea waves, since this paper proves that now MILDwave is able to accurately generate such a kind of wave field.