Effect of Depth-Induced Breaking on Wind Wave Simulations in Shallow Nearshore Waters off Northern Taiwan during the Passage of Two Super Typhoons

: Super Typhoons Maria (2018) and Lekima (2019) were adopted for this case study, although they only passed the northern offshore waters of Taiwan without making landfall. A direct modiﬁcation technique was employed to create the atmospheric conditions for a wave-circulation model to hindcast large typhoon-driven waves. The radius of the modiﬁed scale ( R trs ) for a hybrid typhoon wind plays an important role in the signiﬁcant wave height (SWH) simulations during the passage of typhoons. The maximum increment in peak SWH reached 3.0 m and 5.0 m in the deep ocean for Super Typhoons Maria (2018) and Lekima (2019), respectively if the R trs was increased from 4 × R max (radius of the maximum wind) to 7 × R max . The SWHs induced by the typhoon winds in the surf zone were more sensitive to different wave-breaking formulations used in the wave-circulation model. The maximum difference in peak SWH reached 2.5 m and 1.2 m for Super Typhoons Maria (2018) and Lekima (2019), respectively, when the wave-breaking formulations of BJ78 (proposed by Battjes and Janssen in 1978) and CT93 (proposed by Church and Thornton in 1993) were introduced to the wave-circulation model. The SWH simulations in the surf zone were insensitive to the wave-breaking criterion ( γ ) during the passage of typhoons. In shallow nearshore waters, the utilization of a constant γ for the wave-circulation model always produces peak SWHs that are smaller than those using γ based on local steepness or peak steepness. waves induced by Super Typhoon Maria in 2018 and Super Typhoon Lekima in 2019 via model validations. A series of numerical experiments were conducted using the SCHISM-WWM-III modeling system incorporated with modiﬁed typhoon winds to better understand the wave hydrodynamics in the shallow nearshore waters off northern Taiwan. The results derived from the designed numerical experiments reveal that the wave breaking formulations inﬂuence the hindcast of storm waves in the surf zone of northern Taiwan. The maximum difference in peak SWH could reach 2.5 m and 1.2 m for Super Typhoons Maria (2018) and Lekima (2019), respectively, when the wave-breaking formulations of BJ78 and CT93 were introduced to the SCHISM-WWM-III modeling system. Regarding the wave-breaking criterion on the hindcast of typhoon waves in the surf zone of northern Taiwan, compared with the wave-breaking formulation, the maximum difference in peak SWHs was relatively non-sensitive to the wave-breaking criterion. The maximum difference in peak SWH is only 0.5 m using the constant breaking criterion, the breaking criterion based on local steepness or the breaking criterion based on peak steepness. Another important ﬁnding is that the utilization of the BJ78 wave-breaking formulation usually underpredicts the typhoon-generated SWHs in shallow nearshore waters than other parametrizations (i.e., TG83 and CT93). In future research, it will be important to acquire measurements in the surf zone where wave hydrodynamics are more sensitive to wave-breaking formulations and criteria during the passage or landfall of typhoons.


Introduction
Ocean surface wind waves are a dominant process in coastal, nearshore and offshore regions worldwide. Understanding the characteristics of cyclone-driven extreme waves, their variability and historical and projected future changes are important considerations for the sustainable development of coastal and offshore infrastructures and the management of coastal resources and ecosystems. The energy in ocean surface waves is transmitted from the wind. The wind patterns above the ocean surface have been affected as the upper ocean has warmed, consequently resulting in stronger ocean waves. According to the report from Reguero et al. [1], the energy of ocean waves, i.e., wave height, has grown over the past seven decades, which could have significant implications for coastal communities and ecosystems. Since the contributions of ocean surface waves to extreme total water levels are substantial at open coasts, they have mostly been considered in many local studies [2][3][4][5].
Predicting wave heights accurately in coastal and nearshore areas is essential for a number of human activities there, such as renewable energy applications, aquaculture, maritime transport and infrastructure; furthermore, information on both swells and wind seas is important to coastal applications. For example, the prediction of locally generated 2 of 21 wind seas is essential for high-speed passenger ferries, whereas information about the propagation of low-frequency swells in coastal areas and fjords is critical for the design of coastal structures or in the planning of marine operations.
Depth-induced wave breaking is one of the most dominant hydrodynamic processes occurring in coastal regions. This wave breaking not only controls the amount of wave energy impacting our coastlines and coastal defenses but also plays a crucial role in driving many nearshore processes, such as sediment transport, bottom morphology [6] and turbulence, which has been shown to be very important for local ecology [7]. Wave breaking also induces radiation stresses that drive wave-induced setup and currents [8], both of which are important for coastal engineering design and management. However, despite the importance and relevance toward our knowledge of wave hydrodynamics, depth-induced wave breaking is still poorly understood, which is partially due to its highly nonlinear nature; therefore, it is heavily parameterized in most wind wave models.
Many parameterizations for the wave-breaking index have been proposed and reported in previous studies. These parameterizations include dependencies of the wavebreaking index on the offshore wave steepness [9,10], a dissipation rate based on a normalized surf zone width [11], the offshore wave height and the inverse Iribarren number [12]. Ruessink et al. [13] also introduced a parameterization of the wave-breaking index that linearly increases with the local nondimensional depth based on the peak period.
This study provides a critical and objective assessment of three specialized depthinduced wave-breaking models and wave-breaking criteria, which are based on state-ofthe-art breaking wave formulations. The aim of this paper is to study the effect of the depth-induced wave-breaking formulation and wave-breaking criteria on the hindcasts of SWH in shallow nearshore waters off northern Taiwan during the passage of Super Typhoons Maria in 2018 and Lekima in 2019, as well as to better understand which hindcast is more influential in wave height hindcasting inside the surf zone. The details of the materials and methods are presented in the following section, and the results of model validation and a series of designed model experiments are described in Section 3. A discussion and uncertainties of the present study are given in Section 4. Finally, Section 5 summarizes and concludes the paper.

Description of the Study Cases
Super Typhoon Maria was a powerful tropical cyclone that affected Guam (the United States), the Ryukyu Islands (Japan), Taiwan and East China in early July 2018. Maria became a tropical storm and passed the Mariana Islands on July 4 and rapidly intensified the next day due to favorable environmental conditions. Maria reached its first peak intensity on 6 July, and a second stronger peak intensity with 1-min sustained winds of 270 km/h (equivalent to category 5 super typhoon status on the Saffir-Simpson scale) and a minimum pressure of 915 hPa was reached on 9 July. Maria finally made landfall over the Fujian Province, China, early on 11 July after crossing the Yaeyama Islands and passing the northern offshore waters of Taiwan on 10 July. Figure 1a demonstrates the track and arrival times of Super Typhoon Maria in 2018.
Super Typhoon Lekima originated from a tropical depression that was developed in the eastern Philippines on 30 July 2019. Lekima became a tropical storm and was named on 4 August. Under favorable environmental conditions, Lekima intensified and reached its peak with 1-min sustained winds of 250 km/h (equivalent to a category 4 super typhoon status on the Saffir-Simpson scale) and a minimum pressure of 925 hPa on 8 August. Lekima made landfall in the Zhejiang Province, China, on late 9 August and made its second landfall in the Shandong Province, China, on 11 August after moving across eastern China. The track and arrival times of Super Typhoon Lekima are illustrated in Figure 1b. Super Typhoon Lekima originated from a tropical depression that was developed in the eastern Philippines on 30 July 2019. Lekima became a tropical storm and was named on 4 August. Under favorable environmental conditions, Lekima intensified and reached its peak with 1-min sustained winds of 250 km/h (equivalent to a category 4 super typhoon status on the Saffir-Simpson scale) and a minimum pressure of 925 hPa on 8 August. Lekima made landfall in the Zhejiang Province, China, on late 9 August and made its

Information on Offshore Wave Buoys
Three wave buoys, namely, Fuguijiao, Longdong and Suao, located in the northern and northeastern offshore waters of Taiwan, were selected for model validation because they are closest to the tracks of Super Typhoons Maria in 2018 and Lekima in 2019 (as shown in Figure 1). The sampling frequency of wave buoys is 2 Hz for 10 min at the beginning of each hour with an accuracy of ±10 cm for the SWH measurements according to the annual buoy observation data report from the CWB. Information on the coordinates of the three wave buoys and their corresponding water depths is listed in Table 1.

Direct Modification of Typhoon Winds from ERA5
Since the mid-1960s, many parametric cyclone wind models have been proposed [14,15] and widely used to mimic the wind distribution of typhoons [16][17][18][19] because of their simplicity. Parametric cyclone wind models could be used to accurately reconstruct the wind distributions near the center of the typhoon; however, they are unable to accurately reproduce wind speeds in regions far from the center of the typhoon. In contrast, reanalysis wind data obtained from the dynamical model with data assimilation show a superior performance for hindcasting the winds outside of the typhoon's center but are generally inferior for the hindcasts of maximum typhoon wind speed [20][21][22][23][24]. A direct modification technique recommended by Pan et al. [20] was applied in the present study to take advantage of combining the parametric cyclone wind model and reanalysis wind data and maintain a reliable structure for the entire typhoon wind field.
where W DM is the wind speed at an arbitrary grid within the model domain through direct modification, W ERA5 is the wind speed extracted from ERA5 (the fifth-generation reanalysis of the European Centre for Medium-Range Weather Forecasts for the global climate and weather) at an arbitrary point in the computational grid, W Bmax is the maximum wind speed of the best track typhoon issued by the Regional Specialized Meteorological Center (RSMC) Tokyo-Typhoon Center, W Emax is the maximum wind speed of the typhoon among the hourly ERA5 wind fields, r is the radial distance from an arbitrary grid within the model domain to the eye of the typhoon, R trs is the radius of the modified scale (also known as the radius of the transitional zone) and R max is the radius at the maximum typhoon wind speed. R max can be expressed as a function of W Bmax and the latitude of the typhoon's center: where φ is the latitude of the typhoon's center. In Equation (2), m 0 , m 1 and m 2 were set to 38.0 (in n·mi), −0.1167 (in n·mi·kt −1 ) and −0.0040 (in n·mi o−1 ), respectively, according to the results derived from Knaff et al. [25] for the Western Pacific typhoon basin. R trs is considered to be a key factor in determining the accuracy of wind fields; therefore, various R trs will be employed to create hybrid typhoon wind fields to better understand their effect on wind wave hindcasting.

Configuration of the Wave-Circulation Modeling System
A seamless cross-scale hydrodynamic model based on an unstructured grid and triangular mesh served as the ocean circulation model in the wave-circulation modeling system. The hydrodynamic model is called SCHISM (semi-implicit cross-scale hydroscience integrated system model), which has been implemented by Zhang et al. [26] and other developers around the world. Similar to the SELFE (semi-implicit Eulerian-Lagrangian finite element/volume model, the predecessor of the SCHISM developed by Zhang and Bap-tista [27]), the SCHISM also avoids the severest stability constraints in the numerical model by means of a highly efficient semi-implicit scheme [28]. Hence, the high-performance calculations can be performed even when a very high spatial resolution mesh is used in the SCHISM. The splitting between internal and external modes could derive a numerical error [29], which is eliminated through the no-mode-splitting technique implemented in the SCHISM. A depth-averaged (two-dimensional (2D)) ocean circulation model is sufficient for simulating typhoon-driven hydrodynamics. Additionally, fewer computing resources and execution times are required for a 2D model. Therefore, a 2D model, i.e., SCHISM-2D, is selected for wind wave modeling in the present study. The SCHISM and its predecessor SELFE are multipurpose models that have been widely applied to the simulation of hydrodynamics and water quality transportation in coastal and estuarine environments in Taiwan, e.g., evaluation of storm tide-induced coastal inundation [30,31], assessment of tidal stream energy [32,33], transport of suspended sediment and fecal coliform [34,35]. The Manning coefficient and time step were set as 0.025 and 120 s for the barotropic ocean model, respectively, according to the geological characteristics of the seafloor in the waters surrounding Taiwan and the numerical stability of the SCHISM-2D.
The WWM-III is a derivative work from the original WWM-II (wind wave model version III, developed by Roland [36]). The WWM-III is a third-generation spectral wave model that is able to simulate and predict the ocean surface sea state [37]. The wave action balance equation governing the WWM-III is solved by the fractional step method on an unstructured grid. The number of directional bins is 36 with a minimum of 0 • and a maximum direction of 360 • . The number of frequency bins is 36, the lowest limit of the discrete wave frequency is 0.04 Hz and a highest frequency limit of the discrete wave period is 1.0 Hz. The peak enhancement factor is specified as 3.3 for the JONSWAP (Joint North Sea Wave Project, [38]) spectra, while the wave breaking coefficients for the constant and bottom friction coefficients are 0.78 and 0.067, respectively, in WWM-III.
To enhance the information exchange efficiency between the ocean circulation and wind wave models and to eliminate the interpolation errors from the two models, SCHISM-2D and WWM-III take advantage of sharing the same subdomains and parallelization through the same domain decomposition scheme. Additionally, time steps of 120 s and 600 s were used in the SCHISM-2D and WWM-III models, respectively, to improve the computational performance of the coupled model, SCHISM-WWM-III. The fully coupled SCHISM-WWM-III modeling system has been applied to predict, simulate and hindcast typhoon-driven storm tides and waves [21][22][23][24]39], as well as long-term wave energy resources in the offshore waters of Taiwan [40][41][42].
For a successful storm surge, tide and wave modeling, the size of the computational domain must be large enough to accommodate the full typhoon; otherwise, the simulations would be affected by the boundary conditions [43,44]. The present study developed a computational domain spanning east longitudes from 105 • to 140 • and north latitudes from 15 • to 31 • (as shown in Figure 2). This large domain is composed of 540,510 nonoverlapping triangular elements and 276,639 unstructured grid points. A local-scale bathymetric dataset covering the area from east longitude 100 • to 128 • and north latitude from 4 • to 29 • with a spatial resolution of 200 m was provided by the Department of Land Administration and the Ministry of the Interior, Taiwan. The latest global-scale bathymetric product released from the General Bathymetric Chart of the Oceans (GEBCO), GEBCO_2020 Grid, was employed to incorporate a local-scale bathymetric dataset (as mentioned above) to construct the gridded bathymetric data in the SCHISM-WWM-III modeling system.
Eight main tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 and Q 1 ) extracted from a regional inverse tidal model (China Seas and Indonesia [45]) were utilized to generate the tidal elevation and horizontal velocity at the ocean boundaries of the SCHISM-WWM-III. The inverse barometric effect for tidal elevation at the boundary nodes was also considered. Since wave-generating typhoons were completely within the computational domain in the present study, open boundary conditions for the waves were not always required [42,46]. ric dataset covering the area from east longitude 100° to 128° and north latitude from 4° to 29° with a spatial resolution of 200 m was provided by the Department of Land Administration and the Ministry of the Interior, Taiwan. The latest global-scale bathymetric product released from the General Bathymetric Chart of the Oceans (GEBCO), GEBCO_2020 Grid, was employed to incorporate a local-scale bathymetric dataset (as mentioned above) to construct the gridded bathymetric data in the SCHISM-WWM-III modeling system. Eight main tidal constituents (M2, S2, N2, K2, K1, O1, P1 and Q1) extracted from a regional inverse tidal model (China Seas and Indonesia [45]) were utilized to generate the tidal elevation and horizontal velocity at the ocean boundaries of the SCHISM-WWM-III. The inverse barometric effect for tidal elevation at the boundary nodes was also considered. Since wave-generating typhoons were completely within the computational domain in the present study, open boundary conditions for the waves were not always required [42,46].

Parameterization of Depth-Induced Wave Breaking
Three widely used parameterizations for characterizing depth-induced wave breaking in nearshore shallow waters were adopted in the present study and are briefly reviewed in this section. Battjes and Beji [47] proposed that the total energy dissipation can be distributed over the wave spectrum in proportion to the spectral density based on extensive laboratory experiments; hence, the wave-induced wave breaking (Sbreak) in the wave action density function is calculated as follows:

Parameterization of Depth-Induced Wave Breaking
Three widely used parameterizations for characterizing depth-induced wave breaking in nearshore shallow waters were adopted in the present study and are briefly reviewed in this section. Battjes and Beji [47] proposed that the total energy dissipation can be distributed over the wave spectrum in proportion to the spectral density based on extensive laboratory experiments; hence, the wave-induced wave breaking (S break ) in the wave action density function is calculated as follows: where S break (σ,θ) is the wave-induced wave breaking in spectral space (σ, θ), D total is the total energy dissipation rate due to depth-induced wave breaking, E(σ,θ) is the wave energy in spectral space (σ, θ) and E total is the total wave energy.

BJ78 Model
In the BJ78 model (proposed by Battjes and Janssen [48]), the total energy dissipation due to depth-induced wave breaking, D tot_BJ78 is given as follows: where α BJ is the proportionality parameter, Q b represents the fraction of breaking waves, f is the averaged frequency and H max represents the maximum individual wave height and is defined as a proportion of the water depth in the local area (d): where γ is the wave-breaking index. The fraction of breaking waves, Q b is determined by the Rayleigh distribution truncated at an upper limit with H max (maximum wave height). Therefore, the fraction of breaking waves is implicitly expressed as follows: where H rms is the root-mean-square wave height.

TG83 Model
The TG83 model can be considered a variant of BJ78, and the total energy dissipation due to depth-induced wave breaking in the TG83 model (proposed by Thornton and Guza [49]), D tot_TG83 is formulated as follows: where B is the proportionality coefficient, H is the wave height and p b (H) is the fraction of breaking waves at each wave height. In which, p b (H) can be given as follows: where p(H) is the Rayleigh wave height probability density function and can be expressed as follows: The weighting function, W(H) is defined as in the TG83 model: where the calibration parameter n = 4, γ TG is the wave-breaking index for the TG83 model and is set to 0.42 according to Thornton and Guza [49,50].

CT93 Model
The total energy dissipation due to the depth-induced wave breaking implemented in the CT93 model (proposed by Church and Thornton [51]) is related to H rms through periodic linear bore theory, which can be given as follows: in which M is expressed as follows: 2.6. Parameterization of the Wave-Breaking Index 2.6.1. Constant Wave Breaking Criterion (γ) The parameter γ, given in Equation (4), is an adjustable coefficient that allows for the effects of the bottom slope on the waves and is one of the important parameters in the energy dissipation formulation. Miche [52]  index γ range from 0.60 to 0.83; however, γ = 0.73 has been used in the most operational third-generation wave models, which is taken as the averaged value over a larger data set [9]. In the present study, a value of γ = 0.78 was specified for the SWH simulation following the results from Chen et al. [21] and Hsiao et al. [22][23][24] if the wave-circulation model adopted a constant breaker index.

Wave Breaking Criterion Based on Local Steepness or Peak Steepness
The wave breaking process in shallow nearshore waters is affected by the seafloor profile (i.e., bottom slope) and incident wave steepness. Battjes and Stive [9] found that a hyperbolic tangent function could be used for predicting the SWHs in shallow-water coastal areas and a relationship between the breaking index γ and the wave steepness s 0 : in which s 0 can be expressed as follows: where H rms0 is the root-mean-square wave height in deep waters, λ 0 is the wavelength and is defined as λ 0 = gd/ f , f is the wave period and g is the acceleration of gravity. The wave breaking criterion in Equation (13) is regarded as the local steepness or peak steepness based on f .

Results
The wind fields during Super Typhoon Maria from 1-15 July in 2018 and Super Typhoon Lekima from 1-15 August in 2019 were extracted from the ERA5 reanalysis and were corrected by the direct modification method expressed in Equation (1). The hybrid typhoon winds derived from Equation (1) with various R trs values were imposed in the SCHISM-WWM-III modeling system to compare the performance of the resulting storm wave hindcast. The effect of the wave-breaking formulation and wave-breaking criterion on the simulation of the SWH in the shallow nearshore waters off northern Taiwan during the passage of Super Typhoons Maria and Lekima was investigated by conducting several designed model experiments, as listed in Table 2.

Validation for Typhoon-Driven SWHs with Various R trs
The inverse distance weighting method was employed to convert the hourly ERA5 and hybrid typhoon winds from the structured grid (at a horizontal resolution of 31 km) to the unstructured grid for the SCHISM-WWM-III modeling system. Comparisons of the SWH time series between model hindcasts using the different radii of the modified scale (R trs ) and the corresponding measurements are depicted in Figure 3a-c for the Fuguijiao (Figure 3a), Longdong ( Figure 3b) and Suao (Figure 3c) buoys during the passage of Super Typhoon Maria in 2018. The hindcasted peak wave heights are underestimated by 1.5 m and 2.0 m for the Fuguijiao and Longdong wave buoys, respectively, when the original ERA5 winds are used in the SCHISM-WWM-III modeling system but slightly overestimate the peak wave height within 0.5 m for the Suao buoy. Modified ERA5 winds with various R trs were imposed on the SCHISM-WWM-III modeling system to improve the performance of typhoon wave hindcasts. The hybrid winds with R trs = 3R max , R trs = 4R max , R trs = 5R max , R trs = 6R max and R trs = 7R max are called the H_3R max , H_4R max , H_6R max , H_6R max and H_7R max winds, respectively. As shown in Figure 3a-c, the hindcasted peak wave heights for all three wave buoys increased when using a larger R trs value. For instance, the hindcasted peak wave height was raised to 7.5 m for the Fuguijiao wave buoy by exerting H_4R max winds on the SCHISM-WWM-III modeling system (as shown in Figure 3a). Similar phenomena can be found in Figure 3b (for the Longdong wave buoy) and Figure 3c (for the Suao wave buoy). The hindcasted peak wave heights were always underestimated for the Fuguijiao wave buoy even though H_7R max winds were utilized. However, the SCHISM-WWM-III modeling system overestimated the peak wave height for the Suao wave buoy once the hybrid typhoon winds were used for the meteorological boundary conditions. Figure 4 illustrates the spatial distribution of the difference in the maximum hindcasted SWH between adopting the winds from H_4R max and the original ERA5 for Super Typhoon Maria in 2018. Significant differences can be detected along the track of Super Typhoon Maria, and the extents with differences exceeding 3.0 m occurred in the deep ocean. The spatial distributions of the difference in maximum hindcasted SWH between employing the winds from H_5R max and H_4R max , H_6R max and H_4R max and H_7R max and H_4R max are demonstrated in Figure 5a-c, respectively. The difference in maximum hindcasted SWH using different winds increased when R max was enlarged. The maximal differences were always distributed over the right side of Super Typhoon Maria in 2018, where the wind speed was highest (as shown in Figure 5a-c). The same validation process was conducted to verify the SWH hindcasts for Super Typhoon Lekima in 2019. Figure 6 presents the comparisons of the SWH time series between model hindcasts using the different R trs values and the corresponding measurements or the Fuguijiao (Figure 6a), Longdong ( Figure 6b) and Suao (Figure 6c) buoys during the period of Super Typhoon Lekima in 2019. The hindcasts from the use of the H_3R max and H_4R max winds are more satisfactory for all three wave buoys. The improvements in the hindcasted peak SWH for Super Typhoon Lekima in 2019 are obvious; for example, the difference in the maximum hindcasted SWH between inputting the winds from H_4R max and the original ERA5 were up to 5 m in the deep ocean (as shown in Figure 7). The spatial distributions of the difference in maximum hindcasted SWH between applying the winds from H_5R max and H_4R max , H_6R max and H_4R max and H_7R max and H_4R max for Super Typhoon Lekima in 2019 are shown in Figure 8a,c, respectively. Similar to the hindcasts of Super Typhoon Maria in 2018, maximal differences were also detected on the right side of Super Typhoon Lekima in 2019, where the wind speed was strongest and increased with the increase in R trs . According to the resulting storm wave hindcasts of Super Typhoons Maria in 2018 and Lekima in 2019, the hybrid wind field with R trs = 4R max , i.e., H_4R max winds, was adopted as the atmospheric forcing for the SCHISM-WWM-III modeling system to conduct a series of numerical experiments.
Maria in 2018, maximal differences were also detected on the right side of Super Typhoon Lekima in 2019, where the wind speed was strongest and increased with the increase in Rtrs. According to the resulting storm wave hindcasts of Super Typhoons Maria in 2018 and Lekima in 2019, the hybrid wind field with Rtrs = 4Rmax, i.e., H_4Rmax winds, was adopted as the atmospheric forcing for the SCHISM-WWM-III modeling system to conduct a series of numerical experiments.

Effect of the Wave-Breaking Formulation on the Typhoon-Driven SWH Simulation in the Surf Zone
To assess the effects of different wave-breaking formulations on the wave hydrodynamics in shallow nearshore waters, three depth-induced wave-breaking parameterizations, introduced in Section 2.5, were applied to hindcast the typhoon waves with the constant wave breaking criterion (γ) and same wind forcing (Rtrs = 4Rmax winds). As listed in Table 2, the designed numerical experiments refer to S_NO1 and S_NO3 for the BJ87 and

Effect of the Wave-Breaking Formulation on the Typhoon-Driven SWH Simulation in the Surf Zone
To assess the effects of different wave-breaking formulations on the wave hydrodynamics in shallow nearshore waters, three depth-induced wave-breaking parameterizations, introduced in Section 2.5, were applied to hindcast the typhoon waves with the constant wave breaking criterion (γ) and same wind forcing (R trs = 4R max winds). As listed in Table 2, the designed numerical experiments refer to S_NO1 and S_NO3 for the BJ87 and CT93 wave-breaking models, respectively, with a constant γ of 0.78. A wave-breaking index of 0.42, namely, γ GT , is specified for the TG83 model in the present study, which is labeled S_NO2 in Table 2. The spatial distribution of the difference in the maximum SWH between the scenarios of S_NO1 and S_NO2 and S_NO1 and S_NO3 during passage Super Typhoon Maria in 2018 are depicted in Figure 9. Figure 9a,b illustrate the differences in S_NO1 and S_NO2 and S_NO1 and S_NO2, respectively. Large storm waves were generated in the deep ocean and subsequently dissipated due to the decrease in water depth across a surf zone toward the shore. Although the surf zone usually lies in a shallow area where the water depth ranges from 5 to 10 m below sea level, nearshore areas with water depths greater than −20 m are shown in Figure 9 to effectively distinguish among the difference in hindcasted SWHs by different wave-breaking models. The surf zones along the north and northeast coasts of Taiwan are quite narrow because these areas are characterized by a very steep sloping seafloor. The differences in the maximum SWH between the S_NO1 and S_NO3 scenarios (Figure 9b) are more significant than those between the S_NO1 and S_NO2 (Figure 9a) scenarios. Two points, namely, P1 and P2, located in the north and northeastern shallow nearshore waters, were selected to compare the time series of hindcasted SWHs for the S_NO1, S_NO2 and S_NO3 scenarios. The comparison results are presented in Figure 10a for P1 and Figure 10b for P2. The differences between the three scenarios can only be found during the passage of Super Typhoon Maria in 2018, i.e., from midday on 10 July to midday on 11 July in 2018. The maximal difference in the hindcasted SWH is approximately 2.5 m between the S_NO1 and S_NO3 scenarios for P2 (as shown in Figure 10b). However, the maximal difference in hindcasted SWH is within 1.0 m between the S_NO2 and S_NO3 scenarios for both P1 and P2 (as shown in Figure 10a,b). Similar results are also shown in Figure 11 (spatial distribution) and Figure 12

Effect of the Wave-Breaking Criterion on the Typhoon-Driven SWH Simulation in the Surf Zone
Many parameterizations of the wave-breaking criterion have been implemented within wind wave spectral models to yield substantial improvements in model hindcasts, simulations and forecasts. Hence, three scenarios, called S_NO1, S_NO4 and S_NO5, exploiting the BJ87 wave-breaking model with Rtrs = 4Rmax wind forcing and different wavebreaking criteria (γ), were applied to hindcast storm waves in the shallow nearshore waters of northern Taiwan during the passage of Super Typhoon Maria in 2018 and Lekima in 2019. The spatial distributions of the difference in the maximum SWH hindcasted by S_NO1 and S_NO4 and S_NO1 and S_NO5 during the passage of Super Typhoon Maria in 2018 are illustrated in Figure 13a,b, respectively. The differences caused by different wave-breaking criteria γ are smaller than those caused by different wave-breaking formulations. Figure 14a,b demonstrate the time series of hindcasted SWHs for P1 and P2 using the S_NO1, S_NO4 and S_NO5 scenarios during the passage of Super Typhoon Maria in 2018. The maximal difference is within 1.0 m at P2 (as shown in Figure 14b) but is less

Effect of the Wave-Breaking Criterion on the Typhoon-Driven SWH Simulation in the Surf Zone
Many parameterizations of the wave-breaking criterion have been implemented within wind wave spectral models to yield substantial improvements in model hindcasts, simulations and forecasts. Hence, three scenarios, called S_NO1, S_NO4 and S_NO5, exploiting the BJ87 wave-breaking model with R trs = 4R max wind forcing and different wave-breaking criteria (γ), were applied to hindcast storm waves in the shallow nearshore waters of northern Taiwan during the passage of Super Typhoon Maria in 2018 and Lekima in 2019. The spatial distributions of the difference in the maximum SWH hindcasted by S_NO1 and S_NO4 and S_NO1 and S_NO5 during the passage of Super Typhoon Maria in 2018 are illustrated in Figure 13a,b, respectively. The differences caused by different wave-breaking criteria γ are smaller than those caused by different wave-breaking formulations. Figure 14a,b demonstrate the time series of hindcasted SWHs for P1 and P2 using the S_NO1, S_NO4 and S_NO5 scenarios during the passage of Super Typhoon Maria in 2018. The maximal difference is within 1.0 m at P2 (as shown in Figure 14b) but is less than 0.5 m at P1 (Figure 14a). The same phenomena were also detected in both the spatial distribution (as shown in Figure 15a,b) and time series of the hindcasted SWHs (as shown in Figure 16a,b) during the passage of Super Typhoon Lekima in 2019. Interestingly, the SWHs hindcasted by a constant γ (S_NO1 scenario) are usually smaller than those hindcasted by γ based on local steepness (S_NO4 scenario) and peak steepness (S_NO5 scenario).
Many parameterizations of the wave-breaking criterion have been implemented within wind wave spectral models to yield substantial improvements in model hindcasts, simulations and forecasts. Hence, three scenarios, called S_NO1, S_NO4 and S_NO5, exploiting the BJ87 wave-breaking model with Rtrs = 4Rmax wind forcing and different wavebreaking criteria (γ), were applied to hindcast storm waves in the shallow nearshore waters of northern Taiwan during the passage of Super Typhoon Maria in 2018 and Lekima in 2019. The spatial distributions of the difference in the maximum SWH hindcasted by S_NO1 and S_NO4 and S_NO1 and S_NO5 during the passage of Super Typhoon Maria in 2018 are illustrated in Figure 13a,b, respectively. The differences caused by different wave-breaking criteria γ are smaller than those caused by different wave-breaking formulations. Figure 14a,b demonstrate the time series of hindcasted SWHs for P1 and P2 using the S_NO1, S_NO4 and S_NO5 scenarios during the passage of Super Typhoon Maria in 2018. The maximal difference is within 1.0 m at P2 (as shown in Figure 14b) but is less than 0.5 m at P1 (Figure 14a). The same phenomena were also detected in both the spatial distribution (as shown in Figure 15a,b) and time series of the hindcasted SWHs (as shown in Figure 16a,b) during the passage of Super Typhoon Lekima in 2019. Interestingly, the SWHs hindcasted by a constant γ (S_NO1 scenario) are usually smaller than those hindcasted by γ based on local steepness (S_NO4 scenario) and peak steepness (S_NO5 scenario).

Discussion and Uncertainty
Hindcasting an accurate storm wave height is highly dependent on the accuracy of the typhoon wind fields. Reanalysis wind products underestimate typhoon winds, e.g., Çalışır et al. [53] evaluated the quality of ERA5 and CFSR (Climate Forecast System Reanalysis) winds and the contribution of reanalysis wind products to a wave modeling performance in a semi-closed sea. Their results revealed that ERA5 and CFSR tend to underestimate wind speeds, and ERA5 performs worse than CFSR at higher wind speeds (such as typhoon winds) and better at lower wind speeds. Thus, the utilization of hybrid winds through the superposition method (the combination of reanalysis wind products and parametric cyclone wind models) or the direct modification method (the combination of reanalysis wind products and the maximum wind speeds of typhoons from the best track

Discussion and Uncertainty
Hindcasting an accurate storm wave height is highly dependent on the accuracy of the typhoon wind fields. Reanalysis wind products underestimate typhoon winds, e.g., Çalışır et al. [53] evaluated the quality of ERA5 and CFSR (Climate Forecast System Reanalysis) winds and the contribution of reanalysis wind products to a wave modeling performance in a semi-closed sea. Their results revealed that ERA5 and CFSR tend to underestimate wind speeds, and ERA5 performs worse than CFSR at higher wind speeds (such as typhoon winds) and better at lower wind speeds. Thus, the utilization of hybrid winds through the superposition method (the combination of reanalysis wind products and parametric cyclone wind models) or the direct modification method (the combination of reanalysis wind products and the maximum wind speeds of typhoons from the best track data) is the most advantageous to consider storm wave hindcasting in both near-field and far-field regions of the typhoon's center. However, uncertainty remains to be clarified; that is, the radius of the modified scale (R trs ) cannot be formulated universally for each typhoon. Owing to the W Bmax is higher than W Emax , the area with stronger winds of a typhoon is extended as R trs increases (according to Equation (1)). This expansion allows the hindcasted SWHs to grow earlier and attenuate later (i.e., larger) than the measurements (as shown in Figures 3 and 6). Additionally, as shown in Figure 3b,c, the measured peak SWH at the Suao wave buoy occurred two hours earlier than that at the Longdong wave buoy. However, the occurrence time of hindcasted peak SWH at the Suao wave buoy coincided with that at the Longdong wave buoy. This phenomenon might be due to the low spatial resolution of the original ERA5 winds (at roughly 31 km), and the periphery circulation of a typhoon cannot be resolved with such a coarse spatial resolution. Although the wind field derived from R trs = 4R max is employed for designing a series of model experiments based on limited case studies, further studies are still needed to verify it.
Fully coupled ocean circulation and spectral wave numerical modeling systems are undergoing widespread use for all types of regional applications, such as operational predictions, wave climate evaluations, extreme storm waves and surge analyses. However, the wave-induced hydrodynamics during the period of typhoons in shallow nearshore waters simulated by these modeling systems remain uncertain. The scarcity of field observations for wave parameters in the surf zone to verify the modeling systems is one of the most important factors leading to uncertainty.
The SWH simulations in the surf zone are more sensitive to the various wave-breaking formulations than the various wave-breaking criteria. To reconfirm the result obtained from Section 3.2, the spatial distribution of the difference in the maximum SWH between the scenarios of S_NO1 (wave-breaking formulation of BJ87 with constant wave-breaking criteria) and S_NO2 (wave-breaking formulation of TG83 with constant wave-breaking criteria) and S_NO1 and S_NO3 (wave-breaking formulation of CT93 with constant wavebreaking criteria) in the surf zone (sea areas with water depths greater than −20 m are shown) of southeastern China for Super Typhoons Maria in 2018 and Lekima in 2019 are shown in Figures 17 and 18, respectively. As seen in Figures 17 and 18, the differences in maximum SWH between the S_NO1 and S_NO3 scenarios (Figures 17b and 18b) are higher than those between the S_NO1 and S_NO2 scenarios (Figures 17a and 18a) for both typhoons. Additionally, the surf zones with significant differences in maximum SWH resulting from the various wave-breaking formulations occur on the right side of the typhoons where the wind speeds are stronger. These findings are identical to the result derived from Section 3.2 and the reports from previous studies [54][55][56][57].
ing formulations than the various wave-breaking criteria. To reconfirm the result obtained from Section 3.2, the spatial distribution of the difference in the maximum SWH between the scenarios of S_NO1 (wave-breaking formulation of BJ87 with constant wave-breaking criteria) and S_NO2 (wave-breaking formulation of TG83 with constant wave-breaking criteria) and S_NO1 and S_NO3 (wave-breaking formulation of CT93 with constant wavebreaking criteria) in the surf zone (sea areas with water depths greater than −20 m are shown) of southeastern China for Super Typhoons Maria in 2018 and Lekima in 2019 are shown in Figures 17 and 18, respectively. As seen in Figures 17 and 18, the differences in maximum SWH between the S_NO1 and S_NO3 scenarios (Figures 17b and 18b) are higher than those between the S_NO1 and S_NO2 scenarios (Figures 17a and 18a) for both typhoons. Additionally, the surf zones with significant differences in maximum SWH resulting from the various wave-breaking formulations occur on the right side of the typhoons where the wind speeds are stronger. These findings are identical to the result derived from Section 3.2 and the reports from previous studies [54][55][56][57].

Summary and Conclusions
In this paper, the effects of wave breaking formulations and wave breaking criteria in hindcasting typhoon-driven storm waves are investigated for shallow nearshore waters off northern Taiwan. A fully coupled high-resolution, unstructured grid wave-circulation modeling system, SCHISM-WWM-III, with a large computational domain was applied to hindcast the wind waves caused by the passages of Super Typhoon Maria in 2018 and Super Typhoon Lekima in 2019. The ERA5 reanalysis product were merged with the maximum wind from the best track dataset, and the hybrid typhoon wind created and served as the meteorological conditions for the SCHISM-WWM-III modeling system through a direct modification technique. The hindcasted SWHs during the typhoon period were more sensitive to the radius of the modified scale, Rtrs and this phenomenon was more pronounced at the peak SWH. Rtrs equal to four times Rmax (the radius at the maximum typhoon wind speed) is an adequate radius of the modified scale for simulating the storm waves induced by Super Typhoon Maria in 2018 and Super Typhoon Lekima in 2019 via model validations. A series of numerical experiments were conducted using the SCHISM-WWM-III modeling system incorporated with modified typhoon winds to better understand the wave hydrodynamics in the shallow nearshore waters off northern Taiwan. The results derived from the designed numerical experiments reveal that the wave breaking formulations influence the hindcast of storm waves in the surf zone of northern Taiwan.

Summary and Conclusions
In this paper, the effects of wave breaking formulations and wave breaking criteria in hindcasting typhoon-driven storm waves are investigated for shallow nearshore waters off northern Taiwan. A fully coupled high-resolution, unstructured grid wave-circulation modeling system, SCHISM-WWM-III, with a large computational domain was applied to hindcast the wind waves caused by the passages of Super Typhoon Maria in 2018 and Super Typhoon Lekima in 2019. The ERA5 reanalysis product were merged with the maximum wind from the best track dataset, and the hybrid typhoon wind created and served as the meteorological conditions for the SCHISM-WWM-III modeling system through a direct modification technique. The hindcasted SWHs during the typhoon period were more sensitive to the radius of the modified scale, R trs and this phenomenon was more pronounced at the peak SWH. R trs equal to four times R max (the radius at the maximum typhoon wind speed) is an adequate radius of the modified scale for simulating the storm waves induced by Super Typhoon Maria in 2018 and Super Typhoon Lekima in 2019 via model validations. A series of numerical experiments were conducted using the SCHISM-WWM-III modeling system incorporated with modified typhoon winds to better understand the wave hydrodynamics in the shallow nearshore waters off northern Taiwan. The results derived from the designed numerical experiments reveal that the wave breaking formulations influence the hindcast of storm waves in the surf zone of northern Taiwan. The maximum difference in peak SWH could reach 2.5 m and 1.2 m for Super Typhoons Maria (2018) and Lekima (2019), respectively, when the wave-breaking formulations of BJ78 and CT93 were introduced to the SCHISM-WWM-III modeling system. Regarding the wavebreaking criterion on the hindcast of typhoon waves in the surf zone of northern Taiwan, compared with the wave-breaking formulation, the maximum difference in peak SWHs was relatively non-sensitive to the wave-breaking criterion. The maximum difference in peak SWH is only 0.5 m using the constant breaking criterion, the breaking criterion based on local steepness or the breaking criterion based on peak steepness. Another important finding is that the utilization of the BJ78 wave-breaking formulation usually underpredicts the typhoon-generated SWHs in shallow nearshore waters than other parametrizations (i.e., TG83 and CT93). In future research, it will be important to acquire measurements in the surf zone where wave hydrodynamics are more sensitive to wave-breaking formulations and criteria during the passage or landfall of typhoons.