The semi-implicit cross-scale hydroscience integrated system model (SCHISM by Zhang et al., 2016) was adopted in this study to simulate an abnormal tide that occurred in Hiroshima Bay and the SIS. This model is a derivative product based on the original SELFE model (v3.1dc) (Zhang and Baptista, 2008). The Navier-Stokes equations are solved in hydrostatic form using semi-implicit finite-element/volume methods, which reduce Courant-Friedrichs-Lewy (CFL) stability constraints and increase numerical efficiency. A higher-order scheme for momentum advection with Explicit Locally Adaptive Dissipation (ELAD) filter is applied to the SCHISM to control excess mass with an iterative smoother. Total Variation Diminishing with two limiter functions (TVD²) is a higher-order implicit advection transport scheme using two space and time limiters. This allows the model to conserve baroclinic instability without filtering out, successfully capturing smaller eddies. These characteristics allow seamlessly stable simulations across varying scales from the creek, lake, river, estuary, shelf to the ocean (Zhang et al., 2023). The model domain in this study included a river (approximately 3 m deep) and the Nankai Trough (approximately 4,900 m deep) with narrow straits and complicated geometry. Therefore, the SCHISM will be an appropriate numerical model for the SIS.

The fundamental equations of the SCHISM model consist of three equations: (Equation 1) the continuity equation, (Equation 2) the momentum equation, and (Equation 3) the transport equation.

where ( x ,   y ) are the horizontal coordinates, z is the vertical coordinate, t is the time, u ( x ,   y ,   z ,   t ) is the horizontal velocity, w is the vertical velocity (upward), η ( x ,   y ,   t ) is the free surface elevation, g is the gravitational acceleration, and C is the concentration of the tracer, such as temperature and salinity. ν is the vertical eddy viscosity, κ is the vertical eddy diffusivity for the tracer, and F is the additional momentum forcings such as the Coriolis effect and air pressure. F h is the horizontal diffusion, and Q is the source/sink mass of the tracers.

To enhance momentum stabilization, (Equation 4) a 5-point Shapiro filter (Shapiro, 1970) was adopted in this study. There are different tools, such as the Laplacian viscosity and biharmonic viscosity, which act similarly to the Shapiro filter. However, the filter is suitable for noneddying regime applications, such as nearshore areas, estuaries, and rivers, because of the following equation:

where u 0 is the velocity at side ‘0’, which is defined between two elements (here, triangular), and u n (n = 1, 2, 3, and 4) is the velocity defined at the other two sides of two triangles (2×2 = 4). γ is a nondimensional Shapiro filter parameter used to define the strength of stability. It is recommended that this parameter be set between 0 and 0.5, and different values of the parameter can be assigned to each node depending on the model instability. To investigate mass exchanges through straits and similar narrow passages, when side velocities are converted to a node velocity in the Shapiro filter, the 5-point Shapiro filter is used again to redistribute excess mass. Therefore, excess mass is redistributed and conserved between sub-basins through straits.

The unstructured mesh covering the Seto Inland Sea and adjacent seas, each with a length of 500 km × 420 km, consisted of 210,505 nodes and 390,461 cells (triangular mesh in Figure 2). The horizontal resolution of the mesh varied from approximately 30 m (near the Ota River mouth and Hiroshima city) to 9,800 m (along the open boundary) depending on the geometry of the coastlines, target areas, and gradients of bathymetry. The 30-m resolution was applied only to the areas shallower than 10 m, while other narrow straits had an approximately 100-m resolution. Considering that the average depth of the SIS is around 38 m, we aimed to avoid situations where the model mesh becomes excessively fine compared to the depths. As a result, the unstructured mesh of the SCHISM model can resolve the complicated coastline geometry of the SIS.

The simulation period was set to 61 days from 10 August to 9 October 2011, which included an abnormal tide event in Hiroshima Bay and two typhoon events. The model was calculated with an interval of 80 seconds. The model output has a 20-minute interval. All the descriptions regarding time used in this paper are based on local time and Japan Standard Time (JST).

Spatially varying filtering parameters were applied using the 5-point Shapiro filter for the sake of numerical stability based on the maximum current velocity and horizontal resolution at each node, with higher values applied to narrow straits and river mouth areas where numerical instability can occur due to strong currents. As a result, 0.2 of the parameters for the normal areas and a maximum of 0.5 for the unstable areas were given. This approach allowed us to increase the model resolution near narrow straits with strong currents frequently developed while the model stability was maintained. The detailed descriptions of the model setup are documented in Jeong and Lee (2023).

Thirty vertical layers in σ coordinate were adopted in the SIS model and LCS² grid (Localized Sigma Coordinates with Shaved Cells) was applied in the model to reduce the instability and pressure gradient error at areas with narrow straits and shallow or steep depths. In addition, the current σ-layer adopted in this study didn’t show significant errors during the simulation, as validated well with field measurements (Jeong and Lee, 2023).

The bottom roughness length ( z 0 ) used near the bottom boundary layer was applied nonuniformly in this study. A higher value of the parameter was defined as the depth ( D ) increase (Equation 5), following the approach of Cowles et al. (2008).

SIS modeling considers external forcings, such as tides, sea surface height (SSH), water temperature, currents, surface winds, air pressure, and river discharge, which can affect sea level changes directly or indirectly. Table 1 describes all the datasets used in the model. The Japanese Meteorological Agency (JMA) Meso-Scale Model (MSM) offers high-resolution data (approximately 5 km) near Japan. Therefore, the MSM was adopted for the hourly surface winds, air pressure, and air temperature, while ERA5 was utilized for hourly precipitation and radiation dew point temperature after interpolation to a grid of the MSM. The specific humidity was calculated from the dew point temperature and air pressure fields of ERA5. Detailed descriptions regarding how each dataset was treated are documented in Jeong and Lee (2023).

The MRI.COM-JPN model from the JMA was utilized primarily for the initial field and boundary forcings in the ocean component of our study. Notably, the model employs a structured grid with an approximate resolution of 2 kilometers, as illustrated in Figure 2. This MRI.COM-JPN was developed by the JMA to enhance coastal monitoring and forecasting systems, covering a wide area from the SIS to the entire coastal seas of Japan (Sakamoto et al. (2019). In contrast, the SIS model, implemented using SCHISM, focuses on the SIS and its immediate vicinity, providing much higher resolution, as demonstrated in Figure 2.

The mean absolute error (MAE) and index of agreement (IOA) (Willmott et al., 2012) were calculated to evaluate the model results. The MAE is a dimensioned measure of average model performance error expressed in units of the variable of interest (Willmott, 2005). The MAE is described as in Equation 6:

where P i is the i_th model result,   O i is the i_th observation, and n is the number of data points.

The IOA adopted in this study is a refined version of the original Willmott’s dimensionless IOA (Willmott, 1981). The refined IOA has a lower limit of −1 and an upper limit of +1, which is double the range of the original IOA (from 0 to +1). The refined IOA is calculated via Equation 7:

The upper equation in Equation 7 is used when the IOA is positive. Basically, an IOA value closer to +1 indicates better performance of a numerical model. When IOA = +0.5, the sum of the absolute errors is half of the sum of the deviation magnitudes of the perfect model and observation. When IOA = 0, Equation 7 is equivalent. If the IOA becomes negative (−0.5, for example), the sum of the absolute errors is double the sum of the deviation magnitude of the perfect model and observation.

The SCHISM model utilizes an unstructured triangular mesh system, making it challenging to depict a straight line for cross-sectional analysis. Therefore, model results were interpolated at new grid points, located along the cross section, as illustrated in Figure 3A. Net fluxes (m³/s) were then calculated by multiplying the velocity normal to the cross section (m/s) and the area (m²) at each grid point, as described in Equation 8. The area at each grid point was obtained by calculating the area of the trapezoid, as depicted in Figure 3B, and expressed in Equation 9. The cross-sectional net flux at a specific time was obtained by summing net fluxes calculated at all mesh and at all vertical layers, as presented in Equation 8.

Here, N is the number of grid points at each cross section, z is a vertical layer (maximum 30^th layer), U n is the velocity normal to the cross section, and H is the total depth.

The number of grid points was determined based on the length of each cross section to achieve consistent resolutions for net flux calculations. For instance, in the Bungo Channel (See Table 2), the length of the cross-sectional line is 29.12 km. It was divided into 101 grid points, making a final resolution of approximately 300 m. Similarly, for other channels, the number of grid points was determined in proportion to their respective lengths, resulting in a cross-sectional resolution of approximately 300 m.

The magnitude of the tides changed as they propagated inside the SIS. The maximum tidal range at Kure station (Figure 5B) was approximately 3.9 m with semidiurnal characteristics (Jeong and Lee, 2023), whereas Hosojima (Shirahama) station in the outer SIS presented low amplitudes of approximately 2 m (1.9 m), as shown in Figures 5D, E. However, Osaka station, located in the eastern part of the SIS, had a lower amplitude of 1.5 m, and diurnal tides became dominant during neap tide. This spatial distribution of tides from the observation and the model matched Yanagi et al. (1982).

Elevation data from tide stations can capture not only tidal signals but also other subtidal components. The subtidal is used to express the signals which are obtained after excluding dominant tidal signals (Payandeh et al., 2019; Laurel-Castillo and Valle-Levinson, 2020) such as storm surges, abnormal tides, and non-tidal long waves. It does not necessarily be shorter than semi-diurnal tide in this study. The long waves mainly indicate CTWs such as edge waves, Kelvin waves, and continental shelf waves including the planetary edge waves. Figure 5A shows the wind vectors (magenta) and air pressure (blue) from the JMA MSM model at the center of the model domain (133.65° E, 33.24° N), indicating the representative condition in the SIS. Two typhoons, TY1 and TY2, passed through the SIS during the simulation period. As TY1 approaches the SIS, low air pressure (973 hPa) and strong winds disturb seawater, causing sea level variations, as shown in Figures 5C, E. However, because tidal components are more dominant than these variations are, it is difficult to distinguish subtidal components from tides.

Daily averaged elevations were calculated to investigate subtidal variations, as shown in Figure 6. The observed subtidal elevations at all stations in Figure 6 exhibit a peak when an abnormal tide occurred. This indicates that abnormal tides occurred not only in Hiroshima Bay but also throughout the whole SIS, and they can be observed at the outer stations of the SIS. Furthermore, at the Osaka, Yoshima, and Shirahama stations, located in the relatively eastern part of the SIS, one more peak was observed due to the storm surge when Typhoon 1 approached the SIS.

For comparison, to remove the influence of tides in the subtidal data analysis, the T-Tide and T-Predic software were applied and the result at Kure station is shown in Figure A1. The daily-averaging method and the T-Tide application results did not show significant difference in subtidal components.

The results from the SIS model could not reproduce these peaks at TY1 and abnormal tide. To investigate why the model failed to describe this phenomenon, sea surface height (SSH) from MRI.COM-JPN utilized as initial and open boundary conditions was compared with the observed and SIS model results in this study. In the comparison, we noticed that MRI.COM-JPN (blue lines in Figure 6) could not reproduce subtidal components well. The difference in SSH in MRI.COM still remained at the location of the open boundary of the SIS model (green lines in Figure 6). Therefore, additional simulations were conducted after correcting the SSH bias at the open boundary. Moreover, interestingly, the subtidal elevation at obc1 had the highest correlation coefficient (r = 0.77) with the SIS model results, although obc2 (r = 0.70) and obc3 (r = 0.57) were closer to the Kure, Oita, and Hosojima stations. This difference might be related to CTWs, as stated by Usui et al. (2021).

To improve the SIS model results and reduce the SSH open boundary conditions, elevations defined at the open boundary were corrected using differences between the observations and the results of the SIS model (ad hoc correction). The Hosojima and Shirahama tidal stations, which are the closest stations to the western and eastern sides of the open boundary, respectively, were chosen for bias correction. As tidal waves propagate further, the waves can be transformed and amplified in the SIS. Therefore, the stations closest to the boundary were adopted to correct the SSH at the boundary.

Figure 7 shows the water level differences between the observation and model results at the Hosojima and Shirahama tidal stations, and these differences are identical to the corrected water levels at the boundary. When correcting the boundary conditions, spatially varying modifications depending on the distance from two stations are considered. In detail, there were 106 nodes at the open boundary, and for the most western (eastern) nodes, we added the difference in subtidal elevation calculated at the Hosojima (Shirahama) station. Then, we interpolated the differences for the other open boundary nodes according to the distances to both ends (the most western and eastern nodes). The ten green lines in Figure 7 represent samples of the interpolated gaps at ten open boundary nodes among the 106 nodes.

Another simulation (blue solid lines in Figure 8) used the same model configuration as the previous simulation and was conducted after the correction was applied. The correction process effectively improved the model results, achieving a lower MAE and higher IOA, as shown in Table 3. The MAEs at subtidal elevations at all stations were markedly reduced. Before correction, the average MAE at the six stations was approximately 9.19 cm. However, after the correction, the value decreased to 3.15 cm, representing a reduction of 66%. Especially at the Osaka and Shirahama stations, the MAE was reduced by approximately 7.45 cm, which accounts for 77% of the original MAE at these stations. These stations were located on the eastern side of the model domain, close to the eastern open boundary, where we utilized the subtidal component for the correction.

The IOA increased by 0.55 (from 0.08 to 0.63) at Hosojima station, closest to the western side of the open boundary, and increased by 0.42 (from 0.44 to 0.86) at the Shirahama station, closest to the eastern side, after the boundary correction. As the subtidal elevation propagated into the SIS, the average IOA in the SIS also improved by 0.50 (from 0.21 to 0.71). In particular, significant improvements were observed at the Kure station, which is located close to the Itsukushima Shrine, where flooding occurred. After correction, the IOA at this station reached 0.69. Regarding the surface elevation before extracting the subtidal component, the IOA did not change significantly because tides were dominant compared to the subtidal component (see the difference in the scales between Figures 5, 6). With respect to the interpretation of the IOA values, it has only statistical meanings for goodness-of-fit of the model result rather than providing threshold values for acceptable model performance in this study. Thus, the interpretation of the IOA values rests with the readers’ judgment.

From an eventwise perspective, the correction effectively improved the calculation of abnormal tides compared to that in the original case at all stations. For Typhoon 1, the nonexistent peak that was previously calculated at the Kure, Oita, and Hosojima stations (red solid lines in Figures 8A, B, E, respectively) was corrected and disappeared after the SSH correction was applied. This erroneous peak originated from MRI.COM-JPN data is shown in Figure 6. Conversely, at the Osaka, Yoshima, and Shirahama stations on the relatively eastern side of the SIS, a peak was observed during Typhoon 1, which was overestimated in the original case (and was also overestimated in the MRI.COM-JPN data). However, after the correction, the overestimation decreased, and the model results aligned better with the observations at the Osaka and Shirahama stations. Unfortunately, at Yoshima station, located at the center of the SIS, the model presented a lower peak than both the observation and the original model (Figure 8D). The corrected sea level at the open boundary was lower than that in other periods, especially at Shirahama station (smaller than 3.9 cm, magenta line in Figure 7) during the period of Typhoon 2. However, an elevation of approximately -13.6 cm was corrected at the most western node during Typhoon 2 (black line in Figure 7). Despite these corrections, only 6.4 cm decreased at the Oita and Hosojima stations close to the western node, which is less than half of the corrected elevation.

Net seawater fluxes crossing the Bungo Channel, Kii Channel, Kanmon Strait, and Harima Nada (Figure 9) were analyzed to investigate the water sources that induced sea level variations and changes in the sources after correction. The net fluxes were averaged daily to remove tidal signals. The locations of the cross sections are presented in Figure 1. Since the SIS is a semienclosed sea with only three passages, the Bungo Channel, the Kii Channel, and the Kanmon Strait, it is possible to calculate the net influx of seawater into the SIS (black lines in Figure 9). We also took discharged fresh water from 23 out of 24 rivers into account in the net flux budget calculation. One river located outside the cross section (not shown here) was not considered in the calculation. The solid (dotted) lines represent the net fluxes from the corrected (original) case. Positive (negative) values indicate inflow into the SIS (outflow from the SIS), respectively. Positive (negative) values at Harima Nada indicate eastward (westward) flow in the SIS.

There were some similarities among the net fluxes at the cross sections in Figure 9. First, the most noticeable pattern of the total fluctuation was that of spring/neap tides. The outflow (negative) of the total was calculated for every neap tide, and it became positive as the tidal range increased. The fortnightly pattern was mainly affected by the net flux in the Bungo Channel (blue lines in Figure 9). During the second and third neap tides among the four neap tides within the simulation period, the net flux at Bungo presented sharp drops, resulting in a strong and long outflow of the net value. Except for the TY1 period, the net flux in the Kii Channel was almost negative. This means that if there is no extreme event, seawater will always be discharged along the Kii Channel, and the main inflow of seawater into the SIS will be determined by the amount of inflow at the Bungo.

Although the inflow into the SIS was mainly determined by the net flux in the Bungo Channel, water circulation could not be effectively induced at the center of the SIS (Harima Nada). The net flux at Harima Nada (green lines) did not follow the trend of the main inflow (blue lines) but rather the reverse trend of the main outflow (red lines), as depicted in Figure 9. This means that the seawater at the center flowed eastward when the seawater discharged from the Kii Channel. In other words, although there was an inflow of seawater to the Bungo Channel, the SIS would not flow eastward without discharge through the Kii Channel. As a result, inflow through the Bungo Channel tends to increase the sea level around the western side of the SIS because the inflowed sea water is not easily transferred to the eastern side without outflow through the Kii Channel.

Two countercurrents (westward) were calculated at Harima Nada during both the TY1 and TY2 events, which can be explained by the wind direction (Figure 5A) and elevation distribution across the SIS (Figure 4). The average elevation along a cross-section (black line in Figure 1) gradually decreased as it extended to the Kii Channel. This dominant gradient induced a continuous eastward net flux at Harima Nada. However, when the westward wind occurred due to typhoons (Figure 5A), spatially different patterns were calculated, and the sharp reverse gradients, which are opposite to the direction of the average gradient, appeared at a distance of 187 km (red line and green line in Figure 4). Additionally, at a distance of 270 km, a similar reverse gradient was calculated for Typhoon 1. Both locations (187 km and 270 km distances) are entrances to the Harima Nada Strait, which has very limited widths. Because of the narrow straits, each basin suffered from a storm surge. Therefore, the sea level increased near all the eastern sides of the straits (187+α and 270+α), and the sea level decreased near all the western sides of the straits (187-α and 270-α) due to the effect of the storm surge. This process induced temporary strong countercurrents at Harima Nada. However, this unstable gradient of sea level could be maintained only by the effects of the typhoons during the study period. Just after the typhoons diminished, the restoring force to recover the original gradient of the averaged elevation became dominant. The force induced strong eastward net fluxes at Hariman Nada just after the passages of TY1 and TY2, as calculated at Nada (green lines) in Figure 9.

Two additional factors contribute to the total net flux: the Kanmon Strait and river discharge. The net flux through the Kanmon Strait (magenta lines in Figure 9) was mainly outward. The reason for the constant outflow can be found in the inflow at the Bungo Channel. Because the influx in the Bungo Channel mainly induced a sea level rise around the western side of the SIS, in the Kanmon Strait, which is located on the western side, a negative outward flux is dominant for balancing the net flux. Regarding river discharge (yellow line in Figure 9), two river discharge peaks were calculated according to the timing of the typhoon: a maximum of 26,660 m³/s during TY1 and 25,557 m³/s during TY2. The river discharges were comparable to the net flux through the Kanmon Strait, implying that many rivers along the coastlines of the SIS can affect sea level variations at a limited scale in the SIS to a certain extent.

MAEs between the original model results and the corrected results for each net flux were calculated to quantify how much each net flux changed after the correction (Table 4). The largest change in the net flux after the correction was calculated at the Bungo Channel, with an MAE of 5742.2 m³/s. A strong and positive effect of the correction on the net flux in the channel can be found during TY1 and TY2. On 1 September, before TY1 approached, the corrected model calculated less inflow than the original model in the Bungo Channel. After three days, less outflow was calculated in the corrected model. These weakened flows at the Bungo Channel mitigated a nonexistent peak at the Kure and Oita stations (as described in Figures 8A, B), resulting in a smaller MAE and higher IOA of subtidal components, as shown in Table 3. For TY2, there was an insignificant difference between the original and corrected values on 21 September when TY2 approached. However, from 22 September to 30 September, the net flux of the corrected results exceeded that of the original one, indicating a greater inflow of the corrected model through the Bungo Channel. This additional inflow helped address the abnormal tide, as shown in Figure 8.

Other net fluxes, such as those in Kii, Kanmon, and Harima, presented slighter changes than those in the Bungo Channel, as shown in Table 4. However, these changes also helped improve the model results. For the Kii Channel, close to Osaka Station, a greater outflow was calculated from 25 August to 2 September, resulting in a slow increase in subtidal elevation at Osaka Station (Figure 8C). The reason why the subtidal component increased even under the outflow at the Kii Channel was because of an inflow from Harima Nada. The amount of inflow from Harima Nada to Osaka Bay exceeded the amount of outflow from Osaka Bay to the outside of the SIS (compare the green lines and red lines in Figure 9). Then, from 3 September to 5 September, the larger inflow and less outflow at the Kii Channel were calculated after the correction. This change induced a steep rise in sea level at Osaka station, as depicted in Figure 8C. From 21 September to 30 September, after the TY2 passed, less outflow was calculated from the Kii Channel after correction. This change caused the sea level to increase at the Osaka station (Figure 8C). Although there were modifications in the net fluxes at the Bungo and Kii Channels, the net flux at the Harima Nada did not significantly change (green lines in Figure 9).

Discharge from the Kanmon Strait acted as a buffer against sea level rises. While the net flux in the Kanmon Strait presented negative values (outflow) almost always, some differences were calculated after the correction. During the TY1 season, a relatively stronger outflow was calculated in the original season. This was because the sea level on the western side of the SIS was greater in the original model than in the corrected model at that time (Figures 8A, B). If there was no discharge from the Kanmon Strait, the sea level in the SIS would have increased.

Similarly, the opposite situation occurred during TY2. A large discharge was calculated in the corrected model, resulting from higher sea levels in the corrected model. This means that although the width of the Kanmon Strait is much narrower than that of other channels, it also plays an important role in mitigating rapid sea level rises by discharging seawater depending on sea level.

In this section, how each channel (Bungo, Kii, and Kanmon) affects sea surface heights and circulations in the SIS was analyzed. Figure 10 illustrates the spatial ranges influenced by Bungo Channel and Kii Channel. Regarding sea surface heights, nearly the entire SIS was affected by the Bungo Channel (indicated by blue and purple areas). However, the Kii Channel had a significant effect only in proximity to the channel (purple and red areas). Circulations exhibited a slightly different pattern with currents at the center of the SIS (Harima Nada) being more influenced by the currents at the Kii Channel. Sea water was usually discharged from the SIS to outside through the Kanmon Strait.