Hakata Bay is a semi-enclosed bay located in Fukuoka City at the northern part of Kyushu Island, Japan. The bay is important because many ships originating in Asian countries transit there. The bay has a rich ecosystem because it has tidal flats such as Wajiro Higata and Imazu Higata in the eastern and western parts of the bay. Moreover, it is important to protect the bay because migratory birds from other continents frequent the bay  . The bay has an urban area near it, and is experiencing various problems owing to negative changes in its ecosystem and environmental degradation of coastal water, such as red tides and local anoxic water masses. Therefore, it is important to protect the environment of Hakata Bay, including its tidal flats.
In recent years, global-level climate change has induced frequent torrential rain. Since floods occur and cause great damages on land, various studies have examined them (e.g.,  ). When there is torrential rain in Fukuoka City or in the region around Hakata Bay, negative environmental changes are often induced, and it is thought that major fresh-water discharges from rivers cause these problems. However, the impact of major inflows into semi-enclosed bays has not been investigated. Hence, it is necessary to construct a model that combines a hydrodynamic model with a hydrologic model. There are the models what do this (e.g.,  ), but there is no combined model that targets torrential rain.
Firstly, in this research, a tank model is used as a hydrologic model. This is often done to calculate discharges of watersheds (e.g.,   ). The tank model is applied to calculate river discharges due to torrential rain. Secondly, the convective-dispersive model of salinity is applied to estimate the impact of major fresh- water inflow due to torrential rain. Hakata Bay has narrow and complicated sea areas such as Imazu Bay and the artificial island called “Island City” in the western and eastern parts of the bay, respectively. A numerical solution with high resolution is essential and effective for estimating and predicting such environments. Nevertheless, a substantial amount of time is needed for calculation, dividing the whole area into small meshes. Nesting techniques are suitable for dealing with this difficulty: wide areas can be covered quickly while finer calculations are conducted in specific narrow areas (e.g.,   ).
In this study, a two-way nesting “edge” technique is developed and applied in the two-dimensional Reynolds depth-averaged equations and the convective- dispersive equation of salinity in Hakata Bay. The model is combined with a tank model to analyze the impact of major fresh-water discharges from rivers due to torrential rain.
2. Materials and Methods
2.1. Governing Equations
A two-dimensional Reynolds depth-averaged model for tidal current and diffusion of salinity was utilized to investigate the physical environment in Hakata Bay. Two-dimensional models are easier to program than three-dimensional models, and the calculation time of the former is shorter than that of the latter. Two-dimensional models are often used for seas of shallow depth, small changes of vertical hydraulic quantities, and large scales in the horizontal direction (e.g.,  ). The governing equations of the model are as follows:
where η (m) is the sea surface elevation and U and V (m/s) are the components of velocity in the x- and y-directions, respectively. Furthermore, t (s) is the time; h (m), the water depth; g (m/s2), the acceleration of gravity; f (s−1), the Coriolis parameter; AH and KH (m2/s), the coefficients of the eddy viscosity and the eddy diffusivity in the horizontal direction, respectively; ρ (kg/m3), the density of sea water; τs = (τsx, τsy) (N/m2), the wind shear stress on the sea surface; τb = (τbx, τby ) (N/m2), the bottom shear stress; and S (psu), salinity.
The shear stresses at the sea surface and at the bottom were calculated from the following equations:
where Uw and Vw (m/s) are the wind velocity components at 10 m above the surface. Furthermore, ρa (1.293 kg/m3) is the air density, Cfs (0.0013) the wind friction coefficient, and Cfb (0.0026) the bottom friction coefficient.
The coefficients of the eddy viscosity and the eddy diffusivity were determined using the Smagorinsky model  :
where Cs (0.10) is the model parameter and Δ (m) is the mesh size.
A finite difference method was utilized to solve the governing equations on a staggered grid. For the continuity equation and the momentum equations, a leapfrog method was used. For the convective-dispersive equation, the split-op- erator approach was applied  . In the first step, only convective terms were calculated using the upwind differencing method. In the second step, diffusive terms were calculated using the alternating direction implicit (ADI) method and substituting the value of the result from the first step. To set the boundary conditions between the land and the sea, a Neumann-type boundary and a no-slip boundary were introduced for salinity and velocities, respectively. To determine the open boundary condition at the entrance of the bay, a tidal elevation was given by a tidal harmonic constant in Tsuyazaki. Moreover, a land mask function (LMF) was adopted in the wet-and-dry method to express the emergence of tidal flats  .
2.2. Edge Nesting
Nesting techniques are calculation techniques that allow us to conduct both low- resolution calculation in large areas and high-resolution calculations in local areas. Wide areas can be covered quickly while finer calculations are conducted in specific narrow areas by using nesting techniques. Various nesting techniques have been proposed (e.g.,   ). In this research, edge-type nesting was developed. The process of edge-type nesting is as follows:
1) Calculation of low-resolution areas;
2) Interpolation of the results of low-resolution areas to provide boundary conditions for high-resolution areas;
3) Calculation of high-resolution areas;
4) Interpolation of the results of high-resolution areas to provide boundary conditions for low-resolution areas;
5) Repeat the procedures in steps 1 to 4 at every time step.
If a calculation point adjoins a land area, rates of change between two neighboring points were used for the water surface elevation, salinity, and velocity parallel to the land area. In this research, furthermore, the time steps in calculations of low-resolution areas and high-resolution areas were the same so as to facilitate program construction.
2.3. Tank Model
The amount of river discharge is usually difficult to observe, especially after torrential rain. The tank model was applied to determine the amount of river discharge. This model was proposed by Sugawara  and has been used widely for analysis in hydrology (e.g.,  ). As can be seen in Figure 1, multiple tanks were aligned vertically and orifices were located below and beside each tank. The parameters expressed the size and height of the orifices, and river discharges were expressed by adjusting the parameters.
The tank models were applied to the six main rivers: the Tatara River, Mikasa River, Naka River, Hii River, Muromi River, and Zuibaiji River, which flow into
Figure 1. Tank model of the Zuibaiji River.
Hakata Bay. Figure 1 shows the structure and parameters of the tank model at the Zuibaiji River. The tank model for the Zuibaiji River has three tanks, and the first tank has three outlets. This is because the model was focused on flood analysis. The parameters were decided by trial and error process. The calculation period was from 1 April 2002 to 31 March 2003. The calculated discharge of the Zuibaiji River was compared with the discharge observed during the same period shown in Figure 2. The Nash?Sutcliffe coefficient (NS = 0.72) and relative error (RE = 0.28) indicate that the calculated discharge accurately reproduced the observed data. Since the characteristics of the other main rivers are similar to those of the Zuibaiji River, all tank models can be assumed similarly to have good reproducibility. Discharges of small rivers were calculated by multiplying each catchment area by the average runoff (mm/hour) of the six main rivers.
In this research, the impact of large discharges from rivers due to torrential rain was analyzed by using a hydrodynamic model with the edge-type nesting technique targeting Hakata Bay. The calculation period was from 11 September 2002 to 21 September 2002 because there was torrential rain on September 16 in Fukuoka City (163.5 mm/d). Major discharges from approximately 40 inflow points, varying in scale owing to the torrential rain, were taken into consideration in the hydrodynamic model. The total amount of all river discharges during the calculation period of the hydrodynamic model is shown in Figure 3. Because of heavy rain on September 16th, the total amount increased drastically to about 60 times more than the discharge before the torrential rain.
2.4. Calculation Area
Figure 4 shows the bathymetry of the low-resolution area. The target bay was Hakata Bay, which is located in Fukuoka City, Kyushu Island, Japan. An artificial island called “Island City” lies in the eastern part of the bay and Imazu Bay lies in the western part of the bay. These two regions are shown as high-resolu- tion areas in Figure 5 and Figure 6. Because the area around “Island City” is
Figure 2. Comparison of observed and calculated discharges from Zuibaiji River from 1 April 2002 to 31 March 2003.
Figure 3. Sum of all river discharges from 11 September 2002 to 22 September 2002.
narrow owing to constructions and Imazu Bay is a small bay with complex geographical features, the two areas must be evaluated using small meshes. The locations of observation stations for velocity and salinity are indicated in Figure 4 and Figure 5. The calculation was conducted with Δx, Δy = 100.0 m in the low-resolution area and Δx, Δy = 25.0 m in the high-resolution areas, and Δt = 0.2 s in both areas. At the entrance of the bay, the tide level calculated by tidal harmonic constants in Tsuyazaki was used as an open boundary condition in the calculation.
Before the simulation of torrential rain, validation of this hydrodynamic model was conducted. The calculation period for validation was from 10 July 2007 to 3 August 2007. Discharges of rivers in Figure 4 and Figure 5 using the tank models were taken into consideration as river inflow.
Figure 4. Bathymetry of the low-resolution area and locations of observation stations and rivers.
Figure 5. Bathymetry of the area around “Island City” and locations of observation stations and rivers.
2.5. Validation of the Model
Figure 7 shows the comparison of observed and calculated tidal current vectors at observation stations shown in Figure 4 and Figure 5 from 0:00 to 23:00 on July 30, 2007 in 1-h intervals. The tidal current vectors showed good agreement with observed vectors. Complicated flow in narrow areas due to the constructions of “Island City” was well reproduced. Figure 8 shows the comparison of observed and calculated salinity at observation stations shown in Figure 4 and
Figure 6. Bathymetry of Imazu Bay and location of the Zuibaiji River.
Figure 7. Comparison of observed and calculated tidal current vectors in one-hour intervals from 0:00 to 23:00 on July 30, 2007.
Figure 8. Comparison of observed and calculated salinity in one- hour interval from 0:00 to 23:00 on July 30, 2007.
Figure 5 from 0:00 to 23:00 on July 30, 2007, in 1-h interval. The salinity showed good agreement with observed values. Therefore, this hydrodynamic model using edge nesting calculations is clearly applicable.
3. Results and Discussion
Figure 9 shows the calculation results for contour lines of salinity in the low- resolution area. It shows the distribution of salinity (a) just before the torrential rain, (b) after the torrential rain, and (c) two days after the torrential rain, respectively. Comparing Figure 9(a) and Figure 9(b), fresh water spread out near the mouths of rivers. However, large changes can be found when Figure 9(c) is examined. Low-salinity water covered the whole of the inner part of Hakata Bay, and water of lower salinity than outer sea water (<34.0 psu) spread out to the bay’s mouth. Fresh water from Tatara River and Mikasa River, in particular, spread out broadly, and low-salinity water covered all of area D, although it did not cover it in panels (a) and (b). Figure 10 shows the calculation results in the high-resolution area of the eastern part of Hakata Bay, and panels (a), (b), and (c) are shown at the same times as in Figure 9. Fresh water spread out and high-salinity water (>34.0 psu) was reduced when comparing Figure 10(a) and Figure 10(b), but large differences can be seen in Figure 10(c). The whole area around “Island City” was covered with low-salinity water due to river discharges within the area and fresh water from the low-resolution area. It was found that fresh water takes more than two days to spread out, and stays in the inner part of the bay for a long time. According to  , an anoxic water mass was induced in
Figure 9. Distribution of salinity (a) at 8:00 on September 16, 2002, (b) at 4:00 on September 17, 2002, and (c) at 22:00 on September 19, 2002 in the low-resolution area.
Figure 10. Distribution of salinity (a) at 8:00 on September 16, 2002, (b) at 4:00 on September 17, 2002, and (c) at 22:00 on September 19, 2002 in the high-resolution area around “Island City”.
the inner part of the Ariake Sea after torrential rain in the summer of 2001. Therefore, it is important to track the behavior of fresh water because it may promote the formation of anoxic water massesowing to large fresh water inflow in Hakata Bay.
Figure 11 shows the calculation results for contour lines of salinity in the high-resolution area of the Imazu Bay in the western part of Hakata Bay. It shows the distribution of salinity (a) just before the torrential rain, and (b) one day after the torrential rain. After the torrential rain, fresh water flowed out of the mouth of Imazu Bay and spread out fullyinImazu Bay. Many short-neck clams reside in the tidal flats of Imazu Bay. Because short-neck clams have difficulty living in water of less than 5 psu  , it is possible that some of these clams in Imazu Bay are influenced after torrential rains. In addition, a significant amount of mud flows into Imazu Bay, especially in torrential rain. It is possible for mud to accumulate in the blue area in Figure 11(b), where fresh water spreads. It is also possible that a substantial amount of mud deposited after torrential rain exerts a major influence on the benthoses in Imazu Bay. Therefore, increasing torrential rain in recent years may have had amajor impact on the ecosystem of Imazu Bay.
In this research, the two-dimensional Reynolds depth-averaged and convective- dispersive equation of salinity with the two-way nesting “edge” technique was developed and applied to the case of Hakata Bay. The calculation results for
Figure 11. Distribution of salinity (a) at 8:00 on September 16, 2002 and (b) at 6:00 on September 18, 2002 in the high-resolution area, Imazu Bay.
tidal current vectors and salinity showed good agreement with the observed data in most of the calculation area. Complicated flows in narrow areas due to the constructions of “Island City” were well reproduced.
Moreover, in this research, the impact of major discharge from rivers due to torrential rain was analyzed by using a hydrodynamic model with the edge nesting technique targeting Hakata Bay. To obtain the rough amount of discharges from rivers, tank models were applied. The discharge calculated using the tank model of Zuibaiji River showed good agreement with the observed discharge. According to the calculation results, two findings were noted. First, fresh water takes more than two days to spread out and stays in the inner part of the bay for a long time. Second, fresh water flows out of the mouth of Imazu Bay and spreads out fully in Imazu Bay, and it is possible that a substantial amount of mud flows and is transported into the bay, similar to the behavior of salinity. These phenomena suggest remarkable fresh water behavior over a period of a few days, and it is possible that large discharges from rivers influence the benthos living in the tidal flats near the river.
The observed hydrographic data and topographical data in Hakata Bay were provided by the Department of Environment and the Department of Port and Harbor of Fukuoka City Government. This research was partially supported by JSPS KAKENHI Grant Number 17K15347.
 Akiyama, J., Shigeeda, M. and Nomura, S. (2013) Numerical Investigation of Characteristics of Flood Flow in the Onga River Basin and Inundation Flow in Iizuka City due to Torrential Rain. Japan Society of Civil Engineers, 69, 1579-1584. (In Japanese with English Summary)
 Sato, Y., Komatsu, E., Nagare, H., Uehara, H., Yuasa, T., Okubo, T., Okamoto, T. and Kim, J. (2011) Construction and Validation of Lake Biwa Basin Simulation Model with Integration of Three Components of Land, Lake Flow, and Lake Ecosystem. Journal of Japan Society on Water Environment, 34, 125-141.
 Hayashi, S. and Yamada, T. (2015) Effective Rainfall Amount and Slope Surface Failure Caused by Heavy Precipitation—Study on the Seepage Flow in the Soil La-yer of a Slope. Journal of the Japanese Geotechnical Society, 10, 157-162.
 Rahman, M.M., Paul, G.C. and Hoque, A. (2013) Nested Numerical Scheme in a Polar Coordinate Shallow Water Model for the Coast of Bangladesh. Journal of Coastal Conservation, 17, 37-47.
 Debreu, L., Marchesiello, P., Penven, P. and Cambon, G. (2012) Two-Way Nesting in Split-Explicit Ocean Models: Algorithms, Implementation and Validation. Ocean Modelling, 49-50, 1-21.
 Hu, B.S. and Kot, S.C. (1997) Numerical Model of Tides in Pearl River Estuary with Moving Boundary. Journal of Hydraulic Engineering, 123, 21-29.
 Uchiyama, Y. (2004) Wetting and Drying Scheme Based on a Modified Logarithmic Law for Three-Dimensional Terrain-Following Coastal Ocean Models. Annual Journal of Coastal Engineering, Japan Society of Civil Engineers, 51, 351-355. (In Japanese with English Summary)
 Nihei, Y., Sato, K., Nadaoka, K., Kumano, R. and Nishimura, T. (2003) Development of a New Multi-Nesting Approach for Coastal Current Simulation. Proceedings of Japan Society of Civil Engineers, 740, 171-183.
 Martinsen, E.A. and Engedahl, H. (1987) Implementation and Testing of a Lateral Boundary Scheme as an Open Boundary Condition in a Barotropic Ocean Model. Coatal Engineering, 11, 603-627.
 Sugawara, M. (1985) Tank Model—For the Derivation of River Discharge from Rainfall. Journal of Geography, 94, 209-221. (In Japanese with English Summary)
 Noto, F., Maruyama, T., Hayase, Y., Takimoto, H. and Nakamura, K. (2010) An Evaluation of Snow Storage Depth in the Tedori River Basin Using Tank Model. Japanese Society of Irrigation, Drainage and Rural Engineering, 78, 265-271.
 Tsutsumi, H., Okamura, E., Ogawa, M., Takahashi, T., Yamaguchi, H., Montani, S., Kohashi, N., Adachi, T. and Komatsu, T. (2003) Studies of the Cross Section of Water in the Innermost Areas of Ariake Bay with the Recent Occurrence of Hypoxic Water and Red Tide. Oceanography in Japan, 12, 291-305. (In Japanese with English Summary)
 Nakamura, M., Shinagawa, A., Toda, K. and Nakao, S. (1997) Tolerance of 4 Bivalves from Lakes Shinji and Nakaumi to Environmental Factors. Fish Farming, 45, 179-185. (In Japanese with English Summary)