Spectral Density Analysis for Wave Characteristics in Pohang New Harbor

The Pohang New Harbor (PNH), located at the Yongil bay in the northeastern part of Pohang city, South Korea, has experienced extreme wave hazards of about 3.0–5.0 m in elevation due to the seasonal swell from the far ocean. In this paper, both analytical and numerical studies are performed to investigate the wave-induced oscillations in an arbitrary shaped harbor with corner point consideration. By taking the consideration of the actual topography and bathymetry data, the boundary of PNH is constructed. Our theoretical model is based on the assumptions of inviscid, irrotational fluid, infinitesimal wave amplitude, and finally, the Helmholtz equation and its Weber’s solution. The numerical simulations are conducted to analyze the spectral density of the standing waves in PNH at eight respective synthetic record points. The simulation results are validated with real-time measurement data, which is obtained by wave heights and tide gauges at the specified record points within and outside the PNH. To improve the harbor’s design, a tactic such as building the breakwater at the entrance of the harbor is implemented and then spectral density is estimated in the modified geometry of the PNH. The consequential effects are proposed at the same time, suggesting the feasibility of the improvement measures.


Introduction
The Pohang New Harbor (PNH) was built to support the POSCO steel industry in South Korea, which has been a dominating industrial power since the 1970s. PNH obstructs the incoming waves from the East Sea (also known as the Sea of Japan), which is on the southeastern side of the Korean peninsula.
Occasionally, incident waves toward PNH have high amplitude when they are generated by typhoons, seasonal winds and currents, which are typically observed in Korean peninsula. Once such waves are entrapped inside PNH, only a small portion of the waves are diffracted and refracted by the exterior harbor boundary and structures. Usually, irregular ocean waves are often characterized by a wave spectrum that describes the distribution of wave energy (height) with respect to the frequency (wave period). The continuous frequency domain representation shows the power density variation of the wave frequency and is known as the wave amplitude energy density spectrum, or more commonly referred to as the wave energy spectrum. Incident waves with various angles in all possible directions are considered, but the components of such incident waves are strongly excited when their frequencies match the physical frequencies of the harbor because the wave energy provided from the exterior of the harbor tends to be fully accumulated and developed until it reaches a regime.
A variety of external forces can be responsible for substantial oscillations inside the harbor. But the most reported force is caused by incident waves generated by a tsunami or typhoon, which can cause the wave periods for anywhere from a few minutes to several hours. To understand the physical mechanism of the spectral density response in a harbor with arbitrary geometry, we analyzed the wave field on the free surface in the interior domain of the harbor. We observed that onshore waves are perfectly reflected, refracted and absorbed along the coast. A small portion is, however, diffracted through the entrance into the harbor and reflected repeatedly by the interior structures and boundaries.
Several previous studies in the literature have been done on the harbor resonance problem. Some of the theoretical studies on harbor oscillations are based on the constant water depth in and outside the harbor. MC NOWN (1952) analyzed the response of a circular harbor with a small entrance gap by assuming that the crest of the standing waves occurred at the entrance. The wave-induced oscillations in a rectangular harbor domain connecting to open sea was investigated by MILES and MUNK (1961) and IPPEN and GODA (1963). HWANG and TUCK (1970) presented a significant integral method, which involves the boundary wave function along the harbor boundary to evaluate the oscillations in the harbors of constant depth. The Weber's solution of the Helmholtz equation was introduced by LEE (1971) to analyze the wave oscillations in an arbitrary shaped harbor, which is very useful for practical implementation. BREAKHOFF (1976) used the mild slope equation with partial reflection conditions to evaluate the diffraction and refraction along the boundary of the harbor.
Some realistic models were also put forward through either analytical or numerical solutions. An integral equation with variable depth was formulated by MATTIOLI (1978) to determine the characteristics of wave field in the harbor basin while RAICHLEN and NAHEER (1976) demonstrated a finite difference scheme based on a shallow water equation to compensate for the wave oscillation of a harbor with variable depth and arbitrary harbor geometry. A finite element model to measure the interaction of a tsunami with the Hawaiian Island was introduced by HOUSTON (1978). CHOU and LIN (1986) applied the boundary element method to analyze the waveinduced oscillation in a harbor of arbitrary shape with a rigid wall in variable depths. A boundary element method was proposed by HAN (1993, 1994) to predict the wave height distribution in a harbor with variable depth and partially reflecting boundaries. LEE and WILLIAMS (2002) developed a boundary element modeling of multidirectional random waves in a harbor with partially reflecting boundaries.
In the port of Rotterdam, DE JONG and BATTJES (2004) investigated the triggering mechanism to identify the oscillation of the waves driven by wind and the atmospheric pressure. CHEN et al. (2004) presented the shelf resonance and edge wave effects at the bottom the sea, which play an important factor in harbor resonance. The comparison of several wave spectra of the random wave diffraction by a semiinfinite breakwater is given by LEE and KIM (2006). LEE et al. (2007) investigated diffracted wave field in a harbor for short waves using three dimensional numerical modeling. LEE et al. (2009) presented a numerical model based on the boundary element method for the prediction of wave field due to the diffraction of multidirectional random waves in a harbor with a long rectangular navigation channel.
However, most of the theoretical studies on harbor resonance mentioned above assumed that the water depth is constant in both interior and exterior regions of the harbor. In this paper, the mathematical model based on the Helmholtz equation with corner contribution is proposed to analyze the spectral density in the PNH. The present numerical model consists of two parts.
Firstly, we estimate the wave field under the small amplitude assumption in a harbor with an irregular geometry. The boundary integral equation method is used to solve the Helmholtz equation with corner contribution. In this method, the solution of a Helmholtz equation can be written in terms of an integral equation by using Green's identity formula. This is also known as Weber's solution of the Helmholtz equation. Furthermore, the boundary is discretized regularly and the wave field is numerically computed for the incident waves at various synthetic harbor record points. Using the computed fluid flow information such as surface wave displacement, the spectral densities have been analyzed at eight harbor record points in the PNH with respect to wave period (or frequency) for various directional incident waves. After the field survey and inspections inside the PNH, six record points are chosen to measure the real-time wave heights by using wave height and tide gauge (WTG). Consequently, our numerical simulation results are compared with the real-time data accessed by the POSCO steel corporation, Pohang, South Korea.
Secondly, on the basis of measurement data analysis, we have visualized that the incident waves propagate sharply into the PNH at the entrance corner with high amplification. A tactic has been introduced, which is supposed to add one additional breakwater at the sharp corner of the entrance to restrain the high impact incident waves. Thus the original PNH geometry is modified, in which the length of the harbor entrance is also increased. The numerical simulations are performed for the modified geometry of PNH afterwards. The amplification factor, which will be defined in a later section, and spectral densities of the modified PNH are compared with the original PNH at the given six record points inside the PNH. Overall, it can be concluded that the present numerical model is a competent tool to analyze the spectral densities of any harbor with irregular geometry for several practical applications.

Theoretical Development
In this paper, our theoretical background and its formulation are given to analyze the wave-induced oscillation inside the PNH. The theoretical solution was developed based on the Helmholtz equation for monochromatic incident waves, which propagate through the entrance inside the harbor with perfectly reflecting boundaries with sharp corners. The formulation is developed to solve the Helmholtz equation under the given assumptions and boundary conditions.

Model Description
In order to understand the physical phenomena of wave oscillation in a harbor, the following assumptions are considered on the dynamics of the fluid motion: (1) inviscid fluid, (2) irrotational flow, (3) infinitesimal wave amplitude, (4) relatively long wavelength compared to the water depth, and (5) lateral boundaries are perfectly reflective throughout the sea depth vertically. Wave terminology for small amplitude assumption is described by the Cartesian coordinates, which are employed with an origin set at the entrance of the harbor, as shown in Fig. 1. The x-axis is placed along the entrance of the harbor, the y-axis is directed toward the open sea and the z axis is directed vertically upwards. The exterior boundaries at the entrance are denoted by X 1 E 1 and X 2 E 2 . The domain of interest is divided into two regions. The first region is the open sea region X 1 , which includes the exterior coastal boundary X 1 X 2 . Another region is the bounded region X 2 surrounded by harbor walls oX, including the harbor entrance E 1 E 2 . The incident waves are diffracted, and reflected from the exterior boundary and the interior harbor walls oX. The wave field is denoted, respectively, in the open sea and the bounded region by g open and g. The governing equations are given in the bounded region X 2 and the open sea region X 1 , the matching boundary condition at the entrance is given as well (see Fig. 1). The outward and vector is expressed by n out and n in , respectively.

Governing Equations for Bounded Region
For the mathematical interpretation of the problem, the linear shallow water equation for mass conservation is given as follows (MEI et al. 2005): where f is the surface displacement, h is the water depth from sea floor, and the horizontal velocity vector is denoted by ũ ¼ ðu 1 ; u 2 Þ, respectively. For the momentum conservation, we can rewrite the velocity vector ũ in terms of surface displacement f through the following equation: where g is the acceleration of gravity. The total pressure is hydrostatic, which is expressed by where q is density of the water. From Eqs.
(1) and (2), the following equation can be obtained: which is a hyperbolic partial differential equation with variable coefficients. The incident waves are sinusoidal with respect to the time with radian frequency x. The surface displacement f and horizontal velocity ũ can be written as fðx; y; tÞ ffi gðx; yÞe Àixt ; ũðx; y; tÞ ffi uðx; yÞe Àixt ; ( where g is the wave field on the free surface. Then from Eqs. (1) and (5), we obtained the following expression and from Eq. (6), we get Given the constant water depth, the Eq. (4) is rewritten in the classical wave equation form, such that Finally, Eq. (7) turns into the Helmholtz equation where wave number k is defined by the dispersion relation for shallow water waves, i.e., k ¼ xðghÞ À1=2 : The fluid flow does not penetrate on the solid lateral boundary walls oX, and the bottom boundary condition on the sea floor is given by In our numerical model presented, the interior harbor boundary oX is considered as perfectly reflective. The equation r 2 g þ k 2 g ¼ 0 is solved by using the Green's identity formula. The zeroth order Hankel function of the first kind H ð1Þ 0 ðkrÞ is chosen to be the fundamental solution of this equation, which has singularity at origin as kr approaches zero. The wave field g at any point in the bounded region X 2 is expressed in the following integral form (LEE 1969(LEE , 1971 gðx; yÞ ¼ cðiÞ where ðx; yÞ and ðx 0 ; y 0 Þ are the variables, respectively, inside the harbor and on the boundary oX with function H 1 0 ðkrÞ is described in terms of the Bessel function of the first kind J 0 ðkrÞ and the second kind Y 0 ðkrÞ such that H 1 0 ðkrÞ ¼ J 0 ðkrÞ þ iY 0 ðkrÞ. The imaginary parameter c(i) is expressed as follows: where a is the corner angle.

Governing Equation for Open Sea Region
The complete wave field in the open sea region X 1 is expressed as the collection of three distinct parts: incident, reflected, and radiated wave fields. The incident waves are the incoming waves from the open sea region towards the harbor, and reflected waves are generated due to the reflection of incident waves from the interior and exterior boundaries of the harbor. The radiated waves are the emanating waves from the entrance toward the open sea region due to the reflection and diffraction of the incident waves. Thus the wave field in the open sea region is given by where g inc denotes the wave field due to incident waves, g rfl stands for the reflective wave field due to reflection of the exterior boundary and g R is the radiated wave field due to the waves dispersed into the open sea region through the entrance from the local topography of the harbor, respectively. The sinusoidal incident waves are defined as and the reflective wave field is expressed as where A is the wave amplitude of the incident wave and h i is the angle of incident wave. In the open sea region, the boundary condition on the straight coast line satisfies The normal flux vanishes along the straight coast for the radiated wave field, i.e., Furthermore, the radiation boundary condition due to the scattered wave emanating through the entrance towards the open sea must be satisfied on the far field, i.e., kR approaches 1, which can be written as follows and the matching boundary conditions on the harbor entrance E 1 E 2 are given by In the open sea region, the radiated wave field g R satisfies the following expression: The radiated wave field g R at any point in the open region X 1 is determined by using the Green's identity formula such that where oS is the exterior coastal boundary including entrance E 1 E 2 , r ¼ xÀ x 0 j j is the distance of point 1 ðkRÞ or on out . However, the point x i ¼ ðx i ; 0Þ and boundary point x 0 ¼ ðx 0 ; 0Þ lie on the xaxis, and the term or on out is equal to zero. Thus, the radiated wave field on the entrance E 1 E 2 can be written in the following form: From the matching boundary conditions, the wave field at the entrance is known. Furthermore, the wave field on the harbor boundary oX is determined by Eq. (11). Thus the wave field at any interior point in the bounded region X 2 is obtained by the following expression: where N is the total number of discrete segments on the harbor boundary oX including the entrance segments N 1, and Ds j is the length of jth segment. The surface displacement fðx; y; tÞ can be obtained at any point inside the harbor in terms of wave field gðx; yÞ. The amplification factor A.F. is defined as where f i is free surface displacement of the incident wave and g i is the incident wave amplitude, which is considered unity for convenience.

Spectral Analysis
The spectral analysis is commonly used to analyze the noise-like signal because it provides the Vol. 171, (2014) Spectral Density Analysis for Wave Characteristics frequency decomposition in a harmonic, which can be studied separately as a result. We herein have discussed the ocean surface waves, which are inherently nonlinear. The surfaces of sea waves are a collection of random waves of abundant wave frequencies (or wave periods) and wave lengths (wave length). How can we estimate the surface of sea waves? For that we need to do some simplification, which leads to the concept of wave spectrum (or spectral density). The spectrum provides the information about the distribution of the wave energy among the different wave frequencies (or wave periods) or wave length on the sea surface being analyzed. The spectral density can be defined as an image that identifies the relative wave energy presenting all frequencies or periods at a fixed location or region for a predefined time period, regardless of the energy's directional heading. In other words, it shows at which frequencies variations are strong and at which frequencies variations are relatively weak.
In order to compute the spectral density function from a discrete time record, different numerical methods have been introduced in the society of computing fluid. Among them, the most reliable and efficient method is called fast Fourier transformation (FFT), which is an algorithm for calculating the discrete Fourier transformation (DFT). It calculates autocorrelation function and then transforms it. Our DFT calculation of wave spectra is shown in the next section.

Calculation of Wave Spectrum
The free surface displacement fðx; y; tÞ is numerically calculated at any interior point in the harbor region with respect to the time series data. All the regular waves of the sea surface are estimated numerically, assuming that the direction of the incident waves propagating towards the entrance is known. Thus the discretized form of the free surface displacement f at record point (x, y) inside the harbor is where Dt is the sampling interval between two waves and N is the total number of waves. The length of the record wave period is given such that L t ¼ NDt. To calculate the wave spectrum, the discretized Fourier transform Z n of a surface wave elevation record f j with time t j is expressed as the following: By using the FFT, these equations can be summed up and written in the simple spectrum S n (COOLEY et al. 1970) . . .; ðN=2 À 1Þ; The amplitude is calculated by In the numerical calculation, the FFT algorithm of MATLAB has been used. Thus, the Eq. (27) is rewritten as The numerical calculation of Eq. (30) can be interpreted as the product of the matrix (the exponential part) multiplied by the vector (surface wave elevation). The efficiency of the FFT depends on the proper decomposition of such a matrix and other applicable permutations that slightly reduce the number of operations needed (BRIGHAM 1988). The possible decompositions depend on the length of the vector. Therefore, the execution time and the required memory for the FFT depend on the length of f (in terms of number of data points). The FFT runs the fastest for the record point in which the number of data points is a power of two. In addition, amplitude is estimated as corresponding to the interval frequency Df ¼ 1=Dt. The wave spectrum (or spectral density) is calculated by where f i is the ith frequency recorded with ith amplitude A i , and sampling frequency interval is considered as Df ¼ f i À f iÀ1 .

Location of Port and the Record Points
PNH is located in Pohang City, in the southeast of South Korea. Its exact location in terms of longitudinal and latitudinal coordinates is 36.0071-36.0432 and 129.3800-129.4377, respectively. It was constructed to serve one of the largest steel corporations, named POSCO, on its northwest margin. The loading capacity of PNH is 47 million tons per year with maximum ship size 250,000 DWT (KWAK et al. 2008). The bird's-eye view image of PNH shown in Fig. 2 was taken by the Global Mapping Tool software package. The location of the eight record points W1-W8 in the PNH were used to analyze the spectral density (see Fig. 2). The extreme wave oscillations of around 3-5 m in height have been  Vol. 171, (2014) Spectral Density Analysis for Wave Characteristics observed at these record points during seasonal weather conditions. The behavior of the incident wave directions was also noted at the entrance during real-time measurement data analysis. The location of these record points was investigated during the measurement and found to be among the most sensitive wave oscillation positions in the PNH, because the maximum surface wave height of the long waves (or long wave height) was observed at these record points during the seasonal weather conditions. In Fig. 2, the realistic view of the exact locations of the record points is clearly visualized. Hence, the exact locations of eight record points from W1 to W8 are given in terms of longitudinal and latitudinal coordinates in Table 1. The six synthetic interior record points W3, W4, W5, W6, W7 and W8 are basically the inner port stations where the moored ships are loaded and unloaded. The locations of record points W1 and W2 are placed as reference points outside the harbor.

Spectral Density Analysis for 100-min Time Series Data
The numerical simulations are conducted to analyze the spectral density at eight record points from W1 to W8 in the PNH with respect to the wave period. The longest wave period (or frequency) is considered as L t = 100 min with sampling interval Dt ¼ 30 s: Therefore, in our numerical process, we estimate the spectral density for N = 200 discrete samples at each record point from W1 to W8. In Fig. 3, the spectral density at the record points W1, W2, W3 and W4 is shown with respect to the wave period. Record points W1 and W2 are placed outside the PNH as the reference points to examine the incoming waves with particular directions. In the spectral density analysis, the resonance peaks occurred at the wave periods T 1 = 4-5 min, T 2 = 8-9 min, T 3 = 20-30 min and T 4 = 60-90 min for the record points W1, W2, W3 and W4, respectively, as shown Fig. 3. Similarly, the spectral density at record point-5 (W5), record point-6 (W6), record point-7 (W7) and record point-8 (W8) is shown in Fig. 4 with respect to wave period inside the PNH. The resonance peaks occurred approximately at the wave period (or frequency) T 1 = 4-5 min, T 2 = 8-9 min, T 3 = 20-30 min and T 4 = 60-90 min for the record points from W5 to W8. However, the peak heights of the spectral density vary corresponding to each record point in the PNH, suggesting that the wave energy distribution varies with respect to frequency as well as the position of record points.
The incoming waves with periods from T 1 to T 4 have high amplification inside the PNH. Essentially, they are hazardous to the moored ship and coastal structure of the harbor. In another words, the incident waves with such frequencies (or wave periods) generate higher resonance inside the PNH. There are some relatively smaller peaks in the spectral density graphs, which depict the local resonance due to the complex geometry of the harbor with sharp corners. The spectral density at various record points in the interior region of the PNH provides the crucial information about the intensity of incoming waves with the corresponding frequencies (or wave periods) and incoming directions. Usually in the coastal region like a harbor, it is noticed that the wave period range of the ocean waves is from approximately 1 to 10 min. In the next section, the sampling interval size is reduced for the detailed spectral density analysis.

Spectral Density Analysis for 10-min Time Series Data
For the small sampling interval, the spectral density at eight record points (W1-W8) is numerically estimated with respect to wave period. The length of the record wave period is reduced to L t = 10 min with sampling interval Dt ¼ 10 s: The sampling interval (Dt ¼ 10 s) is refined in order to analyze the wave energy distribution at closely distributed frequencies (or wave periods) for the initial range (until 10 min). In Fig. 5, the numerical simulation results for the spectral density at record points W1, W2, W3 and W4 are displayed with respect to wave period with sampling interval Dt ¼ 10 s: The resonance peaks are determined at wave periods T 1 = 4-5 min and T 2 = 8-9 min at each record point in the PNH.
Similarly, the spectral density analysis is shown in Fig. 6 at record point-5 (W5), record point-6 (W6), record point-7 (W7) and record point-8 (W8) with respect to wave period with sampling interval Dt ¼ 10 s: The spectral density provides the distribution of wave energy among the refined frequencies at these record points. Thus the resonance peaks at these record points (W5-W8) occurred at the frequencies (or wave periods) T 1 = 4-5 min and T 2 = 8-9 in the PNH.
In Figs. 5 and 6, the fluctuation of spectral density can be visualized at various record points in the PNH for small sampling interval Dt ¼ 10 s: In other words, Spectral density with respect to wave period (min) is given at record point-1 (W1), record point-2 (W2), record point-3 (W3) and record point-4 (W4) with sampling interval Dt ¼ 10 s: The resonance peaks are obtained at the frequencies (or wave periods) T 1 = 4-5 min and T 2 = 8-9 min, which are encircled the relative wave energy distribution is determined at all frequencies (or wave periods) at a fixed location in the PNH. Moreover, the incident waves corresponding to the resonance peaks in the spectral density graphs have high amplification in the harbor and these waves generate extreme wave oscillation in the harbor. Usually we have noticed the shorter wave periods in the harbor from 10 s to 10 min. Thus the incoming waves with short periods produced higher resonance in the harbor.

Validation of the Numerical Results with Onsite Measurement
The wave heights were measured at record points (W1-W8) in the PNH with respect to wave periods. One directional wave rider buoy is deployed at record point W1 outside the harbor and a pressure-type WTG is installed at record points W3 through W8 inside the harbor. The wave heights were measured at record points W1 and W2 outside the PNH and six record points inside the PNH. For qualitative comparison, we also took advantage of some real-time field observation data from the PNH, which were provided by the POSCO steel plant. The onsite measurements and instruments were deployed as shown in Fig. 7.
For the verification of the numerical model, our numerical calculation results are compared with the in situ real time measurement data at record point-3 (St. W03) and record point-7 (St. W07). In Fig. 8, the spectral density is plotted with respect to wave period based on in situ measurement data at record point-2 at the wave periods T 1 = 5 min, T 2 = 8 min, T 3 = 21-32 min and T 4 = 60-84 min, respectively. In our numerical simulation, the spectral density is analyzed at record point-3 (St. W03) and record point (St. W07) with respect to the wave periods shown in Fig. 9. In addition, we evaluate the first few resonance peaks that occurred at various frequencies (or wave periods) T 1 = 4-5 min, T 2 = 8-9 min, T 3 = 20-30 min and T 4 = 60-90 min. Thus our simulation results can be compared physically with the in situ measurement data, implying a good agreement between them. In another words, our numerical model almost predicts the resonance modes, which are consistent with the observation data. On the basis of comparison, the numerical accuracy of the analysis of spectral density is validated. Therefore the validation of simulation results show the efficiency and stability of our numerical model, which can implemented on any irregular domain with complex topography to analyze the spectral density at any fixed location in the harbor. On the basis of onsite measurement analysis, a tactic has been introduced to improve the geometry of the original PNH, which is explained in the next section.
6. Improved PNH Design and the Numerical Results

Modified PNH Design
During the in situ measurement data analysis, we observed the direction of the incident waves at the entrance. The incoming waves propagate sharply at the corner E 1 on the entrance of the harbor with high amplification. According to the analysis results mentioned in previous sections, a tactic is proposed to restrain such incoming waves at the entrance, i.e., a breakwater is attached at the entrance as shown in    Fig. 10. The original geometry of the PNH is modified by adding a breakwater at the entrance in our numerical model. The idea behind this modification is to slowly and deliberately obstruct the direct incoming waves from the northeast side of the PNH. The length of the breakwater in the modified PNH is given as 424 m to restrain the incoming waves at the entrance. As a result, the present numerical model is employed on the modified geometry of the PNH to determine the amplification factor and the spectral density. The numerical discretization for the modified PNH has been retrieved to employ the numerical scheme at six record points inside the PNH. In the next section, the simulation results in the modified PNH are compared with the original PNH to analyze the response of the incident waves at six record points (W3-W8) inside the PNH.

Numerical Results of the Modified PNH Design
The numerical simulations were carried out on the modified geometry of the PNH to analyze the amplification factor at six record points (W3-W8).
The boundary was discretized precisely, and the segments at the corner were refined. In the modified geometry of the PNH, the length of entrance is increased compared to that of the original PNH. In order to analyze the response of the incident waves inside the harbor, the amplification factor is defined as the ratio of wave amplitude at any position inside the harbor to the amplitude of incident waves. Thus, the amplification factor in the modified PNH is compared with the original PNH as shown in Fig. 11 at six record points from W3 to W8. On the basis of comparison, it is suggested that the amplification factor is reduced in the modified PNH. The resonance peaks heights are reduced in the modified PNH compared to the original PNH at all six record points inside the PNH. It is concluded that the tactic for improvement such as adding the breakwater at the entrance is significantly effective to reduce the resonance in PNH.
In our study, the spectral density has been estimated for the original and modified PNH with respect to wave period with sampling interval DT ¼ 30 s: The length of wave frequency (or wave periods) is considered L t ¼ 100 min: The comparison of spectral density between the original and modified PNH at record point-3 (W3) through record point-8 (W8) is shown in Fig. 12. The spectral density at various record points for the modified PNH has declined compared to the original PNH, which indicates that the implementation of the breakwater can reduce the harbor oscillation significantly. The implementation of some basic tactics such as remolding the geometry of the harbor, adding breakwaters and widening the harbor entrance can be effective to improve the resonance problem in the harbor. The spectral density provides efficient information about the wave energy distribution at key locations in a harbor with complex geometry. Therefore our proposed method can efficiently predict the intensity of the incident waves with various wave periods in the harbor.

Conclusion and Discussion
The spectral density analysis provides crucial information about the energy distribution (wave height) with respect to the frequency variation at key locations in and outside the PNH. The resonance peaks in the spectral density is obtained at different frequencies (or wave periods). Therefore the incident waves with such frequencies have high amplification, which are hazardous to costal structure and moored vessels in the harbor. Our proposed numerical simulation results are validated through the comparison with the in situ measurement data. In order to analyze the spectral density in a harbor, using a realistic physical model is expensive and time intensive. However, the model proposed in this paper is easily implemented on any arbitrary shaped harbor, effectively and efficiently. The potential tactic is put forward, such as adding a breakwater at the entrance, and the original PNH design is modified accordingly. In the modified PNH, the amplification factor and spectral density were reduced at the interior record points from W3 to W6, suggesting that the incident wave amplitude and direction play an important role in reducing the resonance in the harbor. This model can be utilized as a reliable engineering tool for harbor planning and designing to protect the interior structure of a harbor and moored ships.
For more practical problems, the scope of our discourse here is sufficient, however, far more accurate methods are required, such as taking into consideration the nonlinearity to analyze harbor   resonance. To achieve higher precision and account for these additional parameters, more quantitative numerical models are needed to obtain greater accuracy and a more genuine representation of realistic, dominant harbor conditions.

Acknowledgments
The funding for this work was supported by the Department of Mathematics, POSTECH and BK-21 (Brain Korea-21), Pohang, South Korea. It was also