*2.1. Ocean Circulation Model*

The sea surface elevations and currents of the waters surrounding Taiwan were simulated with the community ocean model Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM). SCHISM is an upgraded and evolved version of the original Semi-implicit Eulerian-Lagrangian Finite-Element/volume (SELFE [31]) model that includes many enhancements and improvements such as a new extension to simulate large-scale eddying regimes and a seamless cross-scale capability ranging from creeks to oceans [32]. SCHISM and SELFE have been widely used to simulate the propagation of tsunami waves [33], to predict water quality and ecosystem dynamics [34–36], analyze oil spill diffusion and transport [37], generate flooding maps [38,39], evaluate the extent of inundation induced by extreme river flows, typhoons or hurricanes [40,41] and to estimate tidal-current power [42]. In Taiwanese waters, the surface circulation feature is highly similar to that on the

bottom during the spring [43], and weak winter stratification is caused by vertical mixing [44]; therefore, a two-dimensional SCHISM model (SCHISM-2D) vertically integrated with the barotropic mode was employed in this study. The governing equations in the Cartesian coordinate system and two-dimensional form are as follows:

$$\frac{\partial \eta}{\partial t} + \frac{\partial uH}{\partial x} + \frac{\partial vH}{\partial y} = 0,\tag{1}$$

$$\frac{Du}{Dt} = fv - \frac{\partial}{\partial x} \left\{ g(\eta - a\hat{\psi}) + \frac{P\_A}{\rho\_0} \right\} + \frac{\tau\_{\rm sx} + \tau\_{rx} - \tau\_{bx}}{\rho\_0 H},\tag{2}$$

$$\frac{Dv}{Dt} = -fu - \frac{\partial}{\partial y} \left\{ g(\eta - a\hat{\eta}) + \frac{P\_A}{\rho\_0} \right\} + \frac{\tau\_{sy} + \tau\_{ry} - \tau\_{by}}{\rho\_0 H},\tag{3}$$

where *η*(*x, y, t*) is the free-surface elevation, *H* = *η* + *h* is the total water depth, *h* is the bathymetric depth, *u*(*x, y, t*) and *ν*(*x, y, t*) are the horizontal velocity components in the *x*, *y* direction, respectively, *f* is the Coriolis factor, *g* is the acceleration due to gravity, *ψ*ˆ is the Earth's tidal potential, *α* is the effective earth elasticity factor, *ρ*<sup>0</sup> is the reference density of water, and *PA*(*x, y, t*) is the atmospheric pressure at the free surface. In Equations (2) and (3), *τsx* and *τsy* are the wind stress components in the *x, y* direction, respectively, and can be expressed as follows:

$$
\pi\_{sx} = \rho\_a \mathbb{C}\_s \sqrt{\mathcal{W}\_x^2 + \mathcal{W}\_y^2} \mathcal{W}\_x \tag{4a}
$$

$$
\pi\_{sy} = \rho\_a \mathbb{C}\_s \sqrt{\mathcal{W}\_x^2 + \mathcal{W}\_y^2} \mathcal{W}\_y \tag{4b}
$$

where *Cs* is the wind drag coefficient; *ρ<sup>a</sup>* is the air density; and *Wx*, *Wy* are the wind-speed components 10 m above the sea surface in the *x, y* directions, respectively. *Cs* is often regarded as an increasing function of wind velocity; however, Powell et al. [45] suggested that *Cs* should be limited at high wind speeds. The formula to calculate *Cs* in SCHISM-2D is:

$$C\_s = 1.0^{-3} \begin{cases} (0.61 + 0.063 \times 6.0), \mathcal{W} < 6.0\\ (0.61 + 0.063 \mathcal{W}), \ 6.0 \le \mathcal{W} \le 50.0\\ (0.61 + 0.063 \times 50.0), \mathcal{W} > 50.0 \end{cases} \tag{5}$$

where *W* is the resultant wind speed at 10 m above the sea surface.

Here, *τbx* and *τby* are the bottom shear stress components in the *x, y* directions and are computed by the following formulas:

$$
\pi\_{\rm bx} = \rho\_0 \mathbb{C}\_d \sqrt{u^2 + v^2} u \tag{6a}
$$

$$
\pi\_{by} = \rho\_0 \mathbb{C}\_d \sqrt{\mu^2 + \upsilon^2} \upsilon \tag{6b}
$$

where *Cd* is the bottom drag coefficient, which can be parameterized as follows:

$$\mathcal{C}\_d = \mathfrak{g}n^2/H^{1/3} \tag{7}$$

where *n* is the Manning coefficient, which was set to 0.025 in the model due to the lack of knowledge of the sea-bottom material type. However, the bottom drag coefficient *Cd* varies with *H* based on Equation (7).

Here, *τrx* and *τry* are the radiation stress components in the *x, y* directions, respectively, and can be computed following Longuet-Higgins and Stewart [46,47], as shown below:

$$\pi\_{\ell\mathfrak{x}} = -\frac{\partial S\_{xx}}{\partial \mathfrak{x}} - \frac{\partial S\_{xy}}{\partial y} \tag{8a}$$

*Energies* **2018**, *11*, 499

$$\pi\_{ry} = -\frac{\partial S\_{xy}}{\partial x} - \frac{\partial S\_{yy}}{\partial y} \tag{8b}$$

where *Sxx*, *Sxy* and *Syy* are the wave radiation stress components and are represented according to Battjes [48]:

$$S\_{\rm xx} = \int\_0^{2\pi} \int\_0^{\infty} N \sigma \frac{\mathbb{C}\_{\mathcal{S}}}{\mathbb{C}\_p} \sin(\theta) \mathrm{d}\theta \mathrm{d}\sigma,\tag{9a}$$

$$S\_{xy} = S\_{yx} = \int\_0^{2\pi} \int\_0^{\infty} \text{N}\sigma \left[ \frac{\mathbb{C}\_{\mathcal{X}}}{\mathbb{C}\_p} \left( \cos^2(\theta) + 1 \right) - \frac{1}{2} \right] \text{d}\theta \text{d}\sigma \tag{9b}$$

$$S\_{yy} = \int\_0^{2\pi} \int\_0^{\infty} \text{N}\sigma \left[ \frac{\mathbb{C}\_\mathbb{S}}{\mathbb{C}\_p} \left( \sin^2(\theta) + 1 \right) - \frac{1}{2} \right] \text{d}\theta \text{d}\sigma \tag{9c}$$

where *N* is the spectral wave action density, and the independent variables are the wave relative frequency, *σ*, and the wave direction, *θ*. Here, *Cg* and *Cp* are the wave group velocity and the wave phase velocity, respectively.

#### *2.2. Spectral Wind Wave Model*

A third-generation spectral wind wave model, Wind Wave Model III (WWM-III, 3.0, College of William & Mary, Williamsburg, VA, USA), was adopted to predict the significant wave heights (SWH) and wave peak period (*Tp*) in the coastal waters of Taiwan. The wave-action equation that governs WWM-III is given as:

$$\frac{\partial N}{\partial t} + \frac{\partial \left(\mathcal{C}\_{\mathcal{S}^x} + u\right)N}{\partial x} + \frac{\partial \left(\mathcal{C}\_{\mathcal{S}^y} + v\right)N}{\partial y} + \frac{\partial \left(\mathcal{C}\_{\sigma}N\right)}{\partial \sigma} + \frac{\partial \left(\mathcal{C}\_{\theta}N\right)}{\partial \theta} = \frac{\mathcal{S}\_{tot}}{\sigma} \tag{10}$$

where *Cgx* and *Cgy* are the wave group velocity components in the *x, y* direction, *C<sup>σ</sup>* and *C<sup>θ</sup>* are the propagation velocities in *σ*, *θ* space, respectively, and *Stot* is the sum of the source terms for wave variance, including non-linear interactions, wind growth and dissipation by white-capping, bottom friction and wave breaking. The bottom friction and peak enhancement factors are set to 0.067 and 3.3, respectively, based on the formulation of the Joint North Sea WAve Project (JONSWAP) [49]. WWM-III computes waves breaking in shallow-water areas using the method proposed by [50]. A constant wave-breaking coefficient of 0.78 was adopted. The maximum wave direction is 360◦, which is discretized into 36 bins in WWM-III. The low- and high-frequency limits of the discrete wave frequency are 0.03 and 1.0 Hz, respectively and are also divided into 36 bins.

SCHISM-2D (College of William & Mary, Williamsburg, VA, USA) computes the depth-averaged current and water surface elevation and then passes the results to the wave model, WWM-III, which sends the radiation stress to SCHISM-2D after solving the wave action equation. More details about the model coupling procedures of SCHISM-WWM-III can be found in [51].

#### *2.3. Meteorological and Boundary Conditions*

The wind forcing data used to drive SCHISM-WWM-III were the wind field at 10 m above sea level and the sea level pressure produced by ERA-Interim, which were acquired from public datasets available from the European Center for Medium-Range Weather Forecasts (ECMWF). ERA-Interim is the latest global atmospheric reanalysis, which is normally updated monthly with a two-month delay [52]. The 10-m wind speed in the *x*, *y* components were extracted from the ERA-Interim reanalysis results with a temporal resolution of 6 h (four analysis fields per day, at 00:00, 06:00, 12:00 and 18:00 UTC) and a spatial resolution of 0.125◦ × 0.125◦ and converted to the unstructured grids of SCHISM-WWM-III using the Inverse Distance Weighting (IDW) method. The IDW formula is:

$$p\_{\ell}(x, y) = \frac{\sum\_{i=1}^{n} \frac{1}{d\_i} p\_m(x\_{i\prime} y\_i)}{\sum\_{i=1}^{n} \frac{1}{d\_i}}\tag{11}$$

where *pe*(*x, y*) is the value to be estimated (the wind speed in SCHISM-WWM-III), *pm*(*x, y*) is the known value (the wind speed in ERA-Interim), and *di* represents the distances from the *n* data points to the point estimated *n.*

The wave boundary conditions of the regional SCHISM-WWM-III, including the significant wave height, the mean wave direction, the mean directional spreading, the peak frequency, and the zero down-crossing frequency, were obtained from a global WaveWatch III (WW III) model developed by the French Research Institute for Exploitation of the Sea (IFREMER) with a spatial resolution of 0.5◦ × 0.5◦ and a temporal resolution of 3 h. To obtain the tide boundary conditions, the harmonic constants of eight tidal constituents (M2, S2, N2, K2, K1, O1, P1, and Q1) were extracted from a regional inverse tidal model, China Seas & Indonesia 2016 [53], with a resolution of 1/30◦ and were set at each node on the open boundary to drive the model.

#### *2.4. Estimation of Wave Power*

The extractable wave power per unit of the wave crest length can be calculated through the spectral output of SCHISM-WWM-III and is given in kW per meter as shown below:

$$P = \rho \mathop{\rm g} \int\_0^{2\pi} \int\_{\mathcal{S}} c\_{\mathcal{S}}(\sigma, d) \mathcal{S}(\sigma, \theta) \mathrm{d}\sigma \mathrm{d}\theta \,\tag{12}$$

where *S(σ*, *θ*) is the directional wave variance density spectrum, *ρ* is the density of seawater, *g* is the acceleration due to gravity, and *d* is the water depth. The wave group velocity, *cg(σ*, *θ*), can be expressed as follows:

$$c\_{\mathcal{S}}(\sigma, d) = \frac{\mathcal{g}}{4\pi\sigma} \left[ 1 + \frac{2\kappa d}{\sin h(2\kappa d)} \right] \tan[h(\kappa d)] \tag{13}$$

where *κ =* 2π/*L* is the number of waves, and *L* is the wave length. In deep water conditions (*d* > 0.5*L*) and *cg(σ*, *θ*) = *g*/4π*σ*, Then:

$$P = \frac{\rho\_{\mathcal{S}}^{\ast}}{4\pi} \int\_{0}^{2\pi} \int\_{0}^{\infty} \sigma^{-1} S(\sigma, \theta) \mathbf{d}\sigma \mathbf{d}\theta \tag{14}$$

The spectral moments of order, *n*, are defined as follows:

$$m\_{\rm nl} = \frac{\rho g^2}{4\pi} \int\_0^2 \int\_0^\infty \sigma^{\rm n} S(\sigma, \theta) \mathrm{d}\sigma \mathrm{d}\theta \tag{15}$$

The energy period *Te* and the significant wave height *Hs* in terms of spectral moment are as follows:

$$T\_t = \frac{m\_{-1}}{m\_0} = \frac{\int\_0^{2\pi} \int \sigma^{-1} S(\sigma, \theta) \mathrm{d}\sigma \mathrm{d}\theta}{\int\_0^{2\pi} \int S(\sigma, \theta) \mathrm{d}\sigma \mathrm{d}\theta} \tag{16}$$

$$H\_s = 4(m\_0)^{1/2} = 4 \left( \int\_0^{2\pi} \int \mathbf{S}(\sigma, \theta) \mathbf{d}\sigma d\theta \right)^{1/2} \tag{17}$$

Therefore, Equation (11) can be simplified to:

$$P = \frac{\rho \mathcal{g}^2}{64\pi} H\_s^2 T\_\varepsilon \tag{18}$$

If *ρ* is assumed to be 1025 kg/m3, Equation (18) becomes:

$$P(\text{kilowatt}/\text{meter}, \text{ kw}/\text{m}) = 0.49 H\_\text{s}^2 T\_\text{g} \tag{19}$$

It has been suggested that *Te* is rarely specified and must be estimated from other variables, while the peak period *Tp* is often specified by measurements of sea states [23,49,50]. *Te* is equal to *Tp* when multiplied by a coefficient *α*:

$$T\_{\mathfrak{E}} = \mathfrak{a}T\_{\mathfrak{p}} \tag{20}$$

in which *α* = 0.9 was adopted from the standard JONSWAP spectrum with a peak enhancement factor of *γ* = 3.3 [23,27,54,55]. Thus, the theoretical wave power at any point can be computed using the following expression:

$$P(\text{kilowatt}/\text{meter}, \text{ kw}/\text{m}) = 0.44 H\_\text{s}^2 T\_p. \tag{21}$$

The present study used Equation (15) to evaluate the wave power resources in Taiwanese waters, even when the deep-water hypothesis was not strictly satisfied [56].

Thus, the total annual wave energy per unit length can be generated for a given time interval and is calculated as follows:

$$E\_l(\text{kilowatt hour} / \text{meter}, \text{ kwh} / \text{m}) = \sum\_{i} P\_{i} \times \Delta t\_{i\prime} \tag{22}$$

where Δ*ti* is the temporal sampling interval; Δ*ti* = 1 h in the present study.

#### *2.5. Model Configuration*

The computational domain covers Taiwan in its entirety as well as its main offshore islands. The region covers the area from a longitude of 114◦ E to 130◦ E and a latitude of 19◦ N to 29◦ N, as shown in Figure 1a. Two gridded bathymetry datasets were employed in this study. One dataset consists of global data obtained from the General Bathymetric Chart of the Oceans (GEBCO) and has a resolution of 30 arc-seconds. It was generated by combining quality-controlled ship depth soundings with interpolations between the sounding points guided by satellite-derived gravitational data. The other dataset is local data provided by the Department of Land Administration, Ministry of the Interior in Taiwan, with a resolution of 200 m and is distributed over a longitude of only 100◦ E to 128◦ E and a latitude of 4◦ N to 29◦ N. The bathymetry data for the coupled model were produced by merging the GEBCO and local bathymetry datasets. The method for combining two bathymetry sets is to merge the values in some order of priority using Surface-water Modeling System (SMS) software. A priority must be adopted because data overlaps exist in the source datasets. In the present study, when overlap occurs, the GEBCO data points in a set are considered as lower priority than the local data points; thus, they are not included in the final merged set. When no overlap occurs, all data points are merged. This approach ensured the preservation of the high-resolution local data. Our modeling domain is composed of 278,630 triangular cells and 142,041 non-overlapping, unstructured grids. Coarse meshes with 30 km resolution were arranged on the open ocean boundaries, while fine meshes with 300 m resolution were used along the coastline of Taiwan and the small offshore islands (Figure 1b). After the unstructured grids were created, the final merged, gridded bathymetry dataset

was interpolated to each node to represent the bottom elevations in the coupled model (Figure 1a) using the linear interpolation method in SMS software.

In the ocean circulation model, a time step of 120 s was used for the present unstructured-grid system with numerical stability. SCHISM and WWM-III exchange the computational results every five hydrodynamic time steps (i.e., the time step of WWM-III is 600 s).

**Figure 1.** (**a**) Bathymetry and (**b**) unstructured grid of the computational domain.

## **3. Model Validation**

The observations adopted to validate SCHISM-WWM-III were measured significant wave heights (SWH) and peak periods (*Tp*) provided by the Harbor and Marine Technology Center (HMTC) in Taiwan at four buoys. The distribution map of these four wave buoys around Taiwanese waters is shown in Figure 2. The Keelung, Suao, and Hualien buoys are situated along the northeastern and eastern coast of Taiwan, while the Penghu buoy lies on the western coast of Taiwan (the Taiwan Strait). Figures 3 and 4 present the comparisons of SWH (Figure 3) and *Tp* (Figure 4) between the model hindcasts and the measurements at the four buoys from 31 January 2016 to 30 December 2016. This period includes normal wind-generated waves (*Tp*), monsoon-induced waves (*Tp*) and the abnormal waves (*Tp*) caused by strong winds from Typhoon Nepartax (a severe typhoon, from 6 July 2016 to 9 July 2016), Typhoon Meranti (a severe typhoon, from 12 September 2016 to 15 September 2016) and Typhoon Megi (a moderate typhoon, from 25 September 2016 to 28 September 2016). The maximum SWH and *Tp* reached approximately 5.5 m and 14 s, 9.5 m and 16 s, 7.5 m and 15 s and 3.5 m and 13 s at the Keelung, Suao, Hualien, and Penghu buoys, respectively, due to the approach of Typhoon Megi. Peak SWH occurred at Keelung, Suao, and Hualien buoys on 28 September 2016, with a lag of almost one day at the Penghu buoy (as shown in Figures 3 and 4). This result occurred because Typhoon Megi crossed Taiwan from east to west.

The statistical errors of the differences between the model hindcasts and the observations at the four wave buoys are also estimated. The correlation coefficients of the SWH and *Tp* are 0.86 and 0.74, 0.83 and 0.84, 0.87 and 0.83 and 0.80 and 0.73 for the Keelung, Suao, Hualien, and Penghu buoys, respectively. Even though the hindcasted SWH and *Tp* agree well with the measurements during both ordinary and extraordinary meteorological conditions, slight discrepancies still exist. It must be emphasized that the numerical wave models fail to capture *Tp* well, especially during normal meteorological conditions, for instance, April, May and June. The possible reasons include the nature of wave models. Some of the setup assumptions and numerical solutions within the models affect their accuracy [57]. In addition, the lower spatial and temporal resolutions of the wind field data from the ERA-Interim dataset affected model hindcasting, which is more sensitive to *Tp*. The adoption of higher spatial and temporal resolution wind data could improve the performance [58]. Yet a third reason for poor *Tp* hindcasts may be due to observational error. However, based on the model-data comparison, the hindcasts of SWH and *Tp* by SCHISM-WWM-III are relatively reliable and can be further applied to assess the distribution of wave power and wave energy in the waters surrounding Taiwan.

**Figure 2.** Distribution of wave buoys around Taiwanese waters. The cyan area represents ocean and the white areas represent land.

**Figure 3.** *Cont.*

**Figure 3.** Mode-data comparison of significant wave height (SWH) at four wave buoys. (**a**) Keelung; (**b**) Suao; (**c**) Hualien and (**d**) Penghu.

**Figure 4.** *Cont.*

**Figure 4.** Mode-data comparison of peak period (*Tp*) at four wave buoys. (**a**) Keelung; (**b**) Suao; (**c**) Hualien and (**d**) Penghu.

#### **4. Results and Discussion**

Sea-state hindcasts covering the period from 2005 to 2016 were conducted to estimate the wave power and wave energy output. Additionally, numerical simulations based on unstructured grids were converted to structured grids with a resolution of 25 × 25 km to locate the energetic sea areas.

#### *4.1. Spatial Distribution of Annual Mean Significant Wave Height and Wave Power*

Figure 5 illustrates the spatial distribution of offshore annual mean SWH in Taiwanese waters for each year from 2005 to 2016 and shows highly similar patterns of mean SWH distribution. The mean SWH gradually increased from the nearshore region (shallower water) to the offshore areas (deeper water). This predictable phenomenon is due to the dissipation of SWH by depth effects [59]. The offshore sea areas with the higher SWH lie off the northeastern, eastern and southern coasts of Taiwan. The SWH values range from 1.2 m to 1.6 m along the northeastern coast and 1.4 m–1.8 m along the eastern and southern coasts. Figure 5 also shows that the mean SWH values are relatively low in the Taiwan Strait, ranging from only 0.6 m to 1.2 m. The spatial distribution of offshore SWH was averaged over 12 years (2005–2016) and is presented in Figure 6. The offshore sea areas southeast of Taiwan exhibit the highest mean SWH values, ranging from 1.2 m to 1.8 m. Although the prevalence

of the strong northeast monsoon plays an important role in the higher SWH values, the deeper waters in this area also contribute to the SWH increases.

The SWH and *Tp* were output hourly from SCHISM-WWM-III to compute the wave power for each grid using Equation (21). Figure 7 demonstrates the spatial distribution of offshore annual mean wave power (in kW/m) in the waters surrounding Taiwan from 2005 to 2016. The more energetic sea areas are along the northeastern, eastern and southern coasts of Taiwan, and these locations are consistent with the distributions of the higher SWH. The mean wave power values were 6–12 kW/m for the northeastern coast and 8–16 kW/m for the eastern and the southern coast. The offshore wave power distribution for the 12-year annual mean is displayed in Figure 8. The most energetic areas are the eastern and the southeastern waters of Lanyu, with mean wave power values of 12–14 kW/m.

**Figure 5.** Spatial distribution map of offshore annual average significant wave height (SWH) for each year from 2005 to 2016 in Taiwanese waters.

**Figure 6.** Spatial distribution map of offshore 12-year (2005 to 2016) annual average significant wave height (SWH) in Taiwanese waters.

**Figure 7.** Spatial distribution map of offshore annual average wave power for each year from 2005 to 2016 in Taiwanese waters.

**Figure 8.** Spatial distribution map of offshore 12-year (2005–2016) annual average wave power in Taiwanese waters.

#### *4.2. Identifying the Optimal Offshore Areas for Wave Energy Converter Deployment*

As shown in Figure 8, the offshore areas with the highest wave power density were observed off the northern, northeastern, southeastern (south of Green Island and southeast of Lanyu) and southern coasts of Taiwan. The annual total wave energy yields based on the 12-year average for each 25 × 25 km square grid were calculated by applying Equation (22) and are shown in Figure 9 (left panel). Even though higher wave energies exist in several sea areas, especially in the southeast (105–120 kW/m), the distance between coast and wave energy converter (WEC) should also be considered because of the cabling costs [19]. The five 25 × 25 km energetic sea areas marked with deep blue lines, with spatial average annual total wave energy density of 60–90 kW/m were selected for further analysis (as shown in the left panel of Figures 9–13). The selected 25 × 25 km grid for the northern coast of Taiwan was downscaled to a resolution of 5 × 5 km. A total of 25 square grids are exhibited in the right panel of Figure 9. The grids are within 5 km of the coasts, have higher wave energies and were determined to be optimal areas for WEC deployments (see S1 in the right panel of Figure 9). The same approach was adopted to select an optimal area for the northeastern coast of Taiwan (S2 in the right panel of Figure 10), for the southern coast of Green Island (S3 in the right panel of Figure 11), for the southeastern coast of Lanyu (S4 in the right panel of Figure 12) and for the southern coast of Taiwan (S5 in the right panel of Figure 13). The spatial average (average over a square of 5 × 5 km) annual total wave energy densities are estimated to be 64.3, 84.1, 84.5, 111.0 and 99.3 MWh/m at S1, S2, S3, S4 and S5, respectively. The central longitudes, latitudes, and average water depths for these five optimal areas are listed in Table 1. In other words, the sea areas within 2.5 km of the central coordinate (as listed in Table 1) are considered to be appropriate sites for WEC deployments.

Tsai et al. [21] utilized the SWAN (Simulating WAves Nearshore) driven by NECP (National Centers for Environmental Prediction) global reanalysis wind field to hindcast waves and periods around Taiwanese waters. They revealed that higher significant wave heights were found near Lanyu. Chiu et al. [22] analyzed the temporal and spatial characteristics of wave power density in the coastal areas of Taiwan wave parameter hindcasts (significant wave height and period) based on the TaiCOMS model and observational data. Their results also indicated that the maximum annual mean wave power was found in Lanyu. Su et al. [23] employed the WWM-III model to assess the distribution of wave power in the waters surrounding Taiwan and suggested that the most energetic sea area is southeast of Lanyu. Thus, our results are consistent with those of previous studies.

Wave direction is another important parameter that should be considered in WEC placements [60]. This is particularly true for attenuator-type and terminator-type Weeks: attenuators should be parallel to the predominant wave direction while the principal axis of terminator devices should be perpendicular to the predominant wave direction [61]. Point absorbers are buoy-type WECs that harvest incoming wave energy from all directions, however, arrays of this type of WEC can be sensitive to wave direction when they are deployed in offshore regions [6]. Figure 14 shows wave rose plots for the 12-year average at the centers of S1, S2, S3, S4 and S5. The dominant wave directions for these five optimal areas lie between east and northeast.

**Figure 9.** Spatial distribution map of offshore 12-year annual total wave energy output with (**a**) a resolution of 25 × 25 km and (**b**) a downscaled 5 × 5 km for optimal area 1 (S1).

**Figure 10.** Spatial distribution map of offshore 12-year annual total wave energy output with (**a**) a resolution of 25 × 25 km and (**b**) a downscaled 5 × 5 km for optimal area 2 (S2).

**Figure 11.** Spatial distribution map of offshore 12-year annual total wave energy output with (**a**) a resolution of 25 × 25 km and (**b**) a downscaled 5 × 5 km for optimal area 3 (S3).

**Table 1.** Central latitudes, longitudes, and average water depths of five optimal areas for wave energy converter (WEC) deployment.


\* Depth is below the sea level.

**Figure 12.** Spatial distribution map of offshore 12-year annual total wave energy output with (**a**) a resolution of 25 × 25 km and (**b**) a downscaled 5 × 5 km for optimal area 4 (S4).

**Figure 13.** Spatial distribution map of offshore 12-year annual total wave energy output with (**a**) a resolution of 25 × 25 km and (**b**) a downscaled 5 × 5 km for optimal area 5 (S5).

**Figure 14.** *Cont.*

**Figure 14.** Directional distributions of average wave power for five optimal areas (**a**) S1; (**b**) S2; (**c**) S3; (**d**) S4; (**e**) S5 for 2005–2016.

#### *4.3. Discussion*

In many coastal areas, nearshore regions have been considered to be the ideal locations for deploying wave energy converters due to their lower costs. However, by comparing wave lengths to wave energy conversion, the nearshore areas can be regarded as shallow-water areas. The effect of water depth on wave propagation is significant and cannot be neglected for wave energy assessment in shallow-water regions [62]. However, the effects of water depths are typically ignored because studies usually employ the formulas for deep-water waves. This simplification may underestimate the wave energy in finite water depths. Sheng and Li [63] conducted wave power assessment at a water depth of 50 m using the deep-water formula and the proposed method in their study and compared the results to the actual values. They found that the deep-water formula underestimated the annual mean wave power by up to 10.18% at the studied finite water depth (50 m). However, when the water depth modification factor (*Ch*) was considered, the underestimation of the annual mean wave power fell to 1.47%. Their results suggested that a modification of the deep-water formula is necessary to accurately assess wave power, especially in cases that involve finite water depth effects. Additionally, tide levels and wind intensity may influence the wave characteristics in shallow-water locations [64]. These effects are considered in our tide-surge-wave fully coupled model.

Additional uncertainties in wave power density assessments might be introduced by using a deep-water assumption (Equations (14)–(22)). However, this assumption causes obvious errors only in shallow water areas (i.e., those with water depths <0.5 × wavelength) [65]. The more energetic sea areas such as S1–S5 are located in the northern, eastern and southern offshore waters of Taiwan and can be regarded as deep-water areas. Therefore, the uncertainty of wave power density using Equation (21) can be ignored.

Guillou [66] employed SWAN (Simulating Waves Nearshore model) and TOAWAC (TELEMAC-based Operational Model Addressing Wave Action Computation) to evaluate wave power for the Iroise Sea. The results revealed that SWAN calculates a 15% lower mean wave power than does TOAWAC in offshore waters. Therefore, inter-model comparisons of wave power computation will be performed for Taiwanese waters in the near future.

The wind field data with a temporal resolution of 6 h and a spatial resolution of 0.125◦ by 0.125◦ from the ERA-interim dataset seems too low to accurately force a tide-surge-wave coupled model, although it is convenient for providing meteorological boundary conditions and producing an acceptable sea-state hindcast. Consequently, a three-layer nested WRF (Weather Research and Forecast, [67]) model with a temporal resolution of 1 h and a spatial resolution of 45 × 45 km for

domain 1, 15 × 15 km for domain 2 and 5 × 5 km for domain 3 will be considered in the future as an alternative source of meteorological information to produce more accurate wave power estimations.

#### **5. Conclusions**

A tide-surge-wave fully coupled model based on an unstructured grid system, SCHISM-WWM-III, was implemented to simulate the sea states in the waters surrounding Taiwan. The wind field 10 m above sea level from the ERA-Interim reanalysis data and harmonic constants of eight tidal constituents extracted from a regional inverse tidal model were used to drive SCHISM-WWM-III. The fully coupled model was verified against available significant wave height and peak period measurements. Twelve-year model hindcasts for 2005–2016 were conducted to evaluate the spatial distributions of wave power and wave energy and to determine the most energetic areas in the offshore seas of Taiwan. Numerical simulations with unstructured grids were converted to 25 × 25 km structured grids to identify the most energetic sea areas. An offshore wave power and wave energy distribution map was created for the 12-year annual mean. Five 25 × 25 km energetic sea areas with spatial average annual total wave energy densities of 60–90 kW/m were selected and subsequently downscaled to resolutions of 5 × 5 km. The 5 × 5 km grids are within 5 km of the coasts and have higher wave energies and were determined to be optimal areas for WEC deployments. Five 5 × 5 km square areas off the northern, northeastern and southern coasts of Taiwan (S1, S2 and S5), the southern coast of Green Island (S3) and the southeastern coast of Lanyu (S4) were selected as the optimal areas for WEC deployments. The spatial average annual total wave energy densities are estimated to be 64.3, 84.1, 84.5, 111.0 and 99.3 MWh/m for the areas S1, S2, S3, S4 and S5, respectively. The dominant wave directions for these five optimal areas lie between east and northeast.

**Acknowledgments:** This research was supported by the Ministry of Science and Technology (MOST), Taiwan, grant Nos. MOST 106-2625-M-865-001. The authors would like to thank the Harbor and Marine Technology Center, Taiwan, for providing the observational data and Joseph Zhang at the Virginia Institute of Marine Science, College of William & Mary, and Aron Roland at the Institute for Hydraulic and Water Resources Engineering, Technische Universität Darmstadt, for kindly sharing their experiences concerning the use of the numerical model.

**Author Contributions:** Hung-Ju Shih, Chih-Hsin Chang, Lee-Yaw Lin and Wei-Bo Chen conceived the study, and Wei-Bo Chen performed the model simulations. The final manuscript has been read and approved by all the authors.

**Conflicts of Interest:** The authors declare no conflict of interest.
