Wave Characteristics over a Dual Porous Submerged Breakwater Using a Fully Nonlinear Numerical Wave Tank with a Porous Domain

: This study developed a two-dimensional fully nonlinear numerical wave tank (FN-NWT) to examine the nonlinear interaction between waves and dual submerged porous structures. Using the FN-NWT, not only reﬂection and transmission coefﬁcients, but also wave deformation/force depending on porosity were investigated. The FN-NWT was developed using the boundary element method (BEM), and consisted of a ﬂuid domain and a porous medium domain. Darcy’s law or the non-Darcy (Forchheimer) ﬂow equation were applied to the ﬂow passing through the porous domain. The wave reﬂection coefﬁcient of the porous submerged structures agreed well with the given experimental data when using Forchheimer ﬂow boundary conditions. Excessive attenuation of the transmitted wave occurred when Darcy’s condition was employed. The difference in each coefﬁcient due to the spacing of the submerged structure was reduced in the porous structure compared with the non-porous structure. The difference according to the incident wave height was clearly revealed in the transmission coefﬁcient. The developed dual-domain FN-NWT can be applied to investigate the nonlinear interaction between waves and porous structures as a ﬁrst-cut design tool.


Introduction
Considerable technical efforts are being made to minimize the impact of incident waves, protect coastal facilities, and mitigate coastal erosion. Many researchers have examined the effects of submerged breakwaters. Wave energy and wave loads can be changed due to submerged objects such as reefs [1]. Moroever, submerged structures such as Articulated Concrete Block Mattresses are employed to maintain the ecosystem while minimizing coastal erosion [2]. Structures with porosity are also utilized as breakwaters. The main advantage of porous breakwaters over non-porous breakwaters is greater wave attenuation [3]. They also have additional benefits, such as better water circulation and calmness upstream [4,5]. The scouring depth decreases because the porous structure reduces wave energy [6].
Studies on the characteristics of surface waves over submerged structures have mainly been conducted experimentally. To physically make submerged porous structures, a plastic ball [7] and a wooden caisson model [8] were used. In addition, a laboratory experiment was conducted on the porous structure installed on a sloping seabed as well as a flat seabed [9].
In computational fluid dynamics (CFD) techniques, the porous structure effect has been considered in numerous ways. The modified Navier-Stokes equation including the resistance of the porous medium was defined as the governing equation, without directly linear wave theory was used in the frequency domain, and time domain-based analysis was rare.
In this study, both fluid and porous regions were numerically modeled in the time domain using nonlinear BEM. The present authors have already performed a time domain analysis based on a linear wave model under permeable seabed conditions [35] and simulated nonlinear wave propagation over a porous sloped seabed [36]. However, those studies were for small seabed permeability (such as soil), without porous breakwaters. Therefore, in this study, we developed a fully nonlinear numerical wave tank (FN-NWT) and used a two-domain technique to numerically analyze nonlinear wave propagation over dual submerged breakwaters with high permeability.
Using the FN-NWT, we analyzed waves reflected and transmitted by a dual submerged porous structure. The computational domain consists of the fluid domain and the porous domain. The Laplace equations for the velocity potential (fluid domain) and pore pressure (porous domain) are applied as the governing equations, respectively. A three-step boundary value problem was solved at each time step in the time domain analysis to calculate the interaction between the fluid domain and the porous domain. Darcy flow conditions are limited for high flow velocity (large permeability) [37]. Therefore, non-Darcy (Forchheimer) flow was also applied as a boundary condition at the fluid-porous interface, and the differences in wave characteristics between the Darcy and non-Darcy flows were compared. The numerical model was verified by comparing the wave reflection coefficient by the submerged dual porous structure with existing research data. The wave reflection, transmission, and dissipation coefficients of non-porous and porous structures were compared. The wave characteristics for various incident wave heights and distances between the structures were also compared.

Methods
The boundary value problem was solved by numerically modeling a dual computational domain consisting of fluid and porous domains, which are briefly described below.

Fluid Domain
The viscosity of fluid can be neglected in wave propagation, and the fluid particle motion was assumed to be potential flow. Therefore, the Laplace equation (Equation (1)) for the velocity potential (φ) was used as the governing equation in the fluid domain. This governing equation was converted to the boundary integral equation (Equation (2)) using Green's 2nd identity.
where α is the solid angle and equals 0.5 when the node is on the boundary; G ij = −(1/2π)lnR ij is the green function, which can be defined by the distance between the source and field points (R ij ) in the fluid domain (Ω F ). The boundary conditions of each boundary required to solve the boundary value problem in the fluid domain are as follows: the velocity profile of the incident wave was used for the left-hand side boundary condition. In fully nonlinear calculations, spatial wave modulation can occur due to the mismatch in particle motion on the free surface near the incident wave boundary, but only if linear wave particle velocity is employed as the boundary condition. To mitigate this spatial modulation of free surface waves, Stokes' 2nd-order wave was applied as the incident wave boundary condition.
where g is the gravitational acceleration; A is the incident wave amplitude; k is the wave number (=2π/λ, λ = wavelength); ω is the wave frequency; h is the water depth; n x means the x-direction component of the normal vector. To generate a stable incident wave and minimize any transient effects, a ramp function (Ramp) was applied at the beginning of the time simulation.
where T is the incident wave period. For free surface boundary conditions, the mixed Eulerian-Lagrangian (MEL) method was applied to treat the nonlinear free surface boundary conditions. The free surface boundary conditions were rearranged through the total derivative ( δ δt = ∂ ∂t + → v · ∇) and material node approach (full Lagrangian approach), assuming the node velocity ( → v ) is equal to fluid particle velocity ( (5) and (6)). Therefore, the node points were updated to match the fluid particle movement. In addition, to minimize the wave reflection from the end of the computation domain (right-hand side wall), artificial damping terms were added to Equations (5) and (6), such as the last term of RHS (right-hand side).
δφ δt δη δt where η is the wave elevation and µ m (m = 1, 2) represents an artificial damping coefficient, defined as follows: where l d is the length of the end-damping zone; x fd means the x-coordinate of the end point of the frontal-damping zone; and x d means the x-coordinate of the starting point of the end-damping zone. The damping coefficients have a relationship, µ 02 = kµ 01 . In this study, the submerged breakwater was placed in the computational domain. Thus, the propagating wave can be reflected by the submerged breakwater, and the reflected wave propagates in the opposite direction of the incident wave. Therefore, when the incident wave entered the computational domain, the free surface artificial damping scheme was applied to the difference between the total and incident waves in front of the left-hand side wall in the computational domain to prevent the dissipation of incident waves as follows: δφ δt where φ* and η* mean incident wave potential and elevation, excluding the reflected wave components. The incident wave can be obtained from Stokes' 2nd-order wave theory. From the last terms of RHS in Equations (8) and (9), only the reflected waves can be attenuated. The symbol µ fm (m = 1, 2) represents a frontal artificial damping coefficient, defined as follows: where l fd is the length of the frontal-damping zone. In this study, the length of the artificial damping zone was set to l d = l fd = 1λ for kh < 1.0 and l d = l fd = 2λ for kh ≥ 1.0. Using the above free surface boundary conditions and the Runge-Kutta 4th-order scheme, the wave elevation and velocity potential on the free surface can be updated at each time step in the Lagrangian manner. Moroever, to avoid the possible saw-tooth instability of the free surface wave, which can occur during the fully nonlinear material node approach and the corresponding node update, the Chebyshev five-point smoothing scheme [38] was applied.
The porous boundary condition was classified into Darcy flow (Equation (11)) and Forcheimer (non-Darcy) flow (Equation (12)). As the flow rate through the porous domain increases, the flow velocity and pressure gradient have a nonlinear relationship [39].
where p s is the pore pressure.
→ u is the discharge velocity and is equal to (−∂φ/∂n) at the interface between the fluid and the porous domain. K is the permeability coefficient and µ is the dynamic viscosity of fluids, equaling 10 −3 Pa·s in this study. β F is the non-Darcy factor and is expressed as β F = c/ √ K (c: Forchheimer coefficient). ρ is the water density. The boundary condition for the impermeable bottom and end wall is defined as follows:

Porous Domain
Assuming that the porous domain satisfied the continuity equation (Equation (14)) and Darcy's law (Equation (11)), the Laplace equation for the pore pressure is used as the governing equation in the porous domain (Equation (15)). The boundary integral equation for pore pressure is Equation (16).
In the porous domain, the wave dynamic pressure was applied to the interface boundary condition between the fluid and the porous domain.
An impermeable boundary condition was applied to the side wall and the bottom boundaries.

Numerical Calculation Process and Coupling of Fluid Domain and Porous Domain
Using the fluid and porous medium model described in the above sections, the waveporous structure interaction was calculated and nodes were rearranged at each time step. The flowchart of the total process is presented in Figure 1

Numerical Calculation Process and Coupling of Fluid Domain and Porous Domain
Using the fluid and porous medium model described in the above sections, th porous structure interaction was calculated and nodes were rearranged at each ti The flowchart of the total process is presented in Figure 1.   Figure 1. A three-step calculation scheme [36] was used to calcu coupling of the fluid and porous domains. First, the boundary value problem (B the time derivative of velocity potential (∂ϕ/∂t) in the fluid domain is solved, and tained values are used as the interface boundary condition (Equation (17)) in the domain. The boundary condition for ∂ϕ/∂t on the interface between the fluid and rous domain is assumed to be Equation (19). Then, the time derivative of ∂ϕ/∂n obtained from sub-step values of ∂ϕ/∂n in each time step using the finite differe mula [40] as shown in Equation (20).   Figure 1. A three-step calculation scheme [36] was used to calculate the coupling of the fluid and porous domains. First, the boundary value problem (BVP) for the time derivative of velocity potential (∂φ/∂t) in the fluid domain is solved, and the obtained values are used as the interface boundary condition (Equation (17)) in the porous domain. The boundary condition for ∂φ/∂t on the interface between the fluid and the porous domain is assumed to be Equation (19). Then, the time derivative of ∂φ/∂n can be obtained from sub-step values of ∂φ/∂n in each time step using the finite difference formula [40] as shown in Equation (20).
∂ ∂t where φ n (l) = ∂φ (l) /∂n and the superscript (l) mean each sub-step value in the Runge-Kutta 4th order. Second, the BVP for the pore pressure in the porous domain is solved, and the resulting values are used as the interface boundary conditions (Equations (11) and (12)) in the fluid domain. At this time, in the non-Darcy flow boundary conditions, the discharge velocity ( → u ) is calculated using an iterative method. Finally, the BVP for the velocity potential is solved again, and the free surface elevation affected by the porous domain can be calculated.
where ϕn (l) = ∂ϕ (l) /∂n and the superscript (l) mean each sub-step value in the Runge-Kutta 4th order. Second, the BVP for the pore pressure in the porous domain is solved, and the resulting values are used as the interface boundary conditions (Equations (11) and (12)) in the fluid domain. At this time, in the non-Darcy flow boundary conditions, the discharge velocity ( u ) is calculated using an iterative method. Finally, the BVP for the velocity potential is solved again, and the free surface elevation affected by the porous domain can be calculated. In developing the FN-NWT, the advantage of solving the BVP for the time derivative of velocity potential independently in addition to the BVP of velocity potential itself was well explained as the acceleration potential method in [41]. Alternatively, the time derivative of the velocity potential can be obtained numerically through the backward difference formula, which is not only less accurate but also may cause numerical instability during the time-marching process.

Validation of Numerical Results
The numerical model of this study was compared with experimental data for verification [27]. Figure 3 presents a schematic diagram of a dual porous breakwater. The detailed setup conditions of the numerical model are as follows: water depth (h) = 0.8 m; incident wave height (HI) = 0.04 m; wave period (T) ranges from 1.29 s to 3.73 s. The submerged breakwater was set up under the following conditions: bottom width (Wb) = 2h; top width (Wt) = 0.5h; height (Hb) = 0.5h; distance between structures (Lb) = 3.6 m. The submerged structure was a tetrapod model in an experimental study [27] and its effective diameter (d50) and porosity (n) were 0.076 m and 0.5 m, respectively [28]. In developing the FN-NWT, the advantage of solving the BVP for the time derivative of velocity potential independently in addition to the BVP of velocity potential itself was well explained as the acceleration potential method in [41]. Alternatively, the time derivative of the velocity potential can be obtained numerically through the backward difference formula, which is not only less accurate but also may cause numerical instability during the time-marching process.

Validation of Numerical Results
The numerical model of this study was compared with experimental data for verification [27]. Figure 3 presents a schematic diagram of a dual porous breakwater. The detailed setup conditions of the numerical model are as follows: water depth (h) = 0.8 m; incident wave height (H I ) = 0.04 m; wave period (T) ranges from 1.29 s to 3.73 s. The submerged breakwater was set up under the following conditions: bottom width (W b ) = 2h; top width (W t ) = 0.5h; height (H b ) = 0.5h; distance between structures (L b ) = 3.6 m. The submerged structure was a tetrapod model in an experimental study [27] and its effective diameter (d 50 ) and porosity (n) were 0.076 m and 0.5 m, respectively [28]. First, a convergence test was performed by increasing the number of nodes for the non-porous structure, as shown in Figure 4. The time step was set to Δt = T/64. From Figure 4, if the number of nodes per wavelength is 25 or more, the wave elevation appears to converge. Therefore, in this study, the calculation was performed by setting the number of nodes to Δx = λ/30. First, a convergence test was performed by increasing the number of nodes for the non-porous structure, as shown in Figure 4. The time step was set to ∆t = T/64. From Figure 4, if the number of nodes per wavelength is 25 or more, the wave elevation appears to converge. Therefore, in this study, the calculation was performed by setting the number of nodes to ∆x = λ/30. First, a convergence test was performed by increasing the number of nodes fo non-porous structure, as shown in Figure 4. The time step was set to Δt = T/64. From ure 4, if the number of nodes per wavelength is 25 or more, the wave elevation appea converge. Therefore, in this study, the calculation was performed by setting the num of nodes to Δx = λ/30.  In the reference experimental paper [27], a tetrapod was used to set up the po breakwater. To obtain the permeability coefficient, the suggested relation between pe ability coefficient (K) and the porosity (n) for the tetrapod was used [42].  (21), the permeability coefficient was set to K = 2.039 × 10 −7 m 2 in study. Additionally, the Forchheimer coefficient was obtained using the empirical form [43]. This value was a function of effective diameter (d50) and porosity (n), as follows: Equations (12), (22), and (23), the non-Darcy factor (βF) of the porous structures ca obtained as Equation (24).  In the reference experimental paper [27], a tetrapod was used to set up the porous breakwater. To obtain the permeability coefficient, the suggested relation between permeability coefficient (K) and the porosity (n) for the tetrapod was used [42].
From Equation (21), the permeability coefficient was set to K = 2.039 × 10 −7 m 2 in this study. Additionally, the Forchheimer coefficient was obtained using the empirical formula [43]. This value was a function of effective diameter (d 50 ) and porosity (n), as follows: from Equations (12), (22), and (23), the non-Darcy factor (β F ) of the porous structures can be obtained as Equation (24).
where h p means the piezometric head. From the above equations, the non-Darcy factor was set to β F = 109.54/m (c = 0.0495). Figures 5 and 6 present snapshots of the incident wave elevations over non-porous (rigid) and porous structures. In the case of kh = 1.0, the reflected wave was almost the same regardless of the porosity of the structure, but the transmitted wave was greatly reduced by the porous structure. In Figure 6 (kh = 1.2), both reflected and transmitted wave elevations were reduced by the porous structure. The wave attenuation in the model with the Darcy flow boundary condition was greater than that under the non-Darcy flow condition. The transmitted waves passing over non-porous structures look more deformed from the linear sinusoidal form than those over the porous structures.
(rigid) and porous structures. In the case of kh = 1.0, the reflected wave was almos same regardless of the porosity of the structure, but the transmitted wave was great duced by the porous structure. In Figure 6 (kh = 1.2), both reflected and transmitted elevations were reduced by the porous structure. The wave attenuation in the model the Darcy flow boundary condition was greater than that under the non-Darcy flow dition. The transmitted waves passing over non-porous structures look more defo from the linear sinusoidal form than those over the porous structures. (rigid) and porous structures. In the case of kh = 1.0, the reflected wave was almos same regardless of the porosity of the structure, but the transmitted wave was greatl duced by the porous structure. In Figure 6 (kh = 1.2), both reflected and transmitted w elevations were reduced by the porous structure. The wave attenuation in the model the Darcy flow boundary condition was greater than that under the non-Darcy flow dition. The transmitted waves passing over non-porous structures look more defor from the linear sinusoidal form than those over the porous structures.  Figure 7 compares the wave reflection coefficients with previous experimenta numerical data (Eigenfunction expansion method, EFEM). The present results from FN-NWT are generally in good agreement with the other datasets for both non-po and porous structure conditions. For porous structures, the wave reflection coeffic were similar to the previous study data, except in the long wave region (kh = 0.5 Notably, at the second maximum peak (around kh = 1.3), the wave reflection coeffic applying the non-Darcy (Forchheimer) flow boundary condition (Equation (12)) sho more similar results compared with the experimental and calculation results than th sult of applying the Darcy flow boundary condition (Equation (11)). The reduction o reflected wave amplitude is noticeable, especially at the respective peaks. While Da method underestimates the reflection coefficient. This is especially the case in the ran kh = 1.2-1.4, where there is a fairly good agreement between the non-Darcy-based NWT and the experimental results. This shows that the non-Darcy approach prod  Figure 7 compares the wave reflection coefficients with previous experimental and numerical data (Eigenfunction expansion method, EFEM). The present results from the FN-NWT are generally in good agreement with the other datasets for both non-porous and porous structure conditions. For porous structures, the wave reflection coefficients were similar to the previous study data, except in the long wave region (kh = 0.5-0.7). Notably, at the second maximum peak (around kh = 1.3), the wave reflection coefficients applying the non-Darcy (Forchheimer) flow boundary condition (Equation (12)) showed more similar results compared with the experimental and calculation results than the result of applying the Darcy flow boundary condition (Equation (11)). The reduction of the reflected wave amplitude is noticeable, especially at the respective peaks. While Darcy's method underestimates the reflection coefficient. This is especially the case in the range of kh = 1.2-1.4, where there is a fairly good agreement between the non-Darcy-based FN-NWT and the experimental results. This shows that the non-Darcy approach produces better results than the Darcy-based approach. more similar results compared with the experimental and calculation results than the re-sult of applying the Darcy flow boundary condition (Equation (11)). The reduction of the reflected wave amplitude is noticeable, especially at the respective peaks. While Darcy's method underestimates the reflection coefficient. This is especially the case in the range of kh = 1.2-1.4, where there is a fairly good agreement between the non-Darcy-based FN-NWT and the experimental results. This shows that the non-Darcy approach produces better results than the Darcy-based approach. Table 1 shows the difference between the present results and previous calculation data [31] for the porous structure case (Figure 7b case). In a low frequency region, the results of the non-Darcy flow condition deviate more from the previous calculation results. However, the error of the results under non-Darcy flow condition was greatly reduced as the wave frequency increased.   Table 1 shows the difference between the present results and previous calculation data [31] for the porous structure case (Figure 7b case). In a low frequency region, the results of the non-Darcy flow condition deviate more from the previous calculation results. However, the error of the results under non-Darcy flow condition was greatly reduced as the wave frequency increased.  Figure 8 compares the wave reflection, transmission, and dissipation coefficients over a submerged non-porous and porous dual breakwater. The results for the porous structures were calculated by applying the non-Darcy flow boundary condition. The wave reflection (K R ) and transmission (K T ) coefficients are defined as Equations (25) and (26), and these coefficients define the wave dissipation coefficient (K D ) [20,30,31].

Comparison of Non-Porous and Porous Structures
where H R and H T mean the wave heights of reflected and transmitted waves, respectively. When the wave reflection by the non-porous structure was large, the attenuation of the reflected wave by the porous structure was also large. In the non-porous structure, there are fluctuations in the transmission coefficient (kh = 0.5-1.5). In the porous structure, the transmission coefficient continues to decrease, but the decrease is moderate with shorter wavelengths. Then, in case of the non-Darcy flow, the transmission coefficient becomes minimal at about kh = 2.0 and increases slightly with shorter wavelengths. The dissipation coefficient increases markedly in the kh = 0.7-1.0 region, with a maximum of about kh = 2.0 for the non-Darcy flow condition.
the transmission coefficient continues to decrease, but the decrease is moderate with shorter wavelengths. Then, in case of the non-Darcy flow, the transmission coefficient becomes minimal at about kh = 2.0 and increases slightly with shorter wavelengths. The dissipation coefficient increases markedly in the kh = 0.7-1.0 region, with a maximum of about kh = 2.0 for the non-Darcy flow condition.
The energy dissipation by the non-porous structure was relatively small; however, the energy dissipation by the porous structure was quite large, as intuitively expected. In addition, the porous structure modeling based on the Darcy flow dissipates more energy, which, in turn, results in smaller reflected and transmitted waves. As a result, the wave dissipation coefficient (KD) in the Darcy flow condition was approximately 0.1-0.15 larger than that in the non-Darcy flow condition. Comparisons with previous results show that the non-Darcy model produced better comparisons with the measured data than the Darcy model.  The energy dissipation by the non-porous structure was relatively small; however, the energy dissipation by the porous structure was quite large, as intuitively expected. In addition, the porous structure modeling based on the Darcy flow dissipates more energy, which, in turn, results in smaller reflected and transmitted waves. As a result, the wave dissipation coefficient (K D ) in the Darcy flow condition was approximately 0.1-0.15 larger than that in the non-Darcy flow condition. Comparisons with previous results show that the non-Darcy model produced better comparisons with the measured data than the Darcy model. Figure 9 shows the wave reflection, transmission, and dissipation coefficients according to the distance (L b ) between dual structures under constant wavelength conditions (kh = 1.0). L b is set from 2.1 to 7.6 m. In the non-porous structure, both reflection and transmission coefficients (K R , K T ) tend to repeat peaks and troughs at the interval of L b /λ = 0.5. In case of porous structure, a similar fluctuation is seen with a significant decrease in amplitudes and mean values. Instead, unlike the non-porous case, significant wave energy dissipation can be observed, as shown in K D . It is of interest that the transmission coefficient of the wave passing over the porous dual structure is insensitive to the gap distance. The wave energy dissipation coefficients of the porous case remain around 0.6 over the kh range considered, showing the effectiveness of the applied dual porous breakwater. J. Mar. Sci. Eng. 2023, 11, 1648 12 of 18 Figure 9 shows the wave reflection, transmission, and dissipation coefficients according to the distance (Lb) between dual structures under constant wavelength conditions (kh = 1.0). Lb is set from 2.1 to 7.6 m. In the non-porous structure, both reflection and transmission coefficients (KR, KT) tend to repeat peaks and troughs at the interval of Lb/λ = 0.5. In case of porous structure, a similar fluctuation is seen with a significant decrease in amplitudes and mean values. Instead, unlike the non-porous case, significant wave energy dissipation can be observed, as shown in KD. It is of interest that the transmission coefficient of the wave passing over the porous dual structure is insensitive to the gap distance. The wave energy dissipation coefficients of the porous case remain around 0.6 over the kh range considered, showing the effectiveness of the applied dual porous breakwater.     In Figure 10a, a slight phase difference was observed with different distances between structures at long wave condition in the region of the reflected waves (x = 10-25 m), but the phases of transmitted waves were almost same. The phase shift of reflected waves was less pronounced in shorter waves because seabed structures are less affected by shorter waves. The dimensionless structure spacing normalized by wave length (Lb/λ) was used in Figure 13 for three representative wave lengths. In Figure 10a, a slight phase difference was observed with different distances between structures at long wave condition in the region of the reflected waves (x = 10-25 m), but the phases of transmitted waves were almost same. The phase shift of reflected waves was less pronounced in shorter waves because seabed structures are less affected by shorter waves. The dimensionless structure spacing normalized by wave length (L b /λ) was used in Figure 13 for three representative wave lengths.  Figure 13 compares their differences according to wavelength (kh). As the wavelength becomes shorter (kh increases), the difference in each coefficient according to the distance between the dual porous structure was not large (Figure 13). This is because the seabed effect decreases as the water depth increases relative to the wavelength (h/λ). Regardless of wavelength, the values of Lb/λ corresponding to the maximum and minimum values of each coefficient are similar. However, the effect of distance between structures decreased in larger kh (i.e., in deeper water or with shorter wave length). As shown in Figure 8, the dissipation coefficient increases with shorter wavelengths.  Figure 14 compares the wave characteristics for various incident wave heights. The transmission coefficient increased as the incident wave height increased (Figure 14b), while the dissipation coefficient decreased (Figure 14c). However, its effect on the wave reflection coefficient was not significant.  Figure 13 compares their differences according to wavelength (kh). As the wavelength becomes shorter (kh increases), the difference in each coefficient according to the distance between the dual porous structure was not large (Figure 13). This is because the seabed effect decreases as the water depth increases relative to the wavelength (h/λ). Regardless of wavelength, the values of L b /λ corresponding to the maximum and minimum values of each coefficient are similar. However, the effect of distance between structures decreased in larger kh (i.e., in deeper water or with shorter wave length). As shown in Figure 8, the dissipation coefficient increases with shorter wavelengths. Figure 14 compares the wave characteristics for various incident wave heights. The transmission coefficient increased as the incident wave height increased (Figure 14b), while the dissipation coefficient decreased (Figure 14c). However, its effect on the wave reflection coefficient was not significant.  Figure 14 compares the wave characteristics for various incident wave heights. The transmission coefficient increased as the incident wave height increased (Figure 14b), while the dissipation coefficient decreased (Figure 14c). However, its effect on the wave reflection coefficient was not significant.  Figure 14 also shows the differences in wave coefficients according to linear and fully nonlinear calculations. As the wave non-linearity increases (the incident wave height increases), the differences in the transmission and dissipation coefficients were observed more clearly, while the reflection coefficients remain similar. When comparing the linear and fully nonlinear results, their differences are noticeable in higher incident wave amplitudes at kh values lower than 0.7.  Figure 14 also shows the differences in wave coefficients according to linear and fully nonlinear calculations. As the wave non-linearity increases (the incident wave height increases), the differences in the transmission and dissipation coefficients were observed more clearly, while the reflection coefficients remain similar. When comparing the linear and fully nonlinear results, their differences are noticeable in higher incident wave amplitudes at kh values lower than 0.7.
The snapshots of wave elevations ( Figure 15) are compared under the condition of kh = 1.2, where the difference in the wave reflection coefficient is the maximum according to Figure 14a. Figure 15a shows the relative wave elevation to the incident wave height (η/H I ).
As the wave height increases, the shapes of the reflected waves show slight differences, but the shapes of the transmitted waves are appreciably deformed compared with the linear sinusoidal wave form. More deformation is seen in higher incident waves, as shown in Figure 15a. The snapshots of wave elevations ( Figure 15) are compared under the condition of kh = 1.2, where the difference in the wave reflection coefficient is the maximum according to Figure 14a. Figure 15a shows the relative wave elevation to the incident wave height (η/HI). As the wave height increases, the shapes of the reflected waves show slight differences, but the shapes of the transmitted waves are appreciably deformed compared with the linear sinusoidal wave form. More deformation is seen in higher incident waves, as shown in Figure 15a. The distribution of wave-induced dynamic pressure in the submerged structures is shown in Figure 16. The pressure was made dimensionless by the magnitude of wave dynamic pressure, a value obtained as shown below: The distribution of wave-induced dynamic pressure in the submerged structures is shown in Figure 16. The pressure was made dimensionless by the magnitude of wave dynamic pressure, a value obtained as shown below: Figure 16 shows the pore pressure distribution in the porous structure for large incident wave height (H I = 0.12 m) and kh = 1.2. Figure 16a,b present linear simulation results, and Figure 16c,d show fully nonlinear simulation results. As the wave passes over submerged structures, the deformation of wave shapes and phases by the submerged structures leads to differences in the pore pressure distribution between linear and fully nonlinear cases. The pressure magnitude/fluctuation of the of the front structure is generally higher than that of the rear structure. The pressure differences between the linear and nonlinear cases are also more noticeable at the rear structure due to more wave deformations after passing the front structure in the fully nonlinear results. Based on previous studies, applying Darcy's flow boundary condition was sufficiently valid under small permeability conditions [36]. In the case of relatively high permeability, however, more accurate results can be obtained by applying the non-Darcy flow boundary condition, especially in short wave regions, as illustrated in Figure 6. When the flow in the highly porous structures was assumed to be Darcy flow, the reflected and transmitted waves were greatly attenuated, which implies greater wave dissipation than in the non-Darcy flow case. In general, when the porous structure was used, the fluctuations of reflection and transmission coefficients corresponding to the wavelengths or the distances between dual structures were reduced. The effect of the porous structure was greater for transmitted waves than reflected waves, and more wave energy was dissipated when using the porous structure. However, the transmission coefficients were relatively insensitive to the gap distances.
In addition, by using the developed dual-domain FN-NWT, the nonlinear deformations of transmitted waves by the porous dual structure were clearly observed. As the incident wave height increased, more nonlinear wave deformations were observed due to higher nonlinear effects. The present dual-domain FN-NWT was based on potential flow theory. More realistic physics can be described by using the CFD model in the porous Based on previous studies, applying Darcy's flow boundary condition was sufficiently valid under small permeability conditions [36]. In the case of relatively high permeability, however, more accurate results can be obtained by applying the non-Darcy flow boundary condition, especially in short wave regions, as illustrated in Figure 6. When the flow in the highly porous structures was assumed to be Darcy flow, the reflected and transmitted waves were greatly attenuated, which implies greater wave dissipation than in the non-Darcy flow case. In general, when the porous structure was used, the fluctuations of reflection and transmission coefficients corresponding to the wavelengths or the distances between dual structures were reduced. The effect of the porous structure was greater for transmitted waves than reflected waves, and more wave energy was dissipated when using the porous structure. However, the transmission coefficients were relatively insensitive to the gap distances.
In addition, by using the developed dual-domain FN-NWT, the nonlinear deformations of transmitted waves by the porous dual structure were clearly observed. As the incident wave height increased, more nonlinear wave deformations were observed due to higher nonlinear effects. The present dual-domain FN-NWT was based on potential flow theory. More realistic physics can be described by using the CFD model in the porous region. However, it is difficult to reproduce the details of the complex flows with many rubles and numerous sub-element structures, even with the advanced CFD programs consuming several orders higher computational effort and time. In this regard, the developed dual-domain FN-NWT can practically be applied as a first-cut design tool.

Conclusions
In this study, a fully nonlinear numerical wave tank (FN-NWT) with a porous domain was developed to analyze wave propagation characteristics over submerged porous breakwaters. A three-step calculation process was devised to accurately reflect the interaction between the dual domains. The Darcy and Forchheimer (non-Darcy) flow conditions were used as the fluid and porous interface boundary conditions in the time domain analysis. In the submerged porous structures, the non-Darcy flow model performed better against the corresponding experimental and previous calculation results than the Darcy flow model, as the wavelength became shortened. The wave transmission coefficient of the submerged porous structure was reduced significantly compared with that of the non-porous structure. In particular, applying the Darcy flow condition resulted in greater excessive wave attenuation than non-Darcy flow condition. In the non-Darcy flow condition, the energy dissipation by the submerged porous breakwater increased, with kh reaching a maximum at kh = 2. As the wavelength became shorter (kh increased), the reflection and transmission coefficients decreased with increased energy dissipation regardless of structure spacing. As the incident wave height increased, the transmission coefficient grew relatively larger. The difference between linear and fully nonlinear results was more pronounced at a low frequency (small kh values). The present FN-NWT can also calculate pressure distribution inside the porous breakwater and the resultant forces against it exerted by nonlinear waves.
The present FN-NWT, coupled with the fluid and porous medium domain at each time step and the application of a three-step calculation scheme, can simulate the nonlinear wave and porous structure interaction, even for a relatively large permeability coefficient. Although the FN-NWT has limitations as a two-dimensional potential flow model, it can sufficiently analyze the characteristics of typical long-crested reflected and transmitted waves, thereby saving computational effort and time. The methodology developed in this study can also be easily extended to the three-dimensional model. The FN-NWT model can be applied to various bottom geometries, tsunamis, and irregular waves.

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