Seasonal Study of the Kako River Discharge Dynamics into Harima Nada Using a Coupled Atmospheric–Marine Model

Abstract

This study developed a coupled atmospheric–marine model using the COAWST model system for the Harima Nada area between spring 2010 and winter 2011 to evaluate the seasonal influence of the Kako River’s discharge in the sea. The Kako River is one of the largest rivers in southwest Japan, contributing almost half of the freshwater discharged in the Harima Nada region in the Seto Inland Sea. Validation was conducted for the entire period, showing a good performance for the atmospheric and marine variables selected. Multiple experiments injecting an inert tracer in the Kako River estuary were performed to simulate the seasonal river water distribution from the estuary into the sea and to analyze the seasonal differences in concentration patterns and mean residence times in Harima Nada. Because the study area is shallow, the results were evaluated at the surface and 10 m depth layers and showed significant seasonal differences in tracer distribution, circulation patterns, and mean residence times for the region. On the other hand, differences seemed to not be significant during the same season at different depths. The obtained results also agreed with the area’s natural water circulation, showing that the Kako River waters tend to distribute towards the west coast of Harima Nada in the warmer seasons but shift towards the east in winter. The influence of the Kako River in the center of the study area is seasonal and strongly dependent on the direction of the horizontal velocities more than their magnitude. The mean residence times varied seasonally from approximately 30 days in spring to 12 days in fall. The magnitude of the horizontal velocity was found to be maximum during summer when circulation patterns at the surface and 10 m depth in the central part of Harima Nada also seem to promote the strongest horizontal and vertical mixes.

1. Introduction

During the past decades, the main environmental problems that coastal management has been focusing on are as follows: (1) the remediation of waters and ecosystems and (2) the prevention and maintenance of its waters’ quality and biodiversity. Both problems are strongly related to human activities and the socio-economic development and growth of a country or region [1,2,3]. All human activities performed on the coast or inland inside a river watershed have some extent of influence in coastal areas [3,4], and the prevention and minimization of environmental risks inland and on the shorelines are essential good practices for attenuating and preventing coastal deterioration [3,5,6].

The Seto Inland Sea is Japan’s largest semi-enclosed sea, with over 700 small islands and many straits, bays, and basins/nadas in an entire range of sizes [7]. Combined with the characteristics above, the mountainous geography of the shoreline is responsible for the special climatological characteristics of the region and its high biodiversity. Around 25% of the Japanese population lives in the surrounding coastal areas, being the most concentrated and industrially developed areas in Osaka and Fukuoka Prefectures [7]. During the past five decades, after the boom in the Japanese economy, the entire area has been affected by severe environmental issues. First, it was the eutrophication of the waters and its associated problems (continuous red tide problems, loss of many pelagic species, changes in phytoplankton populations, etc.); later, after remediation, the oligotrophication processes that developed in some of its regions provoked a sustained reduction of fish catchment and seaweed production, the discoloration of some seaweed species, and biodiversity loss among other issues [7,8,9,10,11,12].

Harima Nada is one of the eastern Seto Inland Sea basins (Figure 1b), where the initial eutrophication problems shifted into oligotrophication issues in the past decades. Because of its enclosed geography and shallow waters (~50 m in the central area), the region is easily affected by the discharge of numerous rivers [8,9]. Rivers link the inland activities of a region with the coastline by transporting in with their waters and sediments a variety of pollutants and nutrients to the sea, and the larger the average river discharge is, the higher the impact of the river on the water quality and ecosystems of the coastal areas.

The most important river in Harima Nada is the Kako River, located in the center of the northern coast. The Kako River accounts for around 40% of the total freshwater discharged in Harima Nada, and it is the higher contributor of major land-derived nutrients like total nitrogen (TN) in the area [12,13,14].

Much research has been performed in the entire Seto Inland Sea [7,8,9,15,16,17,18,19,20], Harima Nada [10,11,12,13,14,21,22,23,24,25,26], and its surrounding areas like Osaka Bay [14,27,28], and different numerical and simulation models and monitoring observations have been proposed, applied, and conducted to describe and study the characteristics of these regions. Until this date, all the reported ocean dynamic models for the area used offline atmospherical forcing, which means that the input of atmospheric variables was either previously calculated using an atmospheric model, belonged to reanalysis datasets, or consisted of a time series of observed data.

The use of reanalysis datasets (that are also datasets for boundary and forcing conditions in dynamic atmospheric simulations) is the most extended offline forcing in marine and ocean simulations but has the main problem of the scale of the dataset, which can be inconvenient in small areas and calculation domains. In addition, reanalysis and single dynamic atmospheric simulations do not consider differences in the sea surface height (SSH) due to tides over the water bodies, which may increase the error in sea surface temperature determinations in the study of specific phenomena. One of the ways this issue can be addressed is by coupling the atmospheric and marine models, which enhances the accuracy of both simulation model results and gives a more realistic approach to reality.

This work used the Coupled Ocean-Atmospheric-Wave-Sediment-Transport (COAWST) Modeling System developed by the United States Geological Survey (USGS) [29,30] to numerically evaluate the dynamics and nutrient distribution discharged from the Kako River into Harima Nada waters using its atmospheric (WRF—Weather Research and Forecasting Model) and ocean (ROMS—Regional Ocean Modeling System) components in a coupled form. This approach has not been reported up to date in the Seto Inland Sea or the Harima Nada region, and it offers the possibility of improving the quality of the results for atmospheric and marine simulation results, as well as the study of particular meteorological phenomena and their effects on the waters of the area.

This work aims to present a coupled WRF-ROMS model and its validation to study the circulation dynamics of riverine waters in Harima Nada as the first step in the development of a biogeochemical model for the area and its surroundings. For analyzing the circulation dynamics of the region, several inert tracer simulations were run to study the Kako River’s discharged water distribution, mean residence time, and surface water age seasonally in the region. Because of the present oligotrophication problems in Harima Nada, estimating these parameters is important in defining and characterizing biological and chemical environmental processes and their parameters for developing numerical models that describe nutrients and pollutant cycles.

2. Methodology

2.1. The Models

The USGS developed the COAWST Modeling System to help study coastal processes and their changes, evaluate the consequences after coastal natural disasters, and study extreme meteorological events and their significance [29,30]. The modeling system distribution combines multiple high-resolution models and tools for enhancing its full capability (models’ coupling tool, models’ grid interpolation, etc.), allowing the user to develop single or coupled model applications. In this work, version 3.7 of COAWST was used, which distributes the Model Coupling Tool (MCT) v 2.6.0, ROMS v 3.9 (Rutgers), WRF Model v 4.2.2, and models for wave and sediment simulations (https://code.usgs.gov/coawstmodel/COAWST/-/releases, accessed on 26 October 2023).

Figure 2 shows the model layout and the flowchart of input, exchanged, and output variables in the presented model setup.

2.1.1. The Atmospheric Model

The WRF model [31] is a 3D terrain-following model used for atmospheric research and forecasting applications in a wide range of horizontal and vertical resolutions. It allows a variety of physical and grid parametrizations to personalize its applications, predicting a wide range of atmospheric variables (wind momentum components, pressure, air temperature, long and shortwave radiative fluxes, rain, etc.). The WRF model distributed with the COAWST Modeling System has some modifications to its source code to better define variables over the ocean surface [30].

2.1.2. The Marine Model

The ROMS is a free-surface terrain-following model that uses a combination of specially adjusted time-stepping algorithms for solving the 3D Reynolds-averaged Navier–Stokes equations using hydrostatic approximations [32,33,34,35]. It allows the personalization of many model components and options for boundary conditions, advection scheme algorithms, biochemical models, etc. [30].

2.2. Model Configuration

Table 1. WRF and ROMS models’ setup and summary of the river forcing and tracer experimental conditions.

WRF Model
Analysis period	March 2010–February 2011
Spin-up	8 months (June 2009–February 2010)
Horizontal grid cells (X × Y)	120 × 120
Horizontal resolution	6.0 km
Vertical layers	30
Time step	25 s
Boundary and initial conditions	ERA5
ROMS model
Analysis period	March 2010–February 2011
Spin-up	8 months (June 2009–February 2010)
Horizontal grid cells (X × Y)	200 × 100
Horizontal resolution	3.0 km
Vertical layers	32
Baroclinic time step	60 s
Boundary and initial conditions	HYCOM
Tidal boundary conditions	TPXO9 Atlas v.5
Rivers in the domain	11 (5 northern coasts, 2 center, 4 in the southern coast)
Rivers considered for calculation	1 (Kako River, largest freshwater contribution)
River input algorithm	Volume vertical influx
River and tracer experimental conditions
Analysis period	March 2010–February 2011
Spin-up	12 months (March 2009–February 2010)
Horizontal resolution	1 km
River discharge temporal resolution	1 h
Vertical levels of the estuary	32 layers
River vertical fractional distribution (Temperature, salinity, tracer)	Homogeneous
River temperature	15 °C, constant for the entire period
River salinity	0.0 PSU
Tracer concentration	- 0.1 kg/m3 (for seasonal distribution analysis) - 0.5 kg/m3 (for mean residence time seasonal analysis)
Tracer discharge injection	- Entire period for seasonal distribution analysis - 24 h for mean residence time in seasonal analysis

summarizes the WRF and ROMS model setup, river forcing, and tracers’ experimental conditions.

The models’ computational domains are shown in Figure 1a. The WRF domain covers Central and West Japan, with a horizontal resolution of 6.0 km and an extension of 120 cells from S-N and W-E, respectively. The number of sigma-pressure coordinated layers selected from the surface to 100 hPa was 30. The WRF model physic options were chosen from the study by Chatani et al. [36] since they have been proven effective in describing the region.

The ROMS computational domain was chosen to be centered on WRF’s domain, covering the entire Seto Inland Sea and extending to the Pacific Ocean approximately 60 km from the Bungo Channel (the western channel that communicates the Seto Inland Sea with the Pacific Ocean, Figure 1a). The selected horizontal resolution was 3.0 km for the entire area, extending 100 and 200 cells from S-N and W-E, respectively. A nudged-open boundary condition following the Orlanski radiation approach [37] was selected for three-dimensional momentum, temperature, and salinity.

Simulation analysis was conducted from March 2010 to February 2011, with eight months of spin-up for the WRF and ROMS models and twelve months for river discharge. Because Japan is a four-season country, seasonal experiments with inert tracers were conducted to estimate the seasonal distribution of Kako River waters and their mean residence time in the Harima Nada area. Seasons were considered as follows: March–May (spring), June–August (summer), September–November (fall), and December–February (winter). This period was chosen to evaluate the model’s performance and validate its results because it is a period of controlled conditions. The availability of observed datasets on streamflow and TN load discharged into Harima Nada, as well as the previous validation of the applied river model [38], reinforces the significance of the results obtained. As an additional consideration, the period of time corresponded to the average precipitation conditions in the area, as detailed in the study by Pintos Andreoli et al. [38].

Ten rivers directly discharge their waters in the study area (five on the northern coast, three in the southern one, and two in the center region), plus the biggest river of Shikoku Island discharging next to the Naruto Strait, one of Harima Nada’s limits (Figure 1b). Only the influence of the Kako River, which is responsible for the larger discharge of freshwater in the study region [13], was considered in this work. In Japan, the seasonal occurrence of heavy rain events is common, and their effects are very significant on Japanese rivers’ nutrient loads and discharged volumes. For the reason above, the Kako River was forced in the ROMS using a volume vertical influx algorithm to avoid high velocities on the estuary point that can lead to CFL instabilities. Because the river estuary is located in a “shallow” area (around 10 m depth), temperature, salinity, and tracers were homogeneously distributed in the vertical water column.

Two different tracer analyses were conducted to study the dynamics of the Kako River waters. The first analysis involved the injection of a tracer in the Kako River estuary throughout the entire simulation period (including the first eight months of spin-up) to study the seasonal distribution of discharged freshwater and nutrients in Harima Nada. To better understand nutrient concentration distribution from the Kako River into the sea, an inert passive tracer with a concentration equal to the average river’s concentration of TN (1.0 kg/m³) was added to the estuary point of the river. The second tracer analysis was performed to study mean seasonal residence times and tracer concentration distribution. A control domain delimited by the points H09, H27, H29, and H30 (Figure 1b) was used to study tracer distribution after 24 h of constant 0.5 kg/m³ concentration injection in the river estuary on the first day of each season. Tracer concentration and water age were followed in the delimiting points of the control domain and H01. The state of total tracer dispersion was defined when its concentration became lower than 1 × 10⁻⁶ kg/m³ (three orders of magnitude smaller than the maximum concentration value discharged). Tracer was followed for 60 days in all the chosen points, but it was found that after the first 30–35 days, it had totally flushed from the control domain.

Because the bathymetry of the entire Harima Nada is particularly shallow (less than 50 m in the deepest areas), it was chosen to evaluate tracer distribution and river water distribution at the surface and 10 m depth. The 10 m layer is deep enough to be separate from the surface layer, and because the average bathymetry is in the order of 10–15 m, it can be used to describe the middle depth in the largest part of the area. The coastal bathymetry in Harima Nada is around 10 m deep or lower at many points; so, in those cases, it was considered that the 10 m layer was representative of the bottom layer.

The seasonal mean residence times were estimated using the mean water age method developed by Zhang et al. [39]. This approach is based on the constituent-oriented age and residence-time theory [40,41,42,43,44] and uses a set of two tracers that are released together to estimate the mean age of the water. The mean tracer age is computed by calculating the mass-weight arithmetic average time since the tracer left the source point and the mean residence time, defined as the time needed for the tracer to leave a control domain of interest. It was demonstrated by Zhang et al. [39] that the mean residence time calculated using this method agrees with the residence time, but their approach and physical significance are different. The mean water age method was selected because this work is the 1st step in the development of a biogeochemical model for the Harima Nada area, and the concept of water age is more significant in the selection of parameters for simulating biogeochemical processes than the residence time (which is considerably more important in the study of ecological dynamics). The mean residence times ($\overline{\tau })$ were calculated using Equation (1), as the relationship between the mean age in a location x (${\alpha }_x)$ and the concentration of tracer at the same location ($C_x)$.

$$\overline{\tau }=\frac{{\alpha }_x}{C_x}$$

(1)

Validation was conducted for the atmospheric and marine models in both domains. Figure 1b shows the study area and the location of the meteorological and sea observatory stations used to validate the model results. Additionally, it shows the selected control domain for studying mean residence time and the points used to analyze tracer concentration peaks and their time distribution.

2.3. Datasets

Input forcing data and boundary conditions used for the WRF were from ERA5, the fifth-generation European Center for Medium-Range Weather Forecast (ECMWF) reanalysis [45]. The ROMS 3D boundary conditions, initial conditions, and climatology fields were extracted from the Hybrid Coordinate Ocean Model (HYCOM) [46,47]. Tidal 2D boundary conditions were processed by extracting tidal harmonic components from the TPXO9 Atlas v.5, a barotropic tide model developed by Egbert and Erofeeva [48]. The rest of the necessary meteorological fields for the ROMS were forced from the WRF results because of the selected coupled model configuration.

The Kako River discharge was calculated using a physically distributed hydrological model in the entire watershed. Pintos Andreoli et al. [38] provided the river model’s description, validation, and verification. Additional hydro-hyetographs are given in Supplementary Materials to show the river model performance (Figure S1 and Table S1).

Observed datasets between March 2010 and February 2011 were used to validate the models’ performance and the obtained results. For the verification of the atmospheric fields, datasets from the Meteorological Observatories and Automated Meteorological Data Acquisition System (AMeDAS) [49] that surround the studied area were used. For the ROMS validation, monthly observed data on sea temperature and salinity at the surface, middle, and bottom layers from the observatory stations of the Hyogo Prefecture [50] (Figure 1b) were used. Additionally, SSH was verified by using biweekly observed data from the JMA coastal observatories [49] at different points of the domain (A, B, C, D, E, F in Figure 1a) to check on the model’s physical results’ consistency.

3. Results and Discussion

3.1. Model Validation

3.1.1. Atmospheric Model Validation

Validation was performed for some of the most common meteorological variables used to force the ROMS model. Because the study region is not particularly large in extension and there are no considerable geographic differences that can affect meteorological variables from one point to another, it was chosen to compare the average values. Observed datasets from three AMeDAS stations [49] (Figure 1b) located in the surroundings of Harima Nada were averaged and compared with simulation results (punctual validation in each station is shown in Table S2 in Supplementary Materials). Additionally, because the study site is comparably smaller than the domain’s size, punctual validation of the simulation results was performed in other meteorological observatories (Table S3). Figure 3 shows the observed and simulated average values of the zonal and meridional wind components, surface pressure, temperature, and specific humidity in the study area, and

Table 2. Statistical results of WRF simulation results in the Harima Nada area.

	Observed Average	Simulated Average	R	pBIAS	MAE	RMSE	IA
u-wind (m/s)	0.74	0.92	0.96	---	0.49	0.62	0.93
v-wind (m/s)	−0.26	−0.41	0.72	---	1.05	1.31	0.62
Wind speed (m/s)	2.9	3.2	0.80	---	0.53	0.73	0.83
Surface pressure (hPa)	1010.3	1010.5	~1.0	<0.1	0.37	0.50	~1.0
Temperature (°C)	--- *	--- *	0.97	0.3	0.65	0.82	0.99
Specific humidity (g/kg)	--- *	--- *	0.97	4.1	0.48	0.64	0.99

* Average values of temperature and specific humidity are not shown due to their strong seasonal variability. R: correlation; pBIAS: percentage of BIAS; MAE: mean absolute error; RMSE: root mean square error; IA: index of agreement. Benchmarks: wind—RMSE < 2.0; 0.6 < IA. Temperature BIAS ≤ ±2.0°; MAE < 2.0; 0.8 < IA. Specific humidity—|MAE| < 2.0; 0.6 < IA.

summarizes the obtained statistics used to validate the results and evaluate the model’s performance. Average values of air temperature and specific humidity are not shown because of their large seasonal variability.

Atmospheric and ocean systems’ variables are intrinsically uncertain, and their turbulent and non-linear characteristics, the limitations of the equations that explain their phenomena, and the wide time scale of the processes involved limit their forecasting and predictability [51]. Some variables, like temperature, can be well described for statistical distributions and considered statistically normal at the considered temporal scales in this work. In contrast, other variables like rain, wind, SSH, or waves have a higher random probability distribution and are affected by more significant errors in their determination. Because of this, a set of statistics considered suitable for the evaluation and description of the results was chosen separately for each variable. Reference benchmarks for wind components, temperature, and specific humidity [52,53] were used to evaluate simulation results. Very good reproducibility was found for the zonal wind component (u-wind), temperature, and specific humidity. The model showed slightly lower performance in reproducing the meridional wind component (v-wind) and the wind speed for the area, which can be associated with the complex geography of the evaluated region. In the case of the surface pressure, the validity of the benchmarks was extended since the model’s results were found to be very good at reproducing the observations. For all the chosen variables, the statistical results of the analyzed variables agreed with the bibliography benchmarks, validating the proposed atmospheric model for the study area.

3.1.2. Marine Model Validation

Daily average water temperature and salinity simulation results were punctually validated in the Harima Nada area in the stations shown in Figure 1b (H28, H01, H11, and H13). These two variables are considered active tracers in the ROMS model, and they are very important parameters in any marine/oceanic system. The reason why it was chosen to validate the results of the tracer experiments was that both types of tracers (active and passive ones) were simulated under the same advection scheme. Figure 4 shows the calculated and observed time series for both variables at four different points at the surface, middle, and bottom layers. The depth of the layers varies from point to point (see

Table 3. Statistical results of ROMS simulated variables in the Harima Nada area.

Temperature
		R	pBIAS	MAE		RMSE		IA
H01	sfc (~0.5 m)	0.99	9.4	1.73		2.05		0.98
	mid (~10.0 m)	~1.00	6.1	1.07		1.48		0.99
	bott (~20.0 m)	~1.00	4.2	0.74		0.88		~1.00
H28	sfc (0.5 m)	0.99	8.6	1.02		1.97		0.98
	mid (~5.0 m)	~1.00	5.3	1.17		1.09		0.99
	bott (~10.0 m)	~1.00	4.3	0.59		0.86		~1.00
H11	sfc (~0.5 m)	~1.00	5.4	0.79		1.24		0.99
	mid (~10.0 m)	0.99	6.5	0.89		1.49		0.99
	bott (~32.0 m)	0.99	1.6	0.73		0.71		~1.00
H13	sfc (0.5 m)	0.99	3.5	1.54		0.97		0.99
	mid (~10.0 m)	~1.00	4.4	0.91		1.08		0.99
	bott (~32.0 m)	0.99	3.7	0.73		1.02		0.99
Salinity
		Observed average (PSU)	Simulated average (PSU)	R	pBIAS	MAE	RMSE	IA
H01	sfc (~0.5 m)	31.0	32.3	0.90	−4.2	1.30	1.62	0.71
	mid (~10.0 m)	31.8	32.9	0.83	−3.4	1.07	1.48	0.61
	bott (~20.0 m)	31.9	33.0	0.78	−3.4	1.07	1.12	0.56
H28	sfc (0.5 m)	30.6	31.8	0.76	−3.9	1.20	1.44	0.71
	mid (~5.0 m)	31.5	32.2	0.82	−2.2	0.71	0.81	0.71
	bott (~10.0 m)	31.6	32.2	0.81	−2.0	0.63	0.73	0.72
H11	sfc (~0.5 m)	32.0	33.2	0.84	−3.6	1.16	1.21	0.54
	mid (~10.0 m)	32.1	33.3	0.86	−3.9	1.24	1.27	0.50
	bott (~32.0 m)	32.3	33.5	0.83	−3.7	1.19	1.23	0.48
H13	sfc (0.5 m)	31.9	33.2	0.91	−4.0	1.27	1.29	0.47
	mid (~10.0 m)	31.9	33.4	0.89	−4.5	1.45	1.48	0.41
	bott (~32.0 m)	31.0	31.0	0.87	−4.0	1.28	1.31	0.40

R: correlation; pBIAS: percentage of BIAS; MAE: mean absolute error; RMSE: root mean square error; IA: index of agreement. Average values for water temperature are not shown due to their strong seasonal variability.

) with the corresponding depths of the observed datasets because of differences in the local bathymetry of the study area [50]. The validation points were selected because of their location and the data availability at three layers of depth (most observatories only collect time series of surface data). Stations H01 and H28 were selected because of their proximity to the Kako River estuary. The observatories H11 and H13 are far from the coast and the river estuary, both in the deepest area of Harima Nada. An additional consideration about the validation of model results is that observed data series with bi-weekly and monthly temporal resolution may be unsuitable for efficiently evaluating the model’s performance and the accuracy of the results for certain variables, in particular those with very high dependence on other variables or with high temporal variability.

The statistical analysis of the results in Table 3 shows that the model’s performance is high and can describe the water temperature in the region. The correlation coefficients at all depths and the indexes of agreement were all found to be close to one. The positive values of pBIAS show (as seen in Figure 4) that the model tends to slightly underestimate water temperature, particularly in the surface layer and near the coast during the months of summer, which can be a consequence of the model setup conditions for the river forcing (Table 1). Rivers in Japan are particularly shallow, and they run very fast, particularly during the months of spring and summer when rainy and typhoon seasons occur. In addition, temperatures during the end of spring and summer are particularly high in West Japan. In this first approach, the water temperature was considered constant during the entire year, as the average between all the months, an assumption that can affect the model’s performance. Also, the effect of this assumption can be observed in the temperature graphs from the closest observatories to the river estuary (H28), where the differences between simulated and observed values are more pronounced. On the other hand, differences are smaller for those observatories in the middle of Harima Nada (H11 and H13), and water temperature seems to not be affected by the river’s freshwater discharge.

The results for salinity showed that the model can effectively describe the region’s conditions. Still, if some extra considerations in the model setup are addressed, the model’s performance can be improved. Harima Nada is a very enclosed area in the Seto Inland Sea, and the strong effect of river discharges has already been studied and demonstrated in previous research [8,9]. The Kako River is the most important contributor of freshwater in the study area (approximately 40% of total freshwater), but it is not the only one. The other rivers that discharge in the area also impact the salinity and nutrient values in Harima Nada, particularly during the rainy and typhoon season from May to September.

The negative values of the pBIAS show that the model overestimates salinity values in the entire region, which can also be seen in the MAE values. This overestimation is minor and in the range of 1 to 2 PSU units, as can be seen in the salinity graphs in Figure 4. The obtained correlation values were good for all the observatories, showing that the model can describe the salinity trends in the entire area at the evaluated levels. However, the values of IA are acceptable for the two observatories near the river estuary (H01 and H28) but show that the model’s performance is lower in the center of the study region, which can also be seen in the salinity graphs of the stations H11 and H13.

The observed values show that the surface layer of water is the one that shows higher variations in salinity. In contrast, the middle and bottom layers behave similarly and are not subjected to strong changes during the year. In Harima Nada’s center region (observatories H11 and H13), there are no significant differences in the salinity at different depth levels, contrary to what happens near the coast where the effect of the Kako River discharge seems to be of considerable importance. During the months of rainy and typhoon seasons (from spring to summer), the salinity in the center of the study region is lower than the calculated value, which indicates that higher volumes of fresh water or lower salinity water enter the study area because of the occurrence of intense precipitations and an increase in freshwater discharged from land to the sea.

Additionally, relative SSH values were compared with observed data in many domain points to evaluate the model’s physical consistency (Figure 1a). For relative SSH values’ comparison, only R and MAE statistics were calculated to evaluate the trends’ behavior.

From the results shown in Figure 5, it can be seen that the model is good at reproducing relative SSH trends (R~0.7) with errors that are in the order of 10 to 15%. Because of the implicit difficulty of simulating waves and tides’ phenomena and considering the complicated geography of the study site, these results validated SSH performance in the Seto Inland Sea in this first approach.

3.2. The Kako River Discharge Dynamics

Figure 6 shows the average seasonal distribution of the inert tracer at the surface level and the average horizontal velocity in each season at two different layers, surface and 10 m depth. The obtained average surface horizontal velocity values agreed with the global water circulation in the Seto Inland Sea (from west to east during the entire year), showing that Harima Nada’s average influx of water comes from the Bisan Strait in all the seasons, while output fluxes mainly occur from the Akashi and Naruto Straits, towards Osaka Bay and the Kii Channel, respectively.

The more significant seasonal differences in the average magnitude and direction of the horizontal velocity can be found in the central zone of Harima Nada, while in the west and east sides, the conditions do not seem to be seasonally affected. The maximum average magnitude values of speed occur during the summer in the central area, reaching values up to 0.3 m/s. On the other hand, the horizontal velocity magnitudes are approximately 20% lower during the rest of the seasons. Changes in the direction patterns of the horizontal velocity occur during summer in the central region and in winter in the entire Harima Nada. In summer, high velocities and circular patterns in the central area promote water mixing at all levels, as can be seen in Figure 6. In winter, it was found that circulation is strongly driven in the direction W-E. During winter, this W-E pattern extends up to the northern coastal region, where the circulation is E-W for the rest of the seasons.

The seasonal distribution of the inert tracer agreed with the calculated horizontal velocity directions and the observed data for many environmental variables in the Harima Nada region. Previous research and routine analysis in the area have shown that the highest concentrations of nutrients and biomass growth rates occur towards the northern coasts of Harima Nada in the months of spring and summer [54], which agrees with the results of the seasonal distribution of the river’s discharged water.

The tracer distributes from the Kako River estuary to the western coast from spring to fall and changes in winter, following the pattern of horizontal velocity. The maximum and minimal dispersions occurred in summer and winter, respectively (Figure 6), which agrees with the higher and lower river discharges. During spring, distribution is primarily centered in the Kako River estuary towards the coast, but in summer, the divergence on horizontal velocities in the central area of Harima Nada strongly promotes the dispersion and mixture of freshwater towards the west and the center of the study region more efficiently than during other seasons. In the months of fall, when water starts cooling, and the vertical mix becomes stronger because of the action of convection mechanisms, tracer distribution is centered in the river estuary and homogeneously mixes in the entire Harima Nada area. This situation changes in winter, when river discharges are the lowest, and discharged water distributes towards the west, dispersing the tracer towards Osaka Bay.

The same analysis was performed at 10 m depth (second column of Figure 6), showing that the tracer distributes mainly in the same manner as on the surface layer. The same trends were observed seasonally for average horizontal velocities and concentration. While higher concentration values were around 5 × 10⁻³ kg/m³ on the surface level, at 10 m depth, concentration was reduced by approximately 30% in all seasons. In addition, it was also verified that vertical mix is more important in summer and fall when the horizontal velocity or convection mechanisms are stronger. The pattern of horizontal velocity directions is almost the same in both layers in the same season, but the magnitude values are 20% smaller when depth increases by 10 m.

Table 4. Results of 60 days of seasonal tracer simulations in Harima Nada at the surface layer.

	Spring	Summer	Fall	Winter
Mean residence time (days)	28	20	12	25
Maximum concentration (kg/m3)	2.4 × 10−3	5.0 × 10−3	4.4 × 10−3	3.3 × 10−3
Location and days where the concentration peak is reached last	H09 9 days	H09 23 days	H30 8 days	H09 Never

summarizes the obtained values of mean residence times, maximum concentration values, location, and days to reach the peak in tracer concentration in the selected control domain to follow the seasonal injection of the inert tracer. The tracer’s temporal concentration distributions in the control domain results shown in Figure 7 coincide with the calculated mean residence times shown in Table 4. It was found that the longer mean residence times (25 and 28 days) occur during winter and spring when water is colder and circulates slower. Additionally, maximum concentration values were found to be the lowest during these two seasons. The mean residence time in summer is 20 days, around 30% lower than in spring. As previously mentioned, the water circulation pattern during summer in Harima Nada is very different than that in the rest of the seasons and is the main cause of horizontal mix in the central area, but also of river water dispersion towards the coast. This particular circulation and the discharge of more significant freshwater in summer can be responsible for the highest seasonal value of tracer concentration during this analysis in the entire control domain. In the fall, the lowest mean residence time in the area (12 days) happens. The temperature difference makes vertical mix mechanisms more important during the fall. Additionally, surface winds have mostly E-S direction too (Figure 3), dispersing tracer faster towards Osaka Bay, same as in winter, as seen in Figure 6.

Figure 7 shows the results of seasonal tracer concentration and the remaining tracer fraction (ratio between maximum tracer concentration and the daily value in the point) at the selected points in the surface layer. It was found that for the points located in the center (H09) and eastern parts of Harima Nada (H27), tracer concentrations were lower in all the seasons (one order of magnitude) than for those points located closer to the coast and towards the west, which agrees with the seasonal tracer distribution shown in Figure 6. For the southern point of the control domain (H09), tracer concentrations were found to be the lowest in all the seasons, with maximum values of approximately 1 × 10⁻⁴ kg/m³ during summer and a minimum of 0.0 kg/m³ during the entire winter months. Additionally, this point showed a more extensive range of time distributions to reach the concentration peak between seasons. Maximum concentrations are reached between 5 and 10 days in the fall and spring, but it takes the trace more than 20 days during summer to reach the peak concentration in the center of Harima Nada. Except for the small peak around the fifth day of summer, the water circulation affects tracer distribution more significantly in the center of Harima Nada. On the other hand, there is no tracer distribution towards the center of the study area in winter, which also confirms the seasonal observations.

In the eastern point of the control domain (H27), it was found that a higher concentration can be found during winter. Also, this point shows a higher time distribution of tracer remains fraction during winter, indicating more river water moving towards Osaka Bay. During spring and summer, concentration is around half the maximum winter value, which also agrees with the seasonal change in water circulation towards the west in Harima Nada. The peaks of maximum concentration occur between the fourth and seventh day after the tracer was injected in all the seasons and follow similar concentration and ratio distribution patterns in cold (winter and spring) and warm (summer and fall) water seasons.

The analysis of point H01, located next to the river estuary, shows higher concentration values (5.0 × 10⁻³ kg/m³ in summer and the fall) and shorter times for reaching the concentration peaks at all seasons (1 to 2 days). In the case of points H29 and H30, located from the estuary towards the west side, the maximum tracer concentration is reached in winter in the shortest time (2–3 days) for both of them. In both points, concentration distribution patterns are very similar in the pairs winter–spring and summer–fall, but concentrations are in the order of 3-times smaller in the western point (H30).

4. Conclusions

In this work, the seasonal description of Harima Nada water dynamics was successfully accomplished by using for the first time a coupled WRF-ROMS model based on the COAWST Model System. The atmospheric–marine model for the Seto Inland Sea was validated regarding its components in Harima Nada, an area of environmental interest because of the oligotrophication processes that affect it, showing good performance at the selected horizontal resolution for atmospheric and marine variables.

The model was used to study the seasonal distribution of the Kako River waters and its mean residence times at two different depths using tracer simulations. The obtained results showed significant differences between seasons in water circulation and mean residence times, but no important differences at two depths were found during the same season, most likely because of the shallow bathymetry of the study area. The seasonal distribution patterns of river water agreed with the biomass production and nutrient distribution in Harima Nada, reinforcing the importance of rivers’ water discharge and nutrient loads in the area.

The obtained mean residence times were found to be consistent with previous research. They showed strong seasonal variations and strong dependence on surface winds, circulation patterns, and water temperature, going from almost a month in spring to a little more than a week in fall. In addition, the distribution and circulation patterns were found to be similar between layers, and a 10 m increment in depth reduced the tracer concentration by around 30% and the magnitude of the horizontal speed by 20% in the entire study region.

The presented approach showed to have some limitations regarding the estimation of marine variables. The model setup needs some adjustments to correct temperature estimation in the coastal areas that are shallow and showed higher errors in the estimation during summer, particularly at the surface layer. Time series for the river water temperature needs to be considered as an additional input due to its high seasonal variability. In the case of salinity, the other rivers of the Harima Nada area need to be considered to improve the estimation in the central area at all depths. The Kako River is the most important river, but the comparison of observations and simulation results showed that the influence of the other rivers is very important during the entire year, and it may be affected considerably by seasonal circulation patterns. Additionally, the vertical distribution of the river variables may need some adjustment for the Kako River and all the rivers of Harima Nada at higher resolutions based on the results obtained. The estuary area has a shallow bathymetry, but a better understanding of the river plume is needed in order to assess this matter.

The extent of this work sets an antecedent for a better understanding and more accurate description of the general conditions of the study area by using a coupled model of these characteristics, but more research needs to be conducted. The horizontal resolution needs to be increased to better describe the coastal regions where biomass production is higher during the spring and summer seasons. Among the future tasks to expand the understanding of the biological dynamics of the region and improve the accuracy of the results are the addition of all the rivers and the study of their influence in the area; the setup of a two-way grid refinement between the WRF and ROMS core models for allowing a higher definition of the coastal circulation patterns and their influence in Harima Nada; and the coupling of waves and sediment modules on COAWST. We are actually working toward the accurate estimation of the discharged freshwater and nutrient loads of the other rivers in the region by applying the same physically distributed model we had applied to the Kako River to get the offline input for new simulations in the area at higher horizontal resolution.

This work is the first step in the development of a biogeochemical model for the region that helps in assessing options for reverting the actual oligotrophication issues and that can forecast particular conditions under future scenarios.