Effect of Open Boundary Conditions and Bottom Roughness on Tidal Modeling around the West Coast of Korea

: The aim of this study was to investigate the effect of open boundary conditions and bottom roughness on the tidal elevations around the West Coast of Korea (WCK) using an open-source computational fluids dynamics tool, the TELEMAC model. To obtain a detailed tidal forcing at open boundaries, three well-known assimilated tidal models—the Finite Element Solution (FES2014), the Oregon State University TOPEX/Poseidon Global Inverse Solution Tidal Model (TPXO9.1) and the National Astronomical Observatory of Japan (NAO.99Jb)—have been applied to interpolate the offshore tidal boundary conditions. A number of numerical simulations have been performed for different offshore open boundary conditions, as well as for various uniform and non-uniform bottom roughness coefficients. The numerical results were calibrated against observations to determine the best fit roughness values for different sub-regions within WCK. In order to find out the dependence of the tidal elevation around the WCK on the variations of open boundary forcing, a sensitivity analysis of coastal tide elevation was carried out. Consequently, it showed that the tidal elevation around the WCK was strongly affected by local characteristics, rather than by the offshore open boundary conditions. Eventually, the numerical results can provide better quantitative and qualitative tidal information around the WCK than the data obtained from assimilated tidal models.


Introduction
The Yellow Sea (YS) is located between the mainland of China and the Korean Peninsula, and its depth is about 44 m on average, with a maximum depth of 152 m. The sea bottom slope gradually rises toward the East China coast and more abruptly toward the Korean Peninsula, and the depth gradually increases from north to south. Southern YS is bounded by the East China Sea (ECS) along a line running northeastward from the mouth of the Changjiang River to the southwestern tip of Korea. ECS covers the area originating around the Taiwan Strait at about 25° N, extending northeastward to the Kyushu Island, Japan, and to the Korea Strait, and bounded by the Ryukyu Island chains on the southeastern open ocean (Figure 1). On the broad and shallow shelf under the Yellow and the East China Sea (YECS), the flow is dominated by strong semidiurnal M2 tidal currents with the superimposition of semidiurnal S2, diurnal K1 and O1 currents [1][2][3][4]. Approaching the WCK, the tidal amplitudes vary from 4 to 8 m, and the speed of the tidal current may increase to more than 1.56 m/s near the coasts.
Along with strong tidal currents, the WCK is surrounded by lots of islands including extremely irregular and indented coastlines, the tidal propagation and shoaling processes on the shore in this broad flat region are very complicated and sensitive to the variation of tidal elevation.
The study on the tidal circulation in this region has been carried out by various authors, such as Le Fevre et al. [5], Kim [6], Naimie et al. [7], and Suh [8] using 2D depth-averaged finite element models on unstructured grids, or Choi [2], Tang [9], Kantha et al. [10], Kang et al. [11] and Lee et al. [12] using 3-D finite difference models on structured grids. However, there are still very few authors using data obtained from different assimilated tidal models to set up the offshore boundary conditions. Most authors have applied a constant boundary forcing condition (BFC) or tidal charts, such as Ogura's chart [1,2,13,14]) and Nishida's chart [3,15,16] or applied a very coarse observation data combined with partly constant boundary forcing conditions [11,17]. Only Le Fevre et al. [5] and Suh [8] have applied tidal charts from FES 94.1 and FES 2004 models for OBC. Using FVCOM, we recently applied a tidal chart from NAO.99Jb. In principle, the tidal information around the WCK can be obtained from assimilated global and regional tidal models [18]. In principle, the tidal information around the WCK can be obtained from assimilated global and regional tidal models. The assimilation model, which is a combination of a dynamical model and observed data, was actively advanced in the research group of the modern global tide model, and pioneered by [19]. By the assimilation model, the problem of the empirical method with regard to the resolution can be compensated by a hydrodynamic model, and any inadequate potential of the hydrodynamic model can be recovered by observed data [20]. Recently, using the unprecedented sea surface height data of satellite altimeter TOPEX/POSEIDON (T/P), several research groups have developed highly accurate assimilation models especially in an open ocean ( [5,[20][21][22][23][24]). T/P was launched on 10 August 1992, which measures sea level relative to the ellipsoid with an overall accuracy of approximately 5 cm [25]. Due to the high accuracy of T/P altimeter data for those assimilation models, it can be used to extract relatively accurate tidal boundary conditions around the marginal sea within YECS.
Among the global assimilation models, FES2014 (https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html) and TPXO9.1 (http://volkov.oce.orst.edu/tides/tpxo9_atlas.html) atlases are open to the public and can be freely downloaded. The FES2014 atlas was produced by Legos and CLS Space Oceanography Division, and distributed by Aviso, with support from CNES (https://datastore.cls.fr/catalogues/fes2014-tide-model/). It is distributed as a tidal chart on grid resolution of 1/16° based on the dynamical finite element tide model CEFMO, and data assimilation code CADOR. The finite element grid of the model CEFMO has about 2,900,000 computational nodes and covered 34 tidal constituents. In addition to above two global tide models, the regional tide model NAO.99Jb (https://www.miz.nao.ac.jp/staffs/nao99/index_En.html) was developed by Matsumoto et al. (2000) for the region around Japan with a resolution of 1/12°, covering from 110 to 165° E and from 20-65° N so that YECS was entirely included in this region. This model was nested with the global model NAO.99b to provide open boundary conditions to the regional model and the assimilation was performed with both T/P data and 219 coastal tidal gauges obtained around Japan and Korea. Most T/P-derived altimetric tides can be inaccurate in shallow coastal oceans because the resolution enforced on data analysis by the spatial resolution of T/P is simply inadequate near ocean margins.
Therefore, we first validated the tidal information of tidal constituents 1 , 1 , 2 and 2 provided by three assimilated tidal models; FES2014, NAO99Jb and TPXO9.1 against the observation at 21 tidal gauges along WCK. Figure 2 shows the comparison of amplitude (left) and phase (right) obtained from the assimilated tidal models with observation data obtained from gauging stations along the WCK. It shows that while the phases are better agreement with observation, the amplitudes obtained from the assimilated tidal models still poorly agreed with the observed data. More detailed quantitative comparisons between the amplitude and phase obtained from the assimilated tidal models and observation data are shown in Table A1 (in Appendix A). Wherein it shows NAO99Jb provides a better agreement of the amplitude for the constituents 2 and 2 , while FES2014 provides a little better fit of the amplitude for 1 and 1 . It also shows the constituent M2 has the largest amplitude in comparison with other tidal constituents; about 5-8 times, 2.5-2.8 times, and 6.6-11.8 times larger than K1, S2 and O1 amplitudes, respectively. In addition, the resolutions (1/6°, 1/12° and 1/16°) of three tidal models are not fine enough along the coastlines, so they cannot provide enough detailed quantitative and qualitative tidal information along the WCK. Therefore, in order to obtain more detail on tidal information along the WCK, an accuracy numerical simulation with higher grid resolutions is needed.
This study is to introduce several numerical investigations using higher grid resolutions and applying different open boundary conditions extracted from three well-known assimilated tidal models FES2014, TPXO9.1 and NAO.99Jb then calibrated with different uniform and non-uniform bottom roughness values for different sub-regions within WCK to investigate the effect of open boundary condition and bottom roughness on the tidal elevations from major constituents (M2, K1, S2 and O1) around the West Coast of Korea.
Based on the form number (F), the ratio of the amplitudes of 1 , 1 , 2 , 2 suggested by [26] is shown in Equation (1).
whereby the mixed tide is classified following the value of the form number F; if the value of F less than 0.25 or larger than 1.25, then it is classified as semidiurnal or diurnal, respectively. Otherwise, if the value of F is located between 0.25 and 1.25, it is considered mixed. Figure 3 shows the values of the form number along the WCK are below 0.25 according to the data obtained from observation at 21 tidal gauges, as well as from three assimilated tidal models, FES2014, TPXO9.1 and NAO.99Jb; therefore, the mixed tide prevails semidiurnal in the WCK.

Numerical Modeling
As mentioned above, since the study region was quite shallow, the numerical simulation was based on an open-source software, the TELEMAC-2D (http://www.opentelemac.org) model. It is well-known software with an abundant history of development and application over many years in fluvial and maritime hydraulics [27,28]. The TELEMAC-2D is capable of capturing the complicated bathymetry and coastlines of the WCK surrounded by lots of islands including extremely irregular and indented coastlines by using unstructured grids.

Basic Equations
The TELEMAC-2D has solved the depth-averaged Navier-Stokes equations with two-equation turbulence closure models (k-ε) using the finite element method, as follows: The continuity equation and momentum equations where h is the water depth; U and V are the depth averaged velocity components; The depth averaged kinetic energy k and its dissipation rate ε are: where ′ is the fluctuating velocity and the overbar represents an average over time; i.e., = + ′ ; U i (=U or V) is an average over time of velocity components; and P is the production term:
In this region, particularly near the shallow coastline, the effect of bottom friction becomes prominent, the bottom friction is calculated by: The value of the coefficient Cf is obtained from the calibration procedure, which is shown in Section 3.1 below.
The boundary condition for turbulence kinetic energy and its dissipation are defined as where δ is a distance from the wall, = 0.1 of a local mesh size; * is friction velocity of the wall.

At the Free Surface
We assume the wind force on the free surface is neglected as mentioned above, i.e., = = 0.

At the Open Boundary
Tidal harmonic parameters obtained from each assimilated tidal models FES2014, NAO99Jb and TPXO9.1 were combined into the open boundary conditions, as follows: = ∑ where Hi is the water elevation, is amplitude; T is period and is phase of each tidal constituent.
In this study, four major tidal constituents M2, K1, S2 and O1, were considered. The harmonic constants of tidal constituents for open boundary conditions were obtained from three assimilated tidal models; FES2014, TPXO9.1 and NAO.99Jb.
Turbulence kinetic energy and its dissipation at the open boundary were set as:

Study Region
The study area covered the entire YS and most of the ECS in spherical coordinate. The bathymetry data of 30 arc-seconds resolution was downloaded from the General Bathymetric Chart of the Ocean (GEBCO2014: https://www.ngdc.noaa.gov/mgg/ibcm/ibcmdvc.html) as shown in  Figure 5. The open boundary was expanded to the shelf edge of YECS which was advantageous to use offshore co-oscillating tidal conditions as a boundary forcing in order to minimize disturbances by the geographical configuration and nonlinear tidal response near-shore. Open boundaries adjacent to the offshore also enabled the numerical results around WCK to minimize the influence of unexpected numerical instabilities at open boundaries.
In this study, the numerical simulation applied a hybrid horizontal resolution of an unstructured grid using a free pre-processing Blue Kenue software tool (https://nrc.canada.ca/en/researchdevelopment/products-services/software-applications/blue-kenuetm-software-tool-hydraulicmodellers). Herein the unstructured triangular grid had a horizontal resolution varying from 0.2 to 1 km in the WCK, 3-5 km in the East coast of China, and 9-12 km in the interior of YS and ECS. After running a number of simulations to check for grid independence, the final grid of 203,927 grid cells was used for the simulation in this study. Bathymetry was interpolated to each nodal point of the horizontal grid using Blue Kenue software as well ( Figure 5).
A total of 21 tidal gauges were collected around the WCK to provide the observed tidal elevation data of semidiurnal M2 and S2 and diurnal K1 and O1 constituents in YECS. As shown in Figure 6, the WCK can be divided into three sub-regions based on the common characteristics of their bathymetry and coastlines-the first region was Kyunggi Bay near Incheon (I), the second region was around the estuary of the Geum River near Gunsan (G), and the third region was around the south-

Numerical Results and Discussion
The numerical simulations were carried out on our high-performance computing cluster (http://cfdlabsnu.com) for the real-time simulation of 60 days with time step of 20 s. The numerical results obtained from the last 30 days simulation were used for the following analyses.

Applying Uniform Roughness Coefficient
Previous tidal calculations in the YECS region were mostly using the quadratic bottom friction law, in which the typical value of roughness coefficient was chosen as a uniform value of 0.0025, as shown in [2,16,29,30]. We again have verified the numerical model with various uniform roughness coefficients that ranged from 0.0015 to 0.03 for different offshore open boundary conditions (OBCs) interpolated from three assimilated tidal models (FES2014, TPXO9.1 and NAO.99Jb).
The WCK was divided into three different sub-regions-Mokpo (M), Gunsan (G), and Incheon (I)-as shown in Figure 7, for applying different bottom friction coefficients. The borderline of this sub-regions was set based on a water depth contour of around 50 m. Different values of bottom friction coefficient were applied with regards to the geographic characteristics.

Applying Non-Uniform Roughness Coefficient
As shown above, applying a uniform bottom roughness coefficient for the entire WCK was inappropriate. According to the bathymetry data, the bottom slope was relatively milder in the Incheon region than in Mokpo and Gunsan regions. In addition, there were lots of islands around the Mokpo region rather than the Incheon region. We can expect that more energy dissipation occurred in the Mokpo and Gunsan regions than in the Incheon region. The larger values of the bottom friction coefficient were applied to Mokpo and Gunsan regions than the Incheon region. In addition, the roughness coefficient values for the Gunsan region were set to be the largest because its bed form was more variable than the Mokpo region. Table A5 shows the list of simulation cases applying non-uniform bottom roughness coefficients to different sub-regions Mokpo, Gunsan, Incheon and other regions within WCK. As shown in Tables A6-A9, once the non-uniform bottom roughness coefficients were applied, the mean RMSE values for M2 from the varied ranged of 14-36 cm down to a value of about 11 cm with OBCs obtained from FES2014 and TPXO9.1. These values varied from 16-41 cm down to 12-13 cm with OBCs obtained from NAO99Jb while, there was no significant improvement for K1, S2 and O1. In the WCK, M2 was the dominant constituent, followed by S2 and K1, as mentioned by Teague et al. (1998). In addition, the semidiurnal tides M2 and S2, as well as the diurnal tides K1 and O1, had similar tendencies responding to the local effects and OBCs. Therefore, M2 was presented for the semidiurnal tide, and K1 presented for diurnal tides in most following numerical results. The comparisons of M2 and K1 amplitudes between the observations and simulation results with uniform and non-uniform bottom roughness coefficients are shown in Figures 8 and 9. It clearly shows the simulation results obtained from non-uniform bottom roughness coefficients have significantly improved the results for M2 and K1. The bottom friction coefficient is a physical parameter characteristic of bottom materials and bathymetries. Therefore, it cannot vary following each tidal constituent. From Tables A5-A9, an acceptable selection from all test cases, which can deliver reasonable results for all tidal constituents (M2, K1, S2 and O1), was the case name "FES/NAO/TPX 2-4" whereby the bottom friction coefficient Cf took a value of 0.0025 for the Mokpo region, 0.0035 for the Gunsan region, and 0.002 for Incheon and other regions. These values are in the range suggested by [31,32] applying also for four tidal constituents (M2, K1, S2 and O1).

Response of the Tide around the WCK to the Open Boundary
As mentioned in Section 2.2 above, the modeling region has been bounded by three different open boundaries; A, B and C, as shown in Figure 3. Averaged tidal amplitudes at each boundary A, B and C obtained from three tidal models are shown in Table A10, in which the nodal boundary amplitudes are averaged with weighting by the length of each segment of boundary cells. The tidal amplitudes along the boundary A were significantly larger than the values along the boundaries B and C. Particularly, the tidal amplitude of M2 along the boundary A were larger than the values along the boundaries B and C; about 3.3 times and 9.3 times, respectively. The tidal amplitudes at the boundary C were the smallest because when the tides propagated from the deep region of the East Sea (its depth was about 2000 m) through boundary C via the shallow region at the Korea Strait (its depth was less than 80 m), the tidal amplitudes were significantly damped due to shoaling process on this shallow shelf region. Table A11 shows the difference of the interpolated tidal amplitudes at three open boundaries (A, B and C) obtained from two global tidal models (FES2014 and TPXO9.1) in comparison with those obtained from the regional tidal model (NAO.99Jb). It shows that the amplitudes of M2, K1, S2 and O1 obtained from three tidal models were not substantially different at the open boundaries A, B and C. Nevertheless, it poses a question whether the open boundary condition at A, B or C will play a substantial role on the tide around the WCK. To identify such influence, we have evaluated the response of coastal tide to open boundary forcing in the WCK, using the sensitivity analysis in the following sections.

Sensitivity Analysis
Response of the Tide around the WCK to Individual Boundary Forcing The first numerical test was performed to find out how the co-oscillation of each boundary effected the tidal elevation along the WCK using uniform forcing at open boundaries A, B and C.
Assuming uniform boundary tide constituent for M2 or K1, an amplitude of 100 cm was forced on each boundary A, B, and C, individually. The resultant amplitudes for all gauge locations obtained from each open boundary are shown in Figure 10. It shows that only the tidal wave from the open boundary B produced much higher coastal amplitude than other incident boundary waves from open boundaries A and C. Whereby the tidal amplitudes of M2 and K1 approached 350 and 160 cm in the Kyunggi Bay within area I, respectively; whereas the force propagated from the boundaries A and C was damped on the way to the WCK before it approached the coast. Defant [33] also noticed that the tidal phenomenon in the entire ECS was almost exclusively conditioned by those water-masses which penetrated through the canals between the Ryukyu Islands. Figure 10 (12) where a1 and a2 are different boundary tidal amplitudes in comparison, and � � =1,2 is a response of coastal amplitude at the i-th gauge location.
An increase of tidal amplitudes of constituents M2 and K1 was averaged over gauge locations within a specified region (I, G or M). As shown in Table A12, once the boundary amplitude of M2 constituent increased from 20 to 40 cm, 40 to 60 cm, 60 to 80 cm, and 80 to 100 cm, the mean coastal amplitude increased 3.25, 2.62, 2.01 and 1.91 cm per 1 cm increase in the open boundary amplitude, respectively. It shows that the higher the increment of the boundary amplitude, the lower the mean increase of the coastal amplitude per 1 cm. Meanwhile, the mean increase of the coastal amplitude of K1 was about 1.0 cm (ranged from 0.91 to 1.10), once the boundary amplitude of K1 constituent increased from 20 to 40 cm, 40 to 60 cm, 60 to 80 cm, and 80 to 100 cm. Figure 11 shows the coastal M2 and K1 amplitudes at gauge locations corresponding to forcing with increasing amplitudes of 20, 40, 60, 80 and 100 cm at open boundary B. It shows that as the amplitudes increased at open boundary B, the M2 amplitude significantly increased once it entered the Mokpo region (M), and held this tendency until it reached to the estuary of Geum River, then it dropped until leaving the Gunsan region (G). Thereafter, it started to increase again when entering the Incheon regions (I) reaching the maximum value within Kyunggi Bay, and then significantly decreased. It means that the M2 amplitude was very sensitive within Kyunggi Bay. While K1 was linearly increased a little when it entered the WCK.
In addition, Figure 11 also shows that once the tidal amplitude of M2 at boundary B was about 20 cm, the resultant amplitudes propagated within Kyunggi Bay had little different tendencies in comparison with those obtained from the amplitude larger than 20 cm, since it was not remarkably damped. Table A12 also shows that regions I and G were most sensitive to the boundary amplitude for M2 simulation, once its amplitude was less than 60 cm. Once the amplitude of M2 increased more than 60 cm at open boundary B, it caused less impact on the increase of M2 amplitude along the WCK, since the value of Da did not significantly increase. Meanwhile, an increase of K1 amplitude at the boundary B with an increment of 20 cm caused less significant change on mean Da within all three regions of WCK.
Due to the nonlinearity in the advection terms of the governing equation, and the tidal dissipation by the nonlinear bottom shear stress on the shallow water, the increasing coastal amplitude was found to not be proportional to the increase of offshore tidal amplitude. For a more detailed analysis of the response of WCK to the tidal amplitude at open boundary B, the nonlinear response of tidal amplitude within WCK was evaluated by the following sensitivity estimation: As shown in Figure 12, until the tide approached the tip of the region M, in which , was close to 1, the M2 amplitude increased proportionally to the boundary amplitude. However, when it propagated along the coastline, the increase of amplitude started to depress. As it was pointed out, a substantial depression of amplitudes occurred at some locations in the region M. The sensitivity , dropped mildly during tide propagation through the region G, then it significantly dropped in the region I. Being different from the M2 constituent, the coastal K1 amplitude was less sensitive in the entire WCK after the tip of region M.
From the sensitivity analysis in this section, we can estimate the response of coastal amplitude in the WCK to the difference in open boundary amplitudes of the three established boundary conditions. As shown in Table A10, the averaged amplitude obtained from three different offshore models on the boundary B was in the range from 55.5 to 56.9 cm for M2, and from 21.4 to 21.6 cm for K1. As mentioned in the previous section, the variation of mean coastal amplitude on total gauges was approximated from 1.91 to 3.25 cm for M2, and from 0.91 to 1.10 cm (about 1 cm) for K1 per 1 cm of boundary amplitude variation shown in Table A12, respectively. In Table A11, the RMSE difference in tidal amplitude between the interpolated boundary conditions by two global models (FES2014 and TPXO9.1) and the regional model (NAO99Jb) was found to be 1.

Comparison between the Numerical Results and Observations
As mentioned above, many authors such as [1][2][3] have carried out the tidal simulations in the Yellow Sea and East China Sea (YECS) region using only four dominant tidal constituents; M2, K1, S2 and O1. Particularly, [4] analyzed 13 significant tidal constituents (M2, K1, S2, O1, MM, P1, MU2, N2, MKS2, L2, K2, M4 and MS4) to conclude that M2 was the dominant constituent, followed by S2 and K1 in this region. In this study, we continued to conduct the tidal modeling for this region, taking into account four dominant tidal constituents. Figure 13 shows a comparison of amplitudes between the observations and numerical results of M2 and K1 at gauge stations. It shows a great improvement in comparison to Figure 2, and the numerical results implying different OBCs obtained from three assimilated tidal models (FES2014, TPXO9.1 and NAO.99Jb) were similar.
A question has arisen regarding the difference between the real water level observed around the WCK in comparison with the water level obtained from the numerical results contributed to by only four dominant tidal constituents. Figure 14 below shows a comparison of water level between numerical results and observations at typical gauge stations in three different regions (Incheon, Gunsan and Mokpo). In detail, Table A14 shows a comparison of water levels between the observations and numerical results at 21 gauge stations. Overall, it shows that the water levels obtained from the simulation were overestimated for the neap tide, and underestimated for spring tide. Most stations show the difference in water level was below 20%; at some specific stations located in the Incheon Region, the difference was up to 30%. Nevertheless, as mentioned by the study of Teague et al. [4], after analysis on the contribution of 13 significant tidal constituents on the water level in this region, they noticed that the tides were found to account for at least 85% of the sea surface height variance in the Yellow Sea. Therefore, the water level obtained from the numerical results was reasonable, since its differences compared to the observations can be partly compensated with the amount of 15% mentioned in [4].

Conclusions
The main goal of this study was to use various numerical investigations to figure out the response of coastal tides in the WCK to the open boundary conditions and bottom roughness. After the application of three different open boundary conditions interpolated from three different assimilated tidal models (FES2014, NAO.99Jb and TPXO9.1), it has been shown that there were no significant differences between the responses of tidal amplitudes in the WCK induced by three open boundary conditions obtained from three assimilated tidal models. In addition, the numerical simulation of the tidal flow in the WCK should not use a uniform bottom roughness coefficient. Due to the complicated bathymetry, indented coastlines and bed variability of the WCK, it caused strong local effects on the tides in this region. Therefore, a non-uniform bottom roughness should be applied to the modeling whereby the smallest value can be applied for Incheon, a larger value for Mokpo, and the largest value for Gunsan. The largest value of the bottom roughness coefficient was applied to the Gunsan region because its bed form was more variable than other regions. The numerical results show that the accuracy of the modeling of the tidal elevation around the WCK was strongly dependent on the bottom roughness rather than the offshore tidal boundary conditions. Moreover, the numerical results can provide not only a better fit to the observations but also higher spatial resolutions in comparison to the results obtained from assimilated tidal models around the WCK.
However, it should be noted that the numerical results obtained from this study were still limited due to the coarse resolution (30 arcs/second) of bathymetry obtained from GEBCO2014, which was not sufficient to capture the real geometries, whose sizes were less than such resolution. Therefore, a further study with a higher resolution is necessary in order to obtain a more precise prediction of the tidal current and its elevation around the WCK. Furthermore, the wind forcing on the sea surface and the tidal energy dissipation should be taken into account.

Acknowledgments:
The authors also would like to thank the anonymous reviewers for their valuable and constructive comments to improve our manuscript.

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