2-D Characteristics of Wave Deformation Due to Wave-Current Interactions with Density Currents in an Estuary

In this study, numerical simulations were conducted in order to understand the role of wave-current interactions in wave deformation. The wave-current interaction mechanisms, wave reflection and energy loss due to currents, the effect of incident conditions on wave-current interactions, the advection-diffusion characteristics of saltwater, and the effect of density currents on wave-current interactions were discussed. In addition, the effect of saltwater–freshwater density on wave-current interactions was investigated under a hypopycnal flow field via numerical model testing. Turbulence was stronger under the influence of wave-current interactions than under the influence of waves alone, as wave-current interactions reduced wave energy, which led to decreases in wave height. This phenomenon was more prominent under shorter wave periods and higher current velocities. These results increase our understanding of hydrodynamic phenomena in estuaries in which saltwater–freshwater and wave-current pairs coexist.


Introduction
In the ocean, various continually interacting physical external forces. In particular, estuaries are characterized by freshwater flows from land and saltwater flows from the ocean, which meet and form a major pathway for material transportation. Estuaries are complex, therefore it is difficult to investigate their hydrodynamic characteristics. Thus, it is critical to analyze the dynamic interactions between complex estuary features, including wave-current interaction mechanisms. Furthermore, it is necessary to understand the effect of a density current, which is generated by density differences between saltwater and freshwater, on wave-current interactions.
In the estuary that exhibits highly complex hydraulic characteristics, the circulation type varies significantly according to vertical salinity distribution, tide, and wave [1,2]. However, if two fluids differing in density meet at the estuary, three types of flow and material transport patterns are largely exhibited [3,4]. As shown in Figure 1, homopycnal flow occurs in the case of no density difference, hypopycnal flow occurs when the influx flow has a lower density, and hyperpycnal flow occurs when the influx flow has a higher density. Hypopycnal flow generally occurs in the estuary when fresh river water flows in to ocean and interacts with ocean waves. However, most studies that attempted to analyze the wave-current interaction mechanism at the estuary region did not consider the density discrepancy between freshwater and saltwater, and were performed under homopycnal flow conditions. Meanwhile, due to the developments of measurement equipment, wave transformation, wave set-up, wave set-up, flow rate, salinity, and temperature at the estuary are being observed [5][6][7][8][9]. In some studies, a model of the estuary was composed in a three-dimensional experimental tank to measure wave variation based on wave-current interaction and vertical salinity distribution influenced by density current behavior at the estuary [10,11]. Onsite investigations or transport of threedimensional experiment to the estuary requires significant time and money and only allows for limited information to be recorded at measurement points. Three types of estuarine hydrodynamics based on the density differences between river and sea/lake water (Modified from Bates [3] and Boggs [4]). (a) Homopycnal flow; (b) hypopycnal flow; (c) hyperpycnal flow.
Therefore, numerical simulations are frequently performed to analyze the hydraulic characteristics of the estuary. Most of studies have used the model based on the depth-averaged model [12][13][14] or the 3-D model based on σ-coordinate system [15,16]. This numerical model cannot directly simulate wave behavior due to a large discrepancy between the calculation grid and time. To consider the wave reaction, the wave radiation stress estimated from Simulating WAves Nearshore (SWAN [17]) or WAve Model (WAM [18]) are substituted into the momentum equations [14,19,20].
Various theoretical, experimental and observational studies on waves propagating with currents have previously been conducted in the field of coastal and ocean engineering [21][22][23][24][25]. In particular, there have been many studies on the effects of wave kinematics (changes in the wave number and frequency due to shoaling and refraction) and dynamics (changes in wave steepness and wave-action conservation) according to wave-current interaction [26][27][28][29]. Wave kinematics include the effects of depth and current in wavelength variation and dispersion (Doppler shift) according to opposing or following current. Wave dynamics include energy and action conservation and the variations in wave height and water level [24,30]. Umeyama [31] and Lee et al. [32] explored the vertical distribution of flow and turbulence structures under wave-current interactions in order to analyze wave-current interaction mechanisms. Smith et al. [10] investigated wave deformation by wave-current interactions in a 3-D experimental basin with an estuary model. However, these experimental studies were not able to define the mechanisms of wave deformation and energy loss due to wave-current interactions. Three types of estuarine hydrodynamics based on the density differences between river and sea/lake water (Modified from Bates [3] and Boggs [4]). (a) Homopycnal flow; (b) hypopycnal flow; (c) hyperpycnal flow. Therefore, numerical simulations are frequently performed to analyze the hydraulic characteristics of the estuary. Most of studies have used the model based on the depth-averaged model [12][13][14] or the 3-D model based on σ-coordinate system [15,16]. This numerical model cannot directly simulate wave behavior due to a large discrepancy between the calculation grid and time. To consider the wave reaction, the wave radiation stress estimated from Simulating WAves Nearshore (SWAN [17]) or WAve Model (WAM [18]) are substituted into the momentum equations [14,19,20].
Various theoretical, experimental and observational studies on waves propagating with currents have previously been conducted in the field of coastal and ocean engineering [21][22][23][24][25]. In particular, there have been many studies on the effects of wave kinematics (changes in the wave number and frequency due to shoaling and refraction) and dynamics (changes in wave steepness and wave-action conservation) according to wave-current interaction [26][27][28][29]. Wave kinematics include the effects of depth and current in wavelength variation and dispersion (Doppler shift) according to opposing or following current. Wave dynamics include energy and action conservation and the variations in wave height and water level [24,30]. Umeyama [31] and Lee et al. [32] explored the vertical distribution of flow and turbulence structures under wave-current interactions in order to analyze wave-current interaction mechanisms. Smith et al. [10] investigated wave deformation by wave-current interactions LES-WASS-3D, which is based on the 3-D Navier-Stokes momentum equations, uses the volume of fluid (VOF) method to reproduce complicated free surfaces, including wave breaking, and the porous body model (PBM) to reflect energy dissipation by a permeable structure; the model was based on Hur and Mizutani [52] and Hur et al. [53]. In addition, to analyze sub-grid scale (SGS) turbulence, Smagorinsky turbulence model (STM [54]) for the large eddy simulation (LES) was adopted in the NWT [55]. By considering inertia [56], turbulence [57,58], and laminar flow resistance [58,59] according to porous medium characteristics (diameter, porosity, and shape), the model can directly analyze interactions among waves, structures, and the seabed [49][50][51].
The water density state equation (refer to Appendix A [60]) and kinematic viscosity coefficient (refer to Appendix B [61]) were applied herein to analyze the density current in view of salinity and temperature. Salinity and temperature were quantitatively calculated by simultaneously applying a 3-D advection-diffusion equation. Due to flow separation around structures, present LES turbulence models do not provide accurate calculations of turbulence, which significantly influences the diffusion coefficient applied in the advection-diffusion equation. Thus, a dynamic modeling procedure involving the application of the dynamic eddy viscosity model suggested by Germano et al. [62] and Lilly [63] was applied to re-evaluate Smagorinsky constant (C s ).

Governing Equations
The basic set of calculations consists of the continuity equation (Equation (1)), which includes the source term that generates waves and currents without reflection in 3-D incompressible and viscous fluids, and a modified Navier-Stokes momentum equation (Equation (2)), in which the fluid resistance of the permeable structure is applied.
where x i is the Cartesian coordinate system (x, y, and z); t is time; ν i is the velocity components (u, ν, and w); ρ is the fluid density; p is the fluid pressure; ν is the kinematic viscosity coefficient; ν t is the eddy viscosity coefficient from the turbulence model; γ ν and γ i are the volume and surface porosities, as shown in Figure 2; R i is the fluid resistance term for the porous media; g i is the acceleration of gravity; E i refers to the wave energy damping term (= −βw; where β is the damping factor which equals 0 except for the added damping zones) for the vertical velocity only; Q i refers to the wave and current source terms (= 2ν 3 ∂q * ∂x k ); and q * is the source term required to generate waves and currents, and it is defined as: where q is the flux density; ∆x s is the grid size at the source position (x = x s ). As suggested by Brorsen and Larsen [64], the flux density of the wave and current generation sources q in Equation (3) is gradually increased during three wave periods starting with wave and current generation, as shown in Equation (4).
where V 0 is the horizontal velocity determined by incident condition (in case of wave, based on the 3rd-order Stokes wave theory); and ζ is the intensity factor at the wave source (= (η s + h)/(η 0 + h)); η s is the water surface elevation at the source position, η 0 is the water surface elevation estimated using the 3rd-order Stokes wave theory, and h is the water depth). The constant "2" accounts for two flows propagating to the left and right sides in NWT. The depth integrated quantity of q is adjusted by ζ to achieve the same quantity in the non-reflection condition [65].
As suggested by Brorsen and Larsen [64], the flux density of the wave and current generation sources q in Equation (3) is gradually increased during three wave periods starting with wave and current generation, as shown in Equation (4).
where 0 V is the horizontal velocity determined by incident condition (in case of wave, based on the 3rd-order Stokes wave theory); and ζ is the intensity factor at the wave source ( η is the water surface elevation at the source position, 0 η is the water surface elevation estimated using the 3rd-order Stokes wave theory, and h is the water depth). The constant "2" accounts for two flows propagating to the left and right sides in NWT. The depth integrated quantity of q is adjusted by ζ to achieve the same quantity in the non-reflection condition [65].
F is the volume of fluid (VOF) in each cell [66], which can be expressed as fluid conservation by applying an assumption of fluid incompressibility and a function volume based on PBM, as in Equation (5).
where F is the VOF function. In the numerical simulation, each cell is assigned as one of three types by the VOF function value: fluid cell ( 1 F = ), empty cell ( 0 F = ), and free surface cell ( 0 1 F < < ).

Advection-Diffusion Equations
In density current analyses, quantitative decisions about influential factors are important in order to accurately calculate the fluid density ρ and kinematic viscosity coefficient v , which are substituted into the governing equations. The state equation used to calculate the density and kinematic viscosity coefficient of water is a function of temperature T and salinity S . Therefore, the 3-D advection-diffusion equations adopted herein include Equation (6) for temperature and Equation (7) for salinity.  Figure 2. Definition of the volume and surface porosities. Each cell is flagged as one of three types according to the γ ν value: an empty cell occupied by fluid (γ ν = 1 ), an impermeable cell without fluid flow (γ ν = 0 ), and a permeable cell containing both solid and fluid (0 < γ ν < 1), where V solid and A solid are the volume and surface area of solids in a cell, respectively. (a) F is the volume of fluid (VOF) in each cell [66], which can be expressed as fluid conservation by applying an assumption of fluid incompressibility and a function volume based on PBM, as in Equation (5).
where F is the VOF function. In the numerical simulation, each cell is assigned as one of three types by the VOF function value: fluid cell (F = 1), empty cell (F = 0), and free surface cell (0 < F < 1).

Advection-Diffusion Equations
In density current analyses, quantitative decisions about influential factors are important in order to accurately calculate the fluid density ρ and kinematic viscosity coefficient ν, which are substituted into the governing equations. The state equation used to calculate the density and kinematic viscosity coefficient of water is a function of temperature T and salinity S. Therefore, the 3-D advection-diffusion equations adopted herein include Equation (6) for temperature and Equation (7) for salinity.
where C is temperature T or salinity S, and ε i is the horizontal and vertical diffusion coefficients. In this study, σ c , which is the Prandtl/Schmidt number, has a value of 1.0, as estimated from experimental results [67] and ocean observations [68,69].

Turbulence Model
In general, to reproduce complete turbulence, the calculation area must be larger than the representative scale of the fluid, and the grid must be set smaller than the minimum turbulence scale. It is essential to consider turbulence when modeling a 3-D phenomenon, therefore the grid number of Re 9/4 (Re; Reynolds Number) is required. However, it is practically impossible conduct numerical analysis using this number of grids. Thus, the use of a turbulence model is practical for reproducing turbulence. To this end, the STM [54] was adopted in this study. However, a model constant (C s ) must be applied in the eddy viscosity model for the fluid. For this, the appropriate C s can be calculated through the dynamic modeling procedure suggested by Germano et al. [62] and Lilly [63].
The dynamic modeling method expresses the physical quantity included in the analytical grid as functions of time and space. This method does not require prior designation of unknown model constants or fine tuning and can accurately predict laminar flow and asymptotic behavior near walls. In addition, the model constant has a negative value showing energy backscattering, as the model constant C s is applied temporally and spatially.

Eddy Viscosity Model
In the STM, length and speed are expressed as L = C s ∆ and V = L|S|, respectively. The speed calculation assumes that the pure dissipation equals the energy transferred from the large-scale to small-scale analytical grids, and that the eddy viscosity coefficient ν t formula is expressed as: where S ij is the strain tensor for the grid size, and ∆ is the filter length scale.

Dynamic Eddy Viscosity Model
In the LES grid, the stress τ ij of the sub-grid scale can be modeled via Equation (12) and Equation (13) using STM.
where |S| and ∆ have the same definition in the STM. Likewise, the stress τ ij of the secondary stress can be modeled as Equations (14) and (15) by applying a secondary filter.
A secondary filter (the test filter) with width∆ larger than ∆ is introduced. It is generally assumed that the width of the test filter is two times of that of the grid filter and that the C s of the two filters are the same. These assumptions are adopted in this study. T ij andτ ij are related by the Germano identity.
Substituting the sub-grid scale stress τ ij from Equation (13) and secondary filter stress T ij from Equation (15) into Equation (16) produces relationships such as those seen in Equation (18).
where α is the ratio of the test and grid filters (=∆/∆). There is only one unknown quantity, C s , in Equation (17). Both L ij and M ij are tensors, so Equation (17) is an overdetermined equation system. Applying a least-squares technique to Equation (17) produces a least error value for C s .
A critical problem exists in the determination of the space-temporal function C s . If C s exhibits severe fluctuations or negative values, numerical instability is induced, and it ultimately diverges. To prevent such numerical instability, a space mean in uniform direction of numerator and denominator is required locally for every calculation. In addition, a clipping of C s is also required so that ν + ν t can always be a positive (+) value. On following the above process, most of the limitations of STM can be overcome. Germano et al. [62] proposed 1/2 as the optimal α value in flow calculation while applying the dynamic eddy viscosity model. We also apply α = 1/2 in this study. Figure 3 shows a schematic diagram of a numerical water basin with a variable grid system and damping factors β in an added damping zone for non-reflected wave generation [70]. Each grid cell in the damping zone increases its size by a factor 1.03 with respect to the previous cell as they approach the open boundaries, as shown in Equation (21). This approach serves to make the analysis zone independent from the influence of the numerical water basin boundaries. In addition, wave reflection is controlled by setting the horizontal differences in physical quantities φ i , such as velocity and the VOF function, to zero at the numerical water basin boundary, as shown in Equation (22).  Figure 3 shows a schematic diagram of a numerical water basin with a variable grid system and damping factors β in an added damping zone for non-reflected wave generation [70]. Each grid cell in the damping zone increases its size by a factor 1.03 with respect to the previous cell as they approach the open boundaries, as shown in Equation (21). This approach serves to make the analysis zone independent from the influence of the numerical water basin boundaries. In addition, wave reflection is controlled by setting the horizontal differences in physical quantities i φ , such as velocity and the VOF function, to zero at the numerical water basin boundary, as shown in Equation (22).

Non-Reflected Boundary
For the boundary condition of the impermeable slope in a quadrilateral mesh system according to spatial discretization based on finite difference method, the impermeable condition is adopted in the normal direction and slip condition is adopted in the tangential force for the inclined structure by applying the reasonable boundary [48,51] based on Petit et al. [71].  Figure 4 presents the computed water surface profiles in the numerical water basin during one wave period. Figure 3 shows that the wave spatial envelope is gradually attenuated in both of the added damping zones by the numerical dispersion caused by the coarse grid system and the wave  For the boundary condition of the impermeable slope in a quadrilateral mesh system according to spatial discretization based on finite difference method, the impermeable condition is adopted in the normal direction and slip condition is adopted in the tangential force for the inclined structure by applying the reasonable boundary [48,51] based on Petit et al. [71]. Figure 4 presents the computed water surface profiles in the numerical water basin during one wave period. Figure 3 shows that the wave spatial envelope is gradually attenuated in both of the added damping zones by the numerical dispersion caused by the coarse grid system and the wave damping term included in the momentum equation; the wave spatial envelope vanishes after x/L i = 2. An added damping zone with a length over 2 L i was applied in this study. From numerical result such as those in Figure 3, one can conclude that the present non-reflected system works efficiently to absorb waves in the damping zone and release energy out of the open boundaries.

Wave-Current Interaction
The hydraulic model experiment in Iwasaki and Sato [72] was simulated to verify the characteristics of wave deformation by wave-current interactions in an open channel. Figure 5 shows the numerical water tank simulated by Iwasaki and Sato [72]. Wave and current generation sources (without reflection) were deployed at each end of the analysis domain, and additional damping zones were included at both ends. For verification, a numerical calculation was conducted with an x direction grid size of 1 cm, a z direction grid size of 0.5 cm, and a time increment of 1/1000 s.   Figure 6. In this study, the experimental results are reproduced well by the calculations. Furthermore, the inclination of the wave height distribution is quite similar in the calculations and the experiments for wave deformations with coexisting wave-current fields. This indicates that the ratios of wave height damping due to wave-

Wave-Current Interaction
The hydraulic model experiment in Iwasaki and Sato [72] was simulated to verify the characteristics of wave deformation by wave-current interactions in an open channel. Figure 5 shows the numerical water tank simulated by Iwasaki and Sato [72]. Wave and current generation sources (without reflection) were deployed at each end of the analysis domain, and additional damping zones were included at both ends. For verification, a numerical calculation was conducted with an x direction grid size of 1 cm, a z direction grid size of 0.5 cm, and a time increment of 1/1000 s.

Wave-Current Interaction
The hydraulic model experiment in Iwasaki and Sato [72] was simulated to verify the characteristics of wave deformation by wave-current interactions in an open channel. Figure 5 shows the numerical water tank simulated by Iwasaki and Sato [72]. Wave and current generation sources (without reflection) were deployed at each end of the analysis domain, and additional damping zones were included at both ends. For verification, a numerical calculation was conducted with an x direction grid size of 1 cm, a z direction grid size of 0.5 cm, and a time increment of 1/1000 s.   Figure 6. In this study, the experimental results are reproduced well by the calculations. Furthermore, the inclination of the wave height distribution is quite similar in the calculations and the experiments for wave deformations with coexisting wave-current fields. This indicates that the ratios of wave height damping due to wavecurrent interaction are similar. These results were used to determine the validity and effectiveness of the numerical model for wave transformation during an interaction with an opposing current.   Figure 6. In this study, the experimental results are reproduced well by the calculations. Furthermore, the inclination of the wave height distribution is quite similar in the calculations and the experiments for wave deformations with coexisting wave-current fields. This indicates that the ratios of wave height damping due to wave-current interaction are similar. These results were used to determine the validity and effectiveness of the numerical model for wave transformation during an interaction with an opposing current.  Figure 7 shows a numerical water tank based on the hydraulic model experiment in Huppert and Simpson [73]. The water tank was 960 cm long and 27 cm wide, and the water depth was 10 cm. A compartment of 30 cm long (x0) was installed in the water tank and filled with the density of saltwater (ρs) of 1.0115 g/cm 3 . The remaining area was filled with the density freshwater (ρf) of 0.9999 g/cm 3 . The density difference (Δρ) was 0.0116 g/cm 3 and the reduced gravity g' (= Δρg/ρf) was 11.4 cm/s 2 . The grid size for verification was Δx = 1 cm, Δy = 1 cm, and Δz = 0.25 cm, and the numerical calculation was carried out for Δt = 1/1000 s.   Figure 7 shows a numerical water tank based on the hydraulic model experiment in Huppert and Simpson [73]. The water tank was 960 cm long and 27 cm wide, and the water depth was 10 cm. A compartment of 30 cm long (x 0 ) was installed in the water tank and filled with the density of saltwater (ρ s ) of 1.0115 g/cm 3 . The remaining area was filled with the density freshwater (ρ f ) of 0.9999 g/cm 3 . The density difference (∆ρ) was 0.0116 g/cm 3 and the reduced gravity g' (= ∆ρg/ρ f ) was 11.4 cm/s 2 . The grid size for verification was ∆x = 1 cm, ∆y = 1 cm, and ∆z = 0.25 cm, and the numerical calculation was carried out for ∆t = 1/1000 s.  Figure 7 shows a numerical water tank based on the hydraulic model experiment in Huppert and Simpson [73]. The water tank was 960 cm long and 27 cm wide, and the water depth was 10 cm. A compartment of 30 cm long (x0) was installed in the water tank and filled with the density of saltwater (ρs) of 1.0115 g/cm 3 . The remaining area was filled with the density freshwater (ρf) of 0.9999 g/cm 3 . The density difference (Δρ) was 0.0116 g/cm 3 and the reduced gravity g' (= Δρg/ρf) was 11.4 cm/s 2 . The grid size for verification was Δx = 1 cm, Δy = 1 cm, and Δz = 0.25 cm, and the numerical calculation was carried out for Δt = 1/1000 s.

Description of Numerical Water Basin and Incident Conditions
A numerical wave tank was defined as shown in Figure 9 in order to numerically investigate the mechanisms of wave-current interactions according to the velocity of the incident wave and current in an open channel. To generate waves based on 3rd-order Stokes waves and currents without reflection, wave and current sources were installed at each end of the analytical domain; damping zones were added at both ends to absorb wave energy. The length of the analytical domain was 5 L i of the incident wavelength, the width was 30 cm, and the water depth was 30 cm. In addition, the length of the damping zones was set to more than double the incident wavelength to absorb wave energy and minimize reflection.
Wave-current interaction numerical simulations were performed for 63 different cases involving three incident wave heights (H i = 3, 5, and 7 cm), three incident wave periods (T i = 1.2, 1.5, and 1.8 s), and seven current velocities (V c = 0, 5, 10, 15, 20, 25, and 30 cm/s). The wave celerity (C i ) is based on 3 rd -order Stokes wave theory. The calculation grids were ∆x = 2 cm, ∆y = 2 cm, and ∆z = 1 cm. In all simulations, the wave was generated after current generation in order to stabilize the flow field.
Water 2020, 12, x FOR PEER REVIEW 12 of 30

Description of Numerical Water Basin and Incident Conditions
A numerical wave tank was defined as shown in Figure 9 in order to numerically investigate the mechanisms of wave-current interactions according to the velocity of the incident wave and current in an open channel. To generate waves based on 3rd-order Stokes waves and currents without reflection, wave and current sources were installed at each end of the analytical domain; damping zones were added at both ends to absorb wave energy. The length of the analytical domain was 5 Li of the incident wavelength, the width was 30 cm, and the water depth was 30 cm. In addition, the length of the damping zones was set to more than double the incident wavelength to absorb wave energy and minimize reflection.
Wave-current interaction numerical simulations were performed for 63 different cases involving three incident wave heights (Hi = 3, 5, and 7 cm), three incident wave periods (Ti = 1.2, 1.5, and 1.8 s), and seven current velocities (Vc = 0, 5, 10, 15, 20, 25, and 30 cm/s). The wave celerity (Ci) is based on 3 rd -order Stokes wave theory. The calculation grids were Δx = 2 cm, Δy = 2 cm, and Δz = 1 cm. In all simulations, the wave was generated after current generation in order to stabilize the flow field.  Figure 10 shows the wave height spatial distributions for various current velocity conditions when the incident wave (Ti) period was 1.2 s. As shown in Equation (23), the wave height was calculated using surface elevations (Equation (24)) from three periods and described via the nondimensional wave height (H/Hi) after the flow and wave field stabilized.

Distribution of Wave Heights
where 2 k = and   Figure 10 shows the wave height spatial distributions for various current velocity conditions when the incident wave (T i ) period was 1.2 s. As shown in Equation (23), the wave height was calculated using surface elevations (Equation (24)) from three periods and described via the non-dimensional wave height (H/H i ) after the flow and wave field stabilized.

Distribution of Wave Heights
where k = 2 and k = k max − 1 represent the bottom and top of the NWT, and ∆z k is the grid size in the z direction. In Figure 10, for all the cases, the wave height is higher than the incident wave height (Hi) at x/Li = 0. The wave height gradually decreases with wave propagation. This phenomenon is more obvious as the ratio (Vc/Ci) between current velocity (Vc) and incident wave celerity (Ci) increases. A partial standing wave field is generally observed due to wave reflection by the wave-current interactions. These results are in agreement with the experimental results [32,39,74]. In addition, an intersecting point with a wave height of H/Hi = 1 is observed in (a-c); this point moves farther from x/Li = 0 as Hi increases because Ci increases, thus decreasing Vc/Ci. In Figure 10, for all the cases, the wave height is higher than the incident wave height (H i ) at x/L i = 0. The wave height gradually decreases with wave propagation. This phenomenon is more obvious as the ratio (V c /C i ) between current velocity (V c ) and incident wave celerity (C i ) increases. A partial standing wave field is generally observed due to wave reflection by the wave-current interactions. These results are in agreement with the experimental results [32,39,74]. In addition, an intersecting point with a wave height of H/H i = 1 is observed in (a-c); this point moves farther from x/L i = 0 as H i increases because C i increases, thus decreasing V c /C i . Furthermore, as V c /C i increases, the partial standing wave wave-height distribution is marginally formed. In this case, the wave height is also relatively high at x/L = 0 in comparison to the point at which the partial standing wave forms. This phenomenon was also reported by Lee et al. [39], who found H i = 5 cm, T i = 1.0 s, and V c = 60 cm/s. In other words, shorter wave periods, lower wave heights, and higher current velocities lead to lower wave orbital motion velocity. Therefore, this phenomenon appears to be a result of the predominance of wave energy loss over wave reflection. Nevertheless, further study is required to verify this mechanism. Figure 11 shows the spatial distribution of the H/H i according changes in the V c under T i = 1.5 s. This figure shows a tendency similar to that in Figure 10, which has T i = 1.2 s. However, the V c /C i is small compared to the case with T i = 1.2 s; thus, the effects caused by a decrease in wave height due to decreases in interactions decrease, and changes in the wave height rate become weaker. For this reason, points (H/H i = 1) generally move far from x/L i = 0. The wave height distribution becomes narrower because the range over which V c /C i fluctuates is relatively narrow. Furthermore, as Vc/Ci increases, the partial standing wave wave-height distribution is marginally formed. In this case, the wave height is also relatively high at x/L = 0 in comparison to the point at which the partial standing wave forms. This phenomenon was also reported by Lee et al. [39], who found Hi = 5 cm, Ti = 1.0 s, and Vc = 60 cm/s. In other words, shorter wave periods, lower wave heights, and higher current velocities lead to lower wave orbital motion velocity. Therefore, this phenomenon appears to be a result of the predominance of wave energy loss over wave reflection. Nevertheless, further study is required to verify this mechanism. Figure 11 shows the spatial distribution of the H/Hi according changes in the Vc under Ti = 1.5 s. This figure shows a tendency similar to that in Figure 10, which has Ti = 1.2 s. However, the Vc/Ci is small compared to the case with Ti = 1.2 s; thus, the effects caused by a decrease in wave height due to decreases in interactions decrease, and changes in the wave height rate become weaker. For this reason, points (H/Hi = 1) generally move far from x/Li = 0. The wave height distribution becomes narrower because the range over which Vc/Ci fluctuates is relatively narrow.   Figure 12 shows the wave height spatial distribution with changes in the Vc under Ti = 1.8 s. This figure also shows a tendency similar to those in Figures 8 and 9, and also shows prolonged wavecurrent interaction period characteristics. As the Ti increases, the range over which Vc/Ci fluctuates decreases, resulting in a narrower wave height distribution. In addition, as the incident wave period and wave height increase, the Ci increases and the decreases in wave height weaken, resulting in gentler inclinations in the wave height distribution.  Figure 12 shows the wave height spatial distribution with changes in the V c under T i = 1.8 s. This figure also shows a tendency similar to those in Figures 8 and 9, and also shows prolonged wave-current interaction period characteristics. As the T i increases, the range over which V c /C i fluctuates decreases, resulting in a narrower wave height distribution. In addition, as the incident wave period and wave height increase, the C i increases and the decreases in wave height weaken, resulting in gentler inclinations in the wave height distribution.  Figure 12 shows the wave height spatial distribution with changes in the Vc under Ti = 1.8 s. This figure also shows a tendency similar to those in Figures 8 and 9, and also shows prolonged wavecurrent interaction period characteristics. As the Ti increases, the range over which Vc/Ci fluctuates decreases, resulting in a narrower wave height distribution. In addition, as the incident wave period and wave height increase, the Ci increases and the decreases in wave height weaken, resulting in gentler inclinations in the wave height distribution. Considering the above discussion of Figures 8-10, the wave orbital motion velocity and current velocity have the same direction under the wave trough and different directions under the wave crest within the wave-current interaction. As confirmed in the existing experiment, more intense wavecurrent interactions occur in the wave crests than in the wave troughs. This leads to wave reflection and wave height decay.

Characteristics of Turbulent Kinitic Energy (TKE)
The mean phase-averaged TKE ( K ) is the average turbulence over ten wave periods, as shown in Equation (25). This TKE is calculated using Equations (26) and (27) where u′ , v′ , and w′ are turbulence components; t v is the eddy viscosity coefficient estimated using STM; s C is a model constant calculated by the dynamic eddy viscosity model; and Δ is the filter length scale. Figure 13 shows the vertical distribution of the K according to the ratio of current velocity and wave celerity (Vc/Ci). This figure shows that the value of the mean phase-averaged TKE is larger with wave-current interactions than with waves alone. The intensity generally increases with higher Vc/Ci values. In particular, still water with higher velocity due to wave orbital motion shows higher TKE. Therefore, wave height decreases with increased TKE. The relationship between TKE and wave height decay in the coexisting wave-current field will be comprehensively discussed later, along with wave energy characteristics. Considering the above discussion of Figures 8-10, the wave orbital motion velocity and current velocity have the same direction under the wave trough and different directions under the wave crest within the wave-current interaction. As confirmed in the existing experiment, more intense wave-current interactions occur in the wave crests than in the wave troughs. This leads to wave reflection and wave height decay.

Characteristics of Turbulent Kinitic Energy (TKE)
The mean phase-averaged TKE (K) is the average turbulence over ten wave periods, as shown in Equation (25). This TKE is calculated using Equations (26) and (27), as suggested by Christensen [75]. Equation (26) can be used to calculate the TKE of a grid scale (K GS ) using the flow velocity calculated by the numerical model, and Equation (27) can be used to calculate TKE of a sub-grid scale (K SGS ) by applying the values calculated by LES. The center of the analysis zone is the measuring point.
where u ,ν , and w are turbulence components; ν t is the eddy viscosity coefficient estimated using STM; C s is a model constant calculated by the dynamic eddy viscosity model; and ∆ is the filter length scale. Figure 13 shows the vertical distribution of the K according to the ratio of current velocity and wave celerity (V c /C i ). This figure shows that the value of the mean phase-averaged TKE is larger with wave-current interactions than with waves alone. The intensity generally increases with higher V c /C i values. In particular, still water with higher velocity due to wave orbital motion shows higher TKE. Therefore, wave height decreases with increased TKE. The relationship between TKE and wave height decay in the coexisting wave-current field will be comprehensively discussed later, along with wave energy characteristics.  The mean phase-averaged TKE increases with increasing Hi, Hi/Li, and Ursell number (Ur). This increment is larger near still-water levels. In addition, with longer incident periods, the mean phaseaveraged TKE increases in sections, except for those near the still-water level. Under the influence of wave-current interactions, TKE is closely related with velocity due to the orbital motion of waves and current velocity. That is, as the incident wave height and period increase, the velocity of water particles increases, and thus the TKE is increased by wave-current interactions.
Nevertheless, it is difficult to directly correlate TKE and the aforementioned wave height decay because they are not proportional. The wave energy conditions differ, therefore wave energy decay and wave height decay are not proportional to TKE. This relationship is closely related to the wave energy described later and will be discussed in more detail in the following section. Wave height decay is caused primarily by turbulence in the coexisting wave-current field, therefore the mechanism of wave attenuation will be analyzed by exploring changes in TKE, wave energy, and wave height decay.
(a)  Figure 14 shows the vertical distribution of the K according to the incident wave period (T i ). The mean phase-averaged TKE increases with increasing H i , H i /L i , and Ursell number (Ur). This increment is larger near still-water levels. In addition, with longer incident periods, the mean phase-averaged TKE increases in sections, except for those near the still-water level. Under the influence of wave-current interactions, TKE is closely related with velocity due to the orbital motion of waves and current velocity. That is, as the incident wave height and period increase, the velocity of water particles increases, and thus the TKE is increased by wave-current interactions.
Nevertheless, it is difficult to directly correlate TKE and the aforementioned wave height decay because they are not proportional. The wave energy conditions differ, therefore wave energy decay and wave height decay are not proportional to TKE. This relationship is closely related to the wave energy described later and will be discussed in more detail in the following section. Wave height decay is caused primarily by turbulence in the coexisting wave-current field, therefore the mechanism of wave attenuation will be analyzed by exploring changes in TKE, wave energy, and wave height decay.  Figure 14 shows the vertical distribution of the K according to the incident wave period (Ti).
The mean phase-averaged TKE increases with increasing Hi, Hi/Li, and Ursell number (Ur). This increment is larger near still-water levels. In addition, with longer incident periods, the mean phaseaveraged TKE increases in sections, except for those near the still-water level. Under the influence of wave-current interactions, TKE is closely related with velocity due to the orbital motion of waves and current velocity. That is, as the incident wave height and period increase, the velocity of water particles increases, and thus the TKE is increased by wave-current interactions.
Nevertheless, it is difficult to directly correlate TKE and the aforementioned wave height decay because they are not proportional. The wave energy conditions differ, therefore wave energy decay and wave height decay are not proportional to TKE. This relationship is closely related to the wave energy described later and will be discussed in more detail in the following section. Wave height decay is caused primarily by turbulence in the coexisting wave-current field, therefore the mechanism of wave attenuation will be analyzed by exploring changes in TKE, wave energy, and wave height decay. (a)

Wave Energy Loss Characteristics
Wave energy ( E ) lost due to wave-current interactions is calculated by applying the simulated data in Equation (28) after a stable wave-current coexisting field is attained. In Equation (29), the mean energy loss ( L E ) is calculated for three periods from the energy ratio of the prior and following waves. Subsequently, the average L E of five wavelengths is used to calculate the wavelengthaveraged energy loss ( L E ).  Figure 15 shows the L E according to the Vc/Ci between current velocity and wave celerity under various wave conditions. Wave energy loss increases as Vc/Ci increases. This phenomenon can be explained with the fact that TKE increases as Vc/Ci increases, which results in increased energy loss. The energy decay rate was expected to increase because the TKE becomes larger as the incident wave height increases. The energy decay due to turbulence increases because higher wave energies accompany higher wave heights, but the general decay rate is relatively low compared to that for low wave heights.

Wave Energy Loss Characteristics
Wave energy (E) lost due to wave-current interactions is calculated by applying the simulated data in Equation (28) after a stable wave-current coexisting field is attained. In Equation (29), the mean energy loss (E L ) is calculated for three periods from the energy ratio of the prior and following waves. Subsequently, the average E L of five wavelengths is used to calculate the wavelength-averaged energy loss (E L ). Figure 15 shows the E L according to the V c /C i between current velocity and wave celerity under various wave conditions. Wave energy loss increases as V c /C i increases. This phenomenon can be explained with the fact that TKE increases as V c /C i increases, which results in increased energy loss. The energy decay rate was expected to increase because the TKE becomes larger as the incident wave height increases. The energy decay due to turbulence increases because higher wave energies accompany higher wave heights, but the general decay rate is relatively low compared to that for low wave heights.
The aforementioned analytical results for wave height distribution, TKE, and wave energy show that the TKE increases with increasing V c /C i rate. This leads to wave energy and height decay. This numerical mechanism can be used to understand wave height decay due to wave-current interactions. The aforementioned analytical results for wave height distribution, TKE, and wave energy show that the TKE increases with increasing Vc/Ci rate. This leads to wave energy and height decay. This numerical mechanism can be used to understand wave height decay due to wave-current interactions.

Wave Reflection Characteristics
Partial standing waves due to wave-current interactions are also found in the aforementioned wave height distribution. Equation (30), which was proposed by Healy [76], was used to calculate the reflection coefficient with the wave-current interactions. In addition, the values calculated for five partial standing wave field wavelengths were averaged to produce the average reflection coefficient.

Wave Reflection Characteristics
Partial standing waves due to wave-current interactions are also found in the aforementioned wave height distribution. Equation (30), which was proposed by Healy [76], was used to calculate the reflection coefficient with the wave-current interactions. In addition, the values calculated for five partial standing wave field wavelengths were averaged to produce the average reflection coefficient.
where H max and H min are the wave heights at the node and the anti-node in a partial standing wave field, respectively. Figure 16 shows the wavelength-averaged reflection coefficient (K R ) according to the ratio (V c /C i ) between current velocity and wave celerity under various wave conditions. The wave reflection rate generally increases as V c /C i increases. Wave reflection due to current is generated by a strong interaction under the wave crest, where the wave orbital motion velocity and current velocity oppose each other. Therefore, the reflection rate increases as the current velocity increases. Meanwhile, (a-c), which feature the same incident wave height, show lower reflection rates with shorter incident wave periods. The wavelength is not long in waves with shorter periods, therefore the exposure time under the wave crest decreases. Therefore, the wave-current interaction time becomes shorter under an opposite flow velocity, resulting in a lower reflection rate. In addition, the current by the incident wave height does not appear to have strong effects on the wave reflection rate, which is similar to the reflective wave characteristics.
Wave reflection by currents, which was reported by an existing theoretical study to be insignificant [37], is therefore confirmed, and the mechanism is elucidated above. The reflection rate is up to 0.068 under the incident conditions applied herein. Such a wave reflection rate cannot be ignored when analyzing hydraulic phenomena in coexisting wave-current fields.
Water 2020, 12, x FOR PEER REVIEW 20 of 30 reflection rate generally increases as Vc/Ci increases. Wave reflection due to current is generated by a strong interaction under the wave crest, where the wave orbital motion velocity and current velocity oppose each other. Therefore, the reflection rate increases as the current velocity increases. Meanwhile, (a-c), which feature the same incident wave height, show lower reflection rates with shorter incident wave periods. The wavelength is not long in waves with shorter periods, therefore the exposure time under the wave crest decreases. Therefore, the wave-current interaction time becomes shorter under an opposite flow velocity, resulting in a lower reflection rate. In addition, the current by the incident wave height does not appear to have strong effects on the wave reflection rate, which is similar to the reflective wave characteristics. Wave reflection by currents, which was reported by an existing theoretical study to be insignificant [37], is therefore confirmed, and the mechanism is elucidated above. The reflection rate is up to 0.068 under the incident conditions applied herein. Such a wave reflection rate cannot be ignored when analyzing hydraulic phenomena in coexisting wave-current fields. 3.1.6. Mean Wave Level Characteristics Figure 17 shows the spatial distribution of the mean water level according to the ratio between the current velocity and wave celerity, Vc/Ci. As shown in Equation (31), the mean water level was calculated using surface elevations (Equation (24)) from three periods and described via the nondimensional mean level ( / i H η ) after the wave and flow field stabilized. 3.1.6. Mean Wave Level Characteristics Figure 17 shows the spatial distribution of the mean water level according to the ratio between the current velocity and wave celerity, V c /C i . As shown in Equation (31), the mean water level was calculated using surface elevations (Equation (24)) from three periods and described via the non-dimensional mean level (η/H i ) after the wave and flow field stabilized.
Water 2020, 12, 183 20 of 30 In Figure 17, The mean water level also increases when V c /C i increases because the wave energy decreases due to wave-current interactions, which show a tendency to increase because the wave energy decreases further when V c /C i increases. Meanwhile, these types of wave level distributions create currents due to wave level differences; currents, in turn, increase the initial current velocity, which will accelerate wave height dissipation. To analyze this further, it is necessary to study the interactions between the flow field and the wave field.
In Figure 17, The mean water level also increases when Vc/Ci increases because the wave energy decreases due to wave-current interactions, which show a tendency to increase because the wave energy decreases further when Vc/Ci increases. Meanwhile, these types of wave level distributions create currents due to wave level differences; currents, in turn, increase the initial current velocity, which will accelerate wave height dissipation. To analyze this further, it is necessary to study the interactions between the flow field and the wave field.

Wave-Current Interactions with Density Current
The effect of current due to density differences between freshwater and saltwater must be considered, as estuaries contain both freshwater and saltwater. In the past, most studies on wavecurrent interactions did not consider density currents. Therefore, wave-current interactions under the influence of density differences were numerically simulated using a newly suggested numerical model in order to calculate density currents.

Description of Numerical Water Basin and Incident Conditions
A numerical water basin was prepared as shown in Figure 18 to evaluate the effect of density current on wave-current interactions in an estuary. To generate waves without reflection, wave and current sources were installed at both ends of the simulation domain, and damping zones were added to absorb the wave energy. The analysis zone featured a length of 5 Li of the incident wavelength, a width of 50 cm, and a water level of 30 cm; a 1:5 slope was installed to a water depth of 10 cm to serve as the current source.
The simulations were performed under fixed incident wave and current conditions and variable water density. The details are described in Table 1. Two major cases were considered: one with density differences under wave-current interactions, and another without density differences (ρw = ρc). The first major case was divided into three sub-experiments, which involved: high density water mass and low density current flows (ρw > ρc); low density water mass and high density current flows (ρw < ρc); and water mass with density stratification and middle density current flows (ρw1 < ρc < ρw2).
Here, the reduced gravity (g') is that the effective change in the acceleration of gravity acting on one fluid in contact with a fluid of different density due to buoyancy forces. For these, the calculation grids were Δx = 2 cm, Δy = 2 cm, and Δz = 1 cm. In all simulations, current was generated first, and waves were generated only after the flow field stabilized.

Wave-Current Interactions with Density Current
The effect of current due to density differences between freshwater and saltwater must be considered, as estuaries contain both freshwater and saltwater. In the past, most studies on wave-current interactions did not consider density currents. Therefore, wave-current interactions under the influence of density differences were numerically simulated using a newly suggested numerical model in order to calculate density currents.

Description of Numerical Water Basin and Incident Conditions
A numerical water basin was prepared as shown in Figure 18 to evaluate the effect of density current on wave-current interactions in an estuary. To generate waves without reflection, wave and current sources were installed at both ends of the simulation domain, and damping zones were added to absorb the wave energy. The analysis zone featured a length of 5 L i of the incident wavelength, a width of 50 cm, and a water level of 30 cm; a 1:5 slope was installed to a water depth of 10 cm to serve as the current source.
The simulations were performed under fixed incident wave and current conditions and variable water density. The details are described in Table 1. Two major cases were considered: one with density differences under wave-current interactions, and another without density differences (ρ w = ρ c ). The first major case was divided into three sub-experiments, which involved: high density water mass and low density current flows (ρ w > ρ c ); low density water mass and high density current flows (ρ w < ρ c ); and water mass with density stratification and middle density current flows (ρ w1 < ρ c < ρ w2 ). Here, the reduced gravity (g') is that the effective change in the acceleration of gravity acting on one fluid in contact with a fluid of different density due to buoyancy forces. For these, the calculation grids were ∆x = 2 cm, ∆y = 2 cm, and ∆z = 1 cm. In all simulations, current was generated first, and waves were generated only after the flow field stabilized.   Figure 19 shows the flow fields for CASEs 1-4 in Table 1, in which current was generated without considering waves. In Figure 19, four different current shapes can be identified due to density differences, including (a) a uniform-depth current due to a lack of density difference, (b) an undercurrent due to a high-density current, (c) an upper current due to a low-density current, and   Figure 19 shows the flow fields for CASEs 1-4 in Table 1, in which current was generated without considering waves. In Figure 19, four different current shapes can be identified due to density differences, including (a) a uniform-depth current due to a lack of density difference, (b) an undercurrent due to a high-density current, (c) an upper current due to a low-density current, and (d) a middle current due to a middle-density current flowing into the density stratification. Based on these results, various currents may arise from various density difference structures. Therefore, the consideration of differences in density will affect the wave-current interactions in co-existing fields.
Water 2020, 12, x FOR PEER REVIEW 23 of 30 (d) a middle current due to a middle-density current flowing into the density stratification. Based on these results, various currents may arise from various density difference structures. Therefore, the consideration of differences in density will affect the wave-current interactions in co-existing fields.
(a)  Figure 20 shows the wave height variations generated by wave-current interactions with density differences, including (a) interactions between the wave and four different current shapes due to density differences, (b) an interaction between the wave and undercurrent (ρw < ρc), and (c) an interaction between the wave and upper current (ρw > ρc).

Wave Height Variation
All cases in Figure 20a show decreases in wave height due to wave energy loss caused by wavecurrent interactions. In particular, CASE 19 (ρw > ρc) shows considerable wave height decay because the inflowing current is concentrated near the water surface, where the velocity caused by wave orbital motion is large. Therefore, there may be stronger interactions near the water surface, which induce wave height decay due to wave energy loss. In addition, the wave height decay seen in CASE 12 (ρw < ρc) and wave height seen in CASE 20 (ρw1 < ρc < ρw2) are lower than those seen in CASE 5 (ρw = ρc) because the current flows in the lower-middle layers, which feature lower velocity caused by orbital motion, lead to less interaction between waves and currents, which results in less wave height decay than that in CASE 5 (ρw = ρc), in which differences in density are not considered. Figure 19. Four types of density currents with varying density differences. (a) Depth-uniform current (ρ w = ρ c ); (b) Undercurrent (ρ w < ρ c ); (c) Upper current (ρ w > ρ c ); (d) Middle current (ρ w1 < ρ c < ρ w2 ). Figure 20 shows the wave height variations generated by wave-current interactions with density differences, including (a) interactions between the wave and four different current shapes due to density differences, (b) an interaction between the wave and undercurrent (ρ w < ρ c ), and (c) an interaction between the wave and upper current (ρ w > ρ c ).

Wave Height Variation
All cases in Figure 20a show decreases in wave height due to wave energy loss caused by wave-current interactions. In particular, CASE 19 (ρ w > ρ c ) shows considerable wave height decay because the inflowing current is concentrated near the water surface, where the velocity caused by wave orbital motion is large. Therefore, there may be stronger interactions near the water surface, which induce wave height decay due to wave energy loss. In addition, the wave height decay seen in CASE 12 (ρ w < ρ c ) and wave height seen in CASE 20 (ρ w1 < ρ c < ρ w2 ) are lower than those seen in CASE 5 (ρ w = ρ c ) because the current flows in the lower-middle layers, which feature lower velocity caused by orbital motion, lead to less interaction between waves and currents, which results in less wave height decay than that in CASE 5 (ρ w = ρ c ), in which differences in density are not considered. In Figure 20b, CASEs 6, 8, and 11 show the same phenomenon and feature less wave height dissipation than in CASE 5, which has a uniform-depth current. Furthermore, CASEs 6, 8, and 11 show more wave height dissipation under smaller differences in salinity. This is because current generation is more similar to that under a uniform-depth current when the difference in density is small. In Figure 20c, there is more wave height decay in CASEs 13, 15, and 18 than in CASE 5. The cases with upper current also show a tendency toward higher wave height decay as density differences increase, although the effect is not very strong.
These results indicate that vertically varying currents affect wave deformation due to wavecurrent interactions. Therefore, density currents must be considered in order to effectively analyze wave-current interactions in areas with density differences. Figure 21 shows the vertical distribution of the mean phase-average TKE ( K ) under various wave-current interactions due to density differences, including (a) interactions between the wave and four different current shapes due to density differences, (b) an interaction between the wave and an undercurrent (ρw < ρc), and (c) an interaction between the wave and an upper current (ρw > ρc). In Figure 20b, CASEs 6, 8, and 11 show the same phenomenon and feature less wave height dissipation than in CASE 5, which has a uniform-depth current. Furthermore, CASEs 6, 8, and 11 show more wave height dissipation under smaller differences in salinity. This is because current generation is more similar to that under a uniform-depth current when the difference in density is small. In Figure 20c, there is more wave height decay in CASEs 13, 15, and 18 than in CASE 5. The cases with upper current also show a tendency toward higher wave height decay as density differences increase, although the effect is not very strong.

Vertical Distribution of TKE
These results indicate that vertically varying currents affect wave deformation due to wave-current interactions. Therefore, density currents must be considered in order to effectively analyze wave-current interactions in areas with density differences. Figure 21 shows the vertical distribution of the mean phase-average TKE (K) under various wave-current interactions due to density differences, including (a) interactions between the wave and four different current shapes due to density differences, (b) an interaction between the wave and an undercurrent (ρ w < ρ c ), and (c) an interaction between the wave and an upper current (ρ w > ρ c ). In Figure 21a, CASE 5 shows more balanced K under the interactions between waves and a uniform-depth current than do the other cases. CASEs 12,19, and 20 show strong K generation at each current point generated by density differences. In particular, K due to the wave-upper current interaction appears to be strong near the free surface in CASE 19, which features the highest wave height decay (Figure 18). In addition, CASE 5, which features a uniform-depth current, possesses higher K than CASE 12, which features undercurrent and less wave height dissipation, and CASE 20, which features a middle current. This turbulence appears to be the major factor in the wave energy loss and height dissipation and can explain the formation mechanisms of the abovementioned wave height distributions. In Figure 21b, CASEs 8 and 11 feature stronger K when larger salinity differences exist in the lower layer. However, CASE 6, which features the lowest salinity difference, shows a K vertical distribution similar to that in CASE 5; this reduces the wave energy, resulting in the wave height distribution described above. In Figure 21c, CASEs 13, 15, and 18 possess stronger K under larger salinity differences near the free surface. Under smaller salinity differences, K vertical distributions similar to that in CASE 5, which features a uniform-depth current, are observed, which reduces wave energy, resulting in the wave height distribution discussed earlier.

Vertical Distribution of TKE
Overall, wave-current interactions in cases with density differences appear different from those found in cases that do not consider density differences. Therefore, consideration of currents due to density differences is strongly recommended during analysis of wave-current interactions in water bodies with density differences, such as estuaries. In Figure 21a, CASE 5 shows more balanced K under the interactions between waves and a uniform-depth current than do the other cases. CASEs 12,19, and 20 show strong K generation at each current point generated by density differences. In particular, K due to the wave-upper current interaction appears to be strong near the free surface in CASE 19, which features the highest wave height decay (Figure 18). In addition, CASE 5, which features a uniform-depth current, possesses higher K than CASE 12, which features undercurrent and less wave height dissipation, and CASE 20, which features a middle current. This turbulence appears to be the major factor in the wave energy loss and height dissipation and can explain the formation mechanisms of the abovementioned wave height distributions. In Figure 21b, CASEs 8 and 11 feature stronger K when larger salinity differences exist in the lower layer. However, CASE 6, which features the lowest salinity difference, shows a K vertical distribution similar to that in CASE 5; this reduces the wave energy, resulting in the wave height distribution described above. In Figure 21c, CASEs 13, 15, and 18 possess stronger K under larger salinity differences near the free surface. Under smaller salinity differences, K vertical distributions similar to that in CASE 5, which features a uniform-depth current, are observed, which reduces wave energy, resulting in the wave height distribution discussed earlier.
Overall, wave-current interactions in cases with density differences appear different from those found in cases that do not consider density differences. Therefore, consideration of currents due to density differences is strongly recommended during analysis of wave-current interactions in water bodies with density differences, such as estuaries.

Conclusions
In the estuary, most of the laboratory experiments and numerical simulations for analysis of the wave-current interaction were performed under homopycnal flow conditions that do not consider the density discrepancy between the river and ocean. However, hypopycnal flow, which refers to river water spreading out to the upper layer of ocean, actually occurs in estuary. In this study, to simulate wave-current interaction under hypopycnal flow conditions, we have revised the NWT to reproduce the density current based on the existing numerical model. To confirm the effectiveness and validity of the improved numerical model, the model was compared with the results of hydraulic experiments [72,73]. These comparisons confirmed the validity and effectiveness of the newly suggested numerical model in wave-current interaction simulations in estuaries that contain density differences.
This study addressed wave attenuation, wave reflection, and wave energy reduction under wave-current interactions, which have not been addressed in previous studies. In addition, this study numerically explored the effects of current due to density differences on wave-current interactions in an open channel. Wave height decay and wave reflection mechanisms were analyzed in the coexisting wave-current field in terms of TKE and wave energy. The hydraulic characteristics of the numerically simulated wave-current interactions can be summarized as follows: • When the ratio (V c /C i ) between current velocity and wave celerity increases, wave height decay increases. Wave reflection by currents forms a partial standing wave field. • When the V c /C i velocity ratio increases, the overall mean phase-average TKE increases. At higher incident wave heights, wave steepness, and Ursell numbers larger increases occur near still-water levels. Longer incident periods lead to larger increases, except near the still-water level.

•
As a result of turbulence caused by wave-current interactions, wave energy increases as the V c /C i increases. When the incident wave height is high, TKE increases and the wave energy decay rate decreases. Higher wave heights have higher wave energy, therefore the energy decay rate is relatively low.

•
The rate of wave reflection due to currents becomes larger as the V c /C i increases. With longer wave periods, the exposure time in wave troughs, in which the wave orbital motion velocity and current velocity are opposite in direction, is longer, and thus the wave reflection rate is increased. As a result, reflection rates of up to 0.068 were found under the incident conditions applied in this study.

•
The mean water level increases more due to wave-current interactions as the velocity ratio V c /C i increases, which further reduces wave energy.
In conclusion, the flow velocity turbulence element is increased in the coexisting wave-current field by wave-current interactions, which decreases the wave energy. Consequently, the wave height decays and the water level increases. Wave reflection by currents occurs under the wave trough. Thus, the results provide an explanation wave height decay and wave reflection mechanisms due to wave-current interactions.
The influence of currents caused by density differences on wave-current interaction in an estuary can be described as follows: • Various currents are generated depending on the density difference between two fluids. • Wave height decay due to various current-wave interactions is highest in cases involving upper currents, in which strong interactions occur near still-water levels. Cases involving middle current and undercurrent feature less wave height decay than do cases with uniform-depth currents, which feature no density current.

•
The mean phase-average TKE vertical distribution is broad near each type of current. Greater wave height decay values are accompanied by higher mean phase-average TKE.
In general, simulated wave-current interactions caused by density differences show different tendencies than cases that do not account for density differences. Therefore, currents arising from density differences must be considered when analyzing wave-current interactions in water bodies with density differences, such as estuaries.

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

Appendix A. Fluid Density Equations
To analyze density current, it is important to estimate the density of a fluid accurately; thus, Equation (A1) was applied to estimate density according to temperature and salinity, as suggested by Gill [60]. Here, ρ 0 is the density of 4 • C freshwater. ∆ρ T is the increase in density with temperature change and is expressed as Equation (A2). ∆ρ S expresses the change in density with changing salinity and is expressed as Equation (A3). The empirical constants used in the density calculation are shown in Table A1.

Appendix B. Kinematic Viscosity Coefficient Equations
The kinematic viscosity coefficient ν for the fluid is calculated using Equation (A4). Here, the calculated value of Equation (A1) is substituted for the density ρ, and the viscosity coefficient µ is calculated using Equation (A5), which considers the water temperature and salinity, as suggested by Riley and Skirrow [61]. ν = µ/ρ (A4) where, µ 0 is the viscosity coefficient for 4 • C freshwater. ∆µ T indicates the change in the viscosity coefficient with a change in temperature and is expressed via Equation (A6). ∆µ S indicates the change in the viscosity coefficient with a change in salinity and is expressed via Equation (A7). The empirical constants used for the water viscosity coefficient are expressed in Table A2. µ 0 = 1.802863 × 10 −2 g/cm·s b 2 = 1.31419 × 10 −5 b 4 = 2.15123 × 10 −5 b 1 = 6.108600 × 10 −4 b 3 = 1.35576 × 10 −7 b 5 = 3.59406 × 10 −10