1 Introduction

Marine environmental risk assessment is vital for ensuring the safety of sea creatures, which are threatened by the presence of persistent organic pollutants (POPs) in seawater and the consequent migration of these pollutants into the bodies of marine life forms [1]. The negative effects of POPs on marine life forms include reproductive abnormalities, deformities, and immune and neurological effects. These negative effects could also impact humans who eat marine life forms [2]. Rivers constitute the major pathways for the dissemination of POPs in coastal areas [3]. Examining the dynamics of effluents from rivers may clarify the modes of transport of POPs, providing valuable input into methods for marine environmental risk assessment.

The region of freshwater influence (ROFI) [4] has been of academic interest because of the diverse interactions between coastal flows, productive ecosystems, sediment transport, and the movement of substances from rivers. Many studies have addressed the ROFI from the perspective of physical oceanography, including phenomena such as low-density water mechanics (e.g., [5,6,7]). The majority of those studies have investigated freshwater plumes from rivers with freshwater discharge rates greater than 1000 m3/s. In contrast, few studies (e.g., [8, 9]) have been conducted on small rivers with freshwater discharge rates of 100 m3/s or less.

Recent studies have revealed the significance of small rivers as contributors to POP loads in coastal regions [10, 11]. The concentrations of microplastics in rivers that absorb POPs [12] are approximately proportional to population density and urban ratio rather than to the basin area [10], with microplastic concentrations in small rivers in Japan being higher than those in large rivers [11]. The transport process of microplastics is similar to that of other pollutants [10]. For example, microplastics have been found along the coastal area of Sagami Bay, suggesting that they were transported by small rivers, such as the Sakawa River, which has an average flow rate of 20 m3/s [13]. Given the magnitude of their potential contribution to seawater pollution, it is crucial to assess the impact of pollutants loaded from small rivers.

The collection of in situ biological data constitutes the most popular method for environmental risk assessment because it offers real information on pollution in the bodies of marine life forms. However, the available data are spatially and temporally limited.

Coastal area utilization (such as large-scale aquaculture) is expected to increase in the future. To manage the increases in coastal area utilization while acknowledging the environmental risk imposed by pollutants, sustainability should be ensured in a manner that accommodates the increasingly utilized coastal areas while remaining sensitive to the presence of risk. Cyberspace is expected to be a key resource for evaluating the sustainability of human activity, as proposed in Society 5.0, a Japanese science and technology policy [14] in which ocean modeling has the potential to play a role. Although previous ocean modeling studies have succeeded in simulating and elucidating material cycles in the ocean [15,16,17], to enable to use the ocean model as a tool for assessing environmental risk, further refinement of its physical aspects is required. Specifically, it is essential to be able to simulate the dynamics of pollutants from small rivers.

River plume length may serve as a practical index for determining the extent to which POPs spread. Therefore, the purpose of this study was to present a method for estimating river plume length using easily obtainable information. To achieve this, a numerical model was developed to simulate the dynamics of the river plume around the mouth of the Sakawa River—a typical small river flowing through an urban area—near Odawara, Sagami Bay, Japan. The proposed model employed numerical computation methods such as the constrained interpolation profile scheme, fractional time-step method, and semi-Lagrangian scheme. Well-known models such as the Princeton ocean model (POM) [18] and the unstructured grid finite volume coastal ocean model (FVCOM) [19] are available and very useful for simulating coastal ocean physics. In this study, an original program was implemented for performing numerical computations because this study has aimed to develop a marine environmental risk assessment method that incorporates physical, biological, and chemical aspects of the marine environment and risk takers’ behaviors. The present paper reports the physical aspects of the simulation results produced by our program.

The remainder of this paper is organized as follows. Section 2 describes the methods used for the numerical simulations, observations, and model validation. The results are presented and discussed in Sects. 3 and 4, respectively. Finally, the conclusions are presented in Sect. 5.

2 Numerical simulation model and observations

2.1 Numerical simulation model

In this subsection, the governing equations and algorithms, boundary conditions, and topographic and forcing data are described.

2.1.1 Governing equations and algorithms

The dynamics of seawater flows on an f-plane were governed by Navier–Stokes equations with Boussinesq and hydrostatic approximations, expressed as:

$$\begin{array}{c}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{{\rho }_{0}}\frac{\partial p}{\partial x}+fv+\frac{\partial }{\partial x}\left({K}_{\text{H}}\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left({K}_{\text{H}}\frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial z}\left({K}_{\text{V}}\frac{\partial u}{\partial z}\right),\end{array}$$
(1)
$$\begin{array}{c}\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{{\rho }_{0}}\frac{\partial p}{\partial y}-fu+\frac{\partial }{\partial x}\left({K}_{\text{H}}\frac{\partial v}{\partial x}\right)+\frac{\partial }{\partial y}\left({K}_{\text{H}}\frac{\partial v}{\partial y}\right)+\frac{\partial }{\partial z}\left({K}_{\text{V}}\frac{\partial v}{\partial z}\right),\end{array}$$
(2)
$$\begin{array}{c}0=-g-\frac{1}{\rho }\frac{\partial p}{\partial z},\end{array}$$
(3)

where \(\left(x,y,z\right)\) denote the eastward, northward, and vertically upward coordinate axes, respectively; \(\left(u,v,w\right)\) denote the flow velocities in the \(x\)-, \(y\)-, and \(z\)-directions, respectively; and \(\left({K}_{\text{H}},{K}_{\text{V}}\right)\) denote the horizontal and vertical viscosity coefficients, respectively. The symbol \({\rho }_{0}\) represents the reference density of the seawater, and \(f\) represents the Coriolis parameter at a latitude of 35° 15ʹ N. The quantity \(g\) is the gravitational acceleration, and \(p\) is the pressure. The dependence of seawater density \(\rho\) on water temperature T, salinity \(S\), and pressure \(p\) is expressed as follows:

$$\begin{array}{c}\rho =\rho \left(T,S,p\right),\end{array}$$
(4)

which was computed by applying the international equation of state of seawater [20].

The continuity of the incompressible fluid was ensured by solving the following equation:

$$\begin{array}{c}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0.\end{array}$$
(5)

The transports of heat and salt in seawater were governed by the following advection–diffusion equations, respectively:

$$\begin{array}{c}\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}=\frac{\partial }{\partial x}\left({D}_{\text{H}}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({D}_{\text{H}}\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left({D}_{\text{V}}\frac{\partial T}{\partial z}\right),\end{array}$$
(6)
$$\begin{array}{c}\frac{\partial S}{\partial t}+u\frac{\partial S}{\partial x}+v\frac{\partial S}{\partial y}+w\frac{\partial S}{\partial z}=\frac{\partial }{\partial x}\left({D}_{\text{H}}\frac{\partial S}{\partial x}\right)+\frac{\partial }{\partial y}\left({D}_{\text{H}}\frac{\partial S}{\partial y}\right)+\frac{\partial }{\partial z}\left({D}_{\text{V}}\frac{\partial S}{\partial z}\right)+\left({\text{source}}\right),\end{array}$$
(7)

where \(\left({D}_{\text{H}},{D}_{\text{V}}\right)\) denote the horizontal and vertical diffusion coefficients, respectively [21]. The term "\(\left({\text{source}}\right)\)" in Eq. 6 represents the salinity variations caused by the freshwater inflow [22].

To perform efficient numerical computations, the abovementioned equations were separated into two systems: one system accommodated rapidly varying phenomena (external mode) and the other accommodated slowly varying phenomena (internal mode). A semi-Lagrangian scheme was applied for the computation of the external mode. For the internal mode, a fractional time-step approach was employed to achieve efficient and accurate computations for each equation term (see Fig. A1-1 for details). The time steps of the external and internal modes were set to 0.1 s and 0.4 s, respectively. These modes interacted by exchanging variables at regular intervals such that the two separated systems of equations were theoretically equivalent to the primitive systems of equations (Eqs. 1–7) (see [22] and [23] for further details).

2.1.2 Boundary conditions

The horizontal flow components normal to the coastline were set to zero, serving as the kinematic boundary conditions in the horizontal direction. The kinematic boundary conditions in the vertical direction were as follows:

$$\begin{array}{c}\frac{\partial \eta }{\partial t}+u\frac{\partial \eta }{\partial x}+v\frac{\partial \eta }{\partial y}-w=0,\end{array}$$
(8)

on the free surface, and

$$\begin{array}{c}u\frac{\partial H}{\partial x}+v\frac{\partial H}{\partial y}-w=0,\end{array}$$
(9)

on the bottom surface, where \(\eta\) denotes the sea surface elevation from the resting free surface. \(H\) denotes the distance from the sea bottom to the resting free surface.

The boundary conditions for the horizontal viscosity terms in Eqs. 1 and 2 had the form of no-slip conditions and were written as follows:

$$\begin{array}{c}\left\{\begin{array}{l}{K}_{\text{V}}\frac{\partial u}{\partial z}=\frac{{\tau }_{x}}{{\rho }_{0}}\\ {K}_{\text{V}}\frac{\partial v}{\partial z}=\frac{{\uptau }_{y}}{{\rho }_{0}}\end{array}\right.,\end{array}$$
(10)

where \(\left({\tau }_{x},{\tau }_{y}\right)\) denote the stresses in the \(x\)- and \(y\)-directions, respectively, and were determined by

$$\begin{array}{c}\left\{\begin{array}{l}{\tau }_{x}={C}_{\text{d}}\hspace{0.17em}{\rho }_{\text{a}}\hspace{0.17em}U\sqrt{{U}^{2}+{V}^{2}}\\ {\tau }_{y}={C}_{\text{d}}\hspace{0.17em}{\rho }_{\text{a}}\hspace{0.17em}V\sqrt{{U}^{2}+{V}^{2}}\end{array}\right.,\end{array}$$
(11)

on the sea surface \(\left(z=0\right)\) and by

$$\begin{array}{c}\left\{\begin{array}{l}{\uptau }_{x}={C}_{\text{f}} {\rho }_{0}\hspace{0.17em}u\sqrt{{u}^{2}+{v}^{2}}\\ {\uptau }_{y}={C}_{\text{f}}\hspace{0.17em}{\rho }_{0}\hspace{0.17em}v\sqrt{{u}^{2}+{v}^{2}}\end{array}\right.,\end{array}$$
(12)

on the sea bottom \(\left(z=-H\right)\), where \({C}_{\text{d}}\) is the bulk coefficient of the wind stress, and \({C}_{\text{f}}\) is the friction coefficient of the sea bottom. The quantity \({\rho }_{\text{a}}\) is the density of air. \(\left(U,V\right)\) denote the wind velocities in the \(x\)- and \(y\)-directions, respectively.

The boundary conditions for the vertical diffusion terms in Eqs. 6 and 7 were:

$$\begin{array}{c}\left\{\begin{array}{l}{D}_{\text{V}}\frac{\partial T}{\partial z}=\frac{Q}{{C}_{\text{p}}^{\text{wtr}}\cdot {\rho }_{0}}\\ {D}_{\text{V}}\frac{\partial S}{\partial z}=0 \end{array}\right.,\end{array}$$
(13)

on the sea surface \(\left(z=0\right)\) and

$$\begin{array}{c}\left\{\begin{array}{l}{D}_{\text{V}}\frac{\partial T}{\partial z}=0\\ {D}_{\text{V}}\frac{\partial S}{\partial z}=0\end{array}\right.,\end{array}$$
(14)

on the sea bottom \(\left(z=-H\right)\), where \(Q\) denotes the heat flux comprising shortwave radiation, longwave radiation, sensible heat, and latent heat (see Appendix 2 for details). \({C}_{\text{p}}^{\text{wtr}}\) denotes the specific heat of the seawater. In Eq. (13), it was assumed that there was no salinity variation due to precipitation and evaporation.

2.1.3 Open boundary conditions

The model domain was bounded by the coastline of Sagami Bay and the eastern and southern open boundaries; these boundaries had to be put for simulating the ocean dynamics in the localized area (Fig. 1). An open-boundary condition that allowed the smooth propagation of waves without reflections there was imposed by employing no-reflection scheme [24, 25] in conjunction with the semi-Lagrangian scheme.

Fig. 1
figure 1

Map of Sagami Bay (left) and grid map of the study area (right). The black-outlined circle and white-bordered square indicate the mouth of the Sakawa River in the two images, respectively

Full size image

The boundary condition imposed on the southern open boundary was expressed as [23],

$$\begin{array}{c}{V}^{+}=\frac{1}{2}f{u}{\prime}\Delta t+\frac{g}{c}\left({h}_{0}+2{\eta }_{{\text{I}}}\right),\end{array}$$
(15)

where \(\Delta t\) is a discrete time step, \(c\) is the phase speed of the long wave excluding the effect of self-rotation of the Earth, \({h}_{0}\) is the water column height at calm state, and \({\eta }_{{\text{I}}}\) is the displacement of the water surface by an incident wave propagating from the outer region. \({V}^{+}\) and \({u}{\prime}\) are functions of several variables, the details of which are provided in Appendix 3. \({\eta }_{{\text{I}}}\) was calculated by summing the astronomical tides and sea surface elevations due to low atmospheric pressure. The former was superimpositions of the eight major tidal constituents (M2, S2, K1, O1, P1, N2, K2, and Sa) determined using the tidal harmonic constants at Odawara, which were available from the Japan Meteorological Agency (JMA) [26]. The latter was supposed to be in proportion to the difference in sea level pressure between reference and measured values [27] (see Sect. 2.1.4 and Appendix 4 for details).

The boundary condition imposed on the eastern open boundary was expressed as [23],

$$\begin{array}{c}{U}^{-}=\frac{1}{2}f{v}{\prime}\Delta t+\frac{g}{c}{h}_{0},\end{array}$$
(16)

where \({U}^{-}\) and \({v}{\prime}\) are functions of several variables (see Appendix 3 for details).

In a default calculation, this boundary condition (Eq. 16) was imposed on the whole of the eastern open boundary. Another calculation was also conducted to consider the counterclockwise circulation in Sagami Bay (e.g., [28]): in this case, the fluid was forced to inflow through the northern half of the eastern open boundary, which was expressed as,

$$\begin{array}{c}{U}^{-}={U}^{+}+f{v}^{\prime}\Delta t-2{u}_{{\text{in}}},\end{array}$$
(17)

where \({U}^{+}\) is a function of several variables (see Appendix 3 for details). \({u}_{{\text{in}}}\) denotes the inflow velocity and was set to 0.15 m/s, as indicated in [28], and the no-reflection condition (Eq. 16) was applied to the southern half of the eastern open boundary.

2.1.4 Topographic and forcing data

The simulation period consisted of three spans of 15 days each, starting on June 2, June 21, and July 7, 2022, respectively. These periods included the three observation periods noted later. The simulation domain (35° 12ʹ 54ʹʹ N–35° 16ʹ 48ʹʹ N, 139° 8ʹ 6ʹʹ E–139° 12ʹ 54ʹʹ E) was discretized by spanning rectangular grids. The horizontal grid scales were set to 250 m in both the \(x\)- and \(y\)-directions. Digital bathymetric data with isobaths at 10 m intervals [29] were interpolated using Akima’s method [30] to set the depths and sea bottom gradients for each horizontal grid. The vertical grid scales were set to 0.10 m on the surface layer and increased with the vertical positions (Fig. 2). The aforementioned grid resolutions were determined as a compromise between avoiding a non-viable computational load and ensuring that the surface thermohaline dynamics were captured with a satisfactory resolution. The initial values of salinity and water temperature were set for each of the thirty-one vertical layers based on the data in [31].

Fig. 2
figure 2

Vertical grid scales. The intervals are 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 20, 50, 100, and 160 m in the 0–1 m, 1–2 m, 2–4 m, 4–6 m, 6–10 m, 10–15 m, 15–35 m, 35–85 m, 85–185 m, and 185–825 m layers, respectively

Full size image

In the simulation, the dynamics was driven by wind stresses, heat fluxes, and riverine water discharge rates. The wind velocities measured at Odawara provided by the JMA [27] (see Fig. A4-1 for details) were used to determine the wind stresses, which were assumed to be applied uniformly on the entire study area. Heat fluxes were computed using meteorological data (air temperature, relative humidity, wind velocity, and hourly sunshine rate) obtained at Odawara, and weather and sea level pressure data obtained at the cities of Yokohama and Mishima, as provided by the JMA [27] (see Appendix 4 for details). These meteorological data with 1-h intervals were linearly interpolated and input into the program. The computed components of the solar radiation included in \(Q\) were validated by comparisons with radiation levels observed at Motomaki (35° 24ʹ 52ʹʹ N, 139° 39ʹ 42ʹʹ E) by the City of Yokohama [32] (Fig. 3).

Fig. 3
figure 3

Temporal variations in solar radiation on the ground produced by the numerical computation (red) and measured at the Motomaki site (Naka ward) at 35° 24ʹ 52ʹʹ N, 139° 39ʹ 42ʹʹ E by the City of Yokohama [32] (blue)

Full size image

The discharge rate of the Sakawa River was assumed to be the sum of the downstream discharge rate from the Miho Dam [33] and the mean discharge rate of 15 m3/s from other tributaries [31] (Fig. 4).

Fig. 4
figure 4

Time series of the river discharge rate of the Sakawa River (solid line) estimated by linear interpolations of the sum of 15 m3/s from other tributaries [31] (dashed line) and the downstream discharge rates from the Miho Dam every hour [33] (circle plots)

Full size image

2.2 In situ observations

In situ observations were conducted on June 9, June 28, and July 14, 2022, to obtain salinity and temperature distributions in seawater for validating the numerical simulation of the thermohaline fields in the surface layers of the ocean. These dates coincided with a neap tide, spring tide, and another spring tide, respectively. The weather conditions on the 3 days were: cloudy with some sunshine, sunny, and cloudy, respectively.

Measurements of seawater conditions were performed along four or five observation lines perpendicular to the coastline at intervals of 1.0–1.5 km (lines A to N in Fig. 5). Individual measurement sites were placed on the observation lines at intervals of 250–500 m (Fig. 5). The total number of measurement sites encompassing all 3 days of observation was eighty-three (twenty-eight sites on June 9, thirty sites on June 28, and twenty-five sites on July 14).

Fig. 5
figure 5

The measurement sites (shaded circles and triangles) and observation lines (solid lines indicated by letters of the alphabet) on June 9 (left), June 28 (middle), and July 14 (right). The profiler could be immersed to a depth of 30 m at the sites indicated by shaded circles. At the sites indicated by shaded triangles, immersion of the profiler to a depth of 30 m was not possible. The large unshaded circle indicated the mouth of the Sakawa River. The sites along each lettered line were numbered from 1 onward, starting with the site nearest to the coastline. Each site was represented by the assigned line letter and its number (for example, the sites on line J were represented as J1, J2, J3, J4, and J5)

Full size image

The research vessel Tachibana (Manazuru Marine Center for Environmental Research and Education, Yokohama National University) was operated along the respective observation lines. The measurement durations (starting at the first site, A1, E1 and J1, and finishing at the last one, D1, I6, and N5) were 4 h 33 min, 4 h 26 min, and 3 h 49 min on June 9, June 28, and July 14, respectively. Vertical profiles of salinity and water temperature in a surface layer 30 m deep were acquired with a resolution of 0.2 m using a conductivity–temperature–depth profiler (CTD, Ocean Seven 316, Serial number 1202287, Idronaut, Brugherio, Italy). Errors in the measurement of depths, which might be caused mainly by time variations in atmospheric pressure, were estimated to be within 5 cm. This error estimation was based on the recorded time variations in atmospheric pressure provided by JMA (Fig. A4-6). At certain sites (as indicated in Fig. 5), the profiler could not descend to a depth of 30 m because the sea was too shallow there.

2.3 Error metrics for model validation

This study used two error metrics to identify whether the model errors applied to an entire vertical profile or only to a partial extent of the profile; mean bias (MB) and centered root mean square error (CRMSE) [34]. The root mean square error (RMSE) was defined as

$$\begin{array}{c}{\text{R}}{\text{M}}{\text{S}}{\text{E}}=\sqrt{\frac{1}{N}{\sum }_{i=1}^{N}{\left({m}_{i}-{o}_{i}\right)}^{2}},\end{array}$$
(18)

where \({m}_{i}\) and \({o}_{i}\) are the simulated and observed values, respectively, at each measurement site (Fig. 5), \(i\) is the sequential number of the vertical positions at a particular observation site where salinity and water temperature were profiled, and \(N\) denotes the total number of vertical measurement positions \(i\) at that site. The simulated results for neighboring model grids were linearly interpolated to obtain \({m}_{i}\).

The MB, which quantified the overall error between the simulated and observed values, was defined as

$$\begin{array}{c}{\text{M}}{\text{B}}=\overline{m }-\overline{o },\end{array}$$
(19)

where \(\overline{m }\left(=1/N{\sum }_{i=1}^{N}{m}_{i}\right)\) and \(\overline{o }\left(=1/N{\sum }_{i=1}^{N}{o}_{i}\right)\) are the means of the simulated and observed values at each site, respectively. The CRMSE was defined as

$$\begin{array}{c}{\text{C}}{\text{R}}{\text{M}}{\text{S}}{\text{E}}=\sqrt{\frac{1}{N}{\sum }_{i=1}^{N}{\left({m}_{i}{\prime}-{o}_{i}{\prime}\right)}^{2}},\end{array}$$
(20)

where \({m}_{i}{\prime}\left(={m}_{i}-\overline{m }\right)\) and \({o}_{i}{\prime}\left(={o}_{i}-\overline{o }\right)\) are the deviations from the means of the simulated and observed values, respectively. It quantified the difference between these two deviations from the mean at each site. Calculating these two error metrics enabled us to validate the simulations of the thermohaline fields in detail.

The error metrics were calculated for the data obtained for the 30-m deep surface layer at each of the seventy-two sites where the profiler could be immersed to a depth of 30 m. The metrics were calculated for the data obtained within the 10-m deep surface layer at all sites.

3 Results

3.1 Water levels

The numerical model produced a seawater level variation modulated by a semidiurnal tidal component (Fig. 6), exhibiting a minimum twice a day. In the following description, the lower of the two low tides in a day is referred to as the lower low tide (LLT) and the higher of the two low tides is referred to as the higher low tide (HLT). The river plumes at different LLTs are examined in Sect. 3.3.

Fig. 6
figure 6

Temporal variations in water levels produced by the numerical computation (red) and measured at Odawara by the JMA [35] (blue). The acronyms LLT and HLT in the upper panel indicated a pair of lower low tides and higher low tides, respectively

Full size image

The tidal range decreased during the period from June 2 to June 9, reaching a minimum of 0.6 m on June 9. Subsequently, the range increased to 1.4 m on June 15. During the other two 15-day spans (June 21 to July 5 and July 7 to July 21), the tidal ranges also increased and decreased over a period of approximately 7 days. Water levels at the LLTs were maximally 1.0 m lower than those at the HLTs during spring tides. The differences in the water levels at the two high tides in a day were within 0.25 m on all days. The average absolute error between the computed and observed [35] water levels was 0.0141 m, and the correlation coefficient between the computed and observed water levels was 0.986.

3.2 Vertical profiles and error metrics

Decreases in salinity within the surface 1.0-m layer at sites B2, G2, and L1 were reproduced by numerical simulation (Fig. 7). At sites where the observed salinities decreased to less than 25 psu (e.g., B2 and H4), the decrease in salinity computed by the simulation was insufficient. The CRMSEs in the surface 10-m layer at these sites exceeded 1.5 psu (Table 1) and the percentage of sites with CRMSEs exceeding 1.5 psu was 6% (Fig. 8). The computed salinities in the 10–30 m layer were in the range of 32–33 psu, which are very close to the observed salinities at all sites (Fig. 7 and Fig. A5-1). The MBs of the salinities in the surface 30-m layer were within \(\pm\) 0.5 psu at all seventy-two sites (Fig. 8).

Fig. 7
figure 7

Vertical profiles of computed (red) and observed (blue) salinity

Full size image
Table 1 Mean bias (MB) and centered root mean square error (CRMSE) of salinity at six selected sites
Full size table
Fig. 8
figure 8

Scatter plots of mean bias (MB) and centered root mean square error (CRMSE) of salinity in the surface a 30-m and b 10-m layers. Circle, triangle, and cross plots represent the results for the data on June 9, June 28, and July 14, respectively

Full size image

The increasing water temperature toward the surface was captured to a satisfactory degree by the numerical simulation (Fig. 9 and Fig. A5-2). Some of the observed vertical profiles of the water temperatures displayed stepwise patterns (e.g., G2 and H4), which involved several layers, each approximately 5 m in thickness. These were not reproduced by the numerical simulation. The CRMSEs in the surface 30-m layer at most sites on June 28, which included stepwise patterns, exceeded 1.5 ℃ (Table 2 and Fig. 10).

Fig. 9
figure 9

Vertical profiles of computed (red) and observed (blue) temperature

Full size image
Table 2 Mean bias (MB) and centered root mean square error (CRMSE) of temperature at six selected sites
Full size table
Fig. 10
figure 10

Scatter plots of mean bias (MB) and centered root mean square error (CRMSE) of temperature in the surface a 30-m and b 10-m layers. Circle, triangle, and cross plots represent the results for the data on June 9, June 28, and July 14, respectively

Full size image

3.3 Sea surface salinity maps

The characteristics of freshwater transport from the river were identified using computed surface salinity maps (Figs. 11 and 12). Freshwater traveled from the mouth of the river along the shore toward both the southwest and northeast. The directions of the river plume motion correlated well with water levels: at LLTs, the freshwater was predominantly transported to the southwest, whereas at HLTs and high tides, it was predominantly transported to the northeast.

Fig. 11
figure 11

Simulated surface salinity maps during the spring tide on July 14 at a high tide (3:45), b lower low tide (11:00), c high tide (18:30), and d higher low tide (23:30). The 31.0 psu contour is represented by the outermost line. The contour interval is 1.0 psu

Full size image
Fig. 12
figure 12

Simulated surface salinity maps at lower low tides during the a middle, b spring, c middle, and d neap tide on July 11, 14, 17, and 20, respectively. The 31.0 psu contour is represented by the outermost line. The contour interval is 1.0 psu

Full size image

The salinity contour lines corresponding to 30 and 27 psu traveled 3.5 km and 2.4 km, respectively, toward the southwest at LLT on July 14 (Fig. 11). At HLT and high tides, the distances from the river mouth to the farthest points on these contour lines were less than 2.2 km and 0.8 km, respectively. During the neap tide, the distances of the river plume’s edge from the river mouth at LLT were 0.5–2.5 km (Fig. 12), while during the spring tide, these distances were approximately 5.0 km or less.

To determine the extent to which substances from the river were transported and what factors determined that extent, the river plume lengths at LLTs (the times when freshwater from the river was transported the farthest) were analyzed in detail. Hereafter, river plume length, denoted by \({l}_{\text{p}}\), is defined as the distance from the river mouth to the farthest extent of the 30 psu salinity contour line.

River plume length was negatively and positively correlated with the water level (Fig. 13) and river discharge rate (Fig. 14), respectively. The relationship between river plume length and water level from July 8 to July 16, when the river discharge rates were almost constant, was examined. The river plume length correlated more negatively with the water level from July 8 to July 16 than during any of the other periods (Fig. 15).

Fig. 13
figure 13

Scatter plot of river plume length versus water level. The correlation coefficient is – 0.60

Full size image
Fig. 14
figure 14

Scatter plot of river plume length versus river discharge rate. The correlation coefficient is 0.044

Full size image
Fig. 15
figure 15

Scatter plot of river plume length \({l}_{\text{p}} \left[{\text{km}}\right]\) versus water level \(\upeta \left[{\text{m}}\right]\) from July 8 to July 16. The correlation coefficient is – 0.970. The equation of the regression line is \({l}_{\text{p}} \left[{\text{km}}\right]=-2.90\upeta \left[{\text{m}}\right]+0.872\). The coefficient of determination is 0.94, and the p value is less than 0.001

Full size image

A regression line was used to estimate the river plume length from the water level. The intercept of the regression line depended on the river discharge rate. Hereafter, the intercept is referred to as reference river plume length, denoted by \({l}_{0}\). Considering the time elapsed since the riverine water discharge, \({l}_{0}\) was presumed to be closely related to the river discharge rate earlier than the time at which \({l}_{0}\) was determined. The strongest correlation was detected for a river discharge rate 25 h earlier, as denoted by \(R\) (Figs. 16 and 17).

Fig. 16
figure 16

Scatter plot of the reference river plume length \({l}_{0} \left[{\text{km}}\right]\) versus the river discharge rates 25 h earlier than the time at which the river plume lengths were determined \(R \left[{\text{m}}^{3}\text{/}{\text{s}}\right]\). The equation of the regression line is \({l}_{0} \left[{\text{km}}\right]=0.143R \left[{\text{m}}^{3}\text{/}{\text{s}}\right]-0.935\). The coefficient of determination is 0.44, and the p value is less than 0.001

Full size image
Fig. 17
figure 17

Correlation coefficients between the reference river plume length \({l}_{0}\) and river discharge rates. The definition of time lag is the time at which the river plume length was determined minus that at which the river discharge rate was determined

Full size image

A multiple regression relation for the river plume length as a function of both the water level and river discharge rate was constructed, utilizing the regression lines described in the captions of Figs. 15 and 16. This relation was used to obtain an estimated river plume length \({l}_{\text{p}}^{\text{est}}\) from the preceding river discharge rate. The values of \({l}_{\text{p}}^{\text{est}}\) at LLTs were calculated for 39 days (Fig. 18). Comparisons of \({l}_{\text{p}}^{\text{est}}\) with river plume lengths determined from the numerically simulated salinity maps demonstrated that the multiple regression equation accurately generated the river plume lengths. The maximum absolute error was 1.25 km on July 18. The estimations exhibited errors of less than 0.5 km on 22 days, and errors of less than 1.0 km on 35 days (Fig. 18).

Fig. 18
figure 18

River plume lengths determined from the salinity maps generated by numerical calculation (shaded circles), those generated by numerical calculation considering the counterclockwise circulation (shaded triangles), and from the regression line results calculated by \({l}_{\text{p}}^{\text{est}} [{\text{km}}]=-2.90\upeta \left[{\text{m}}\right]+0.143R \left[{\text{m}}^{3}\text{/}{\text{s}}\right]-0.935\) (unshaded circles)

Full size image

3.4 River plume length considering the counterclockwise circulation

When imposing the inflow from the eastern open boundary to consider the counterclockwise circulation in Sagami Bay, the river plume length was on average 0.60 km shorter than those without it (Fig. 18). Compared to the river plume lengths estimated by the regression analysis, the average and maximum absolute error were 1.23 km and 3.10 km on July 17, respectively.

The river plume length at LLT on July 14 (Fig. 19) was 6.0 km toward the southwest, which is longer compared to the river plume length for the case without the inflow (Fig. 11). The river plume lengths at most LLTs (e.g., Fig. 20a and d) were shorter than the river plume lengths for the case without the inflow (Fig. 18). When comparing salinity contour lines, the 30 psu line with the counterclockwise circulation considered is shorter, while the 31 psu line is longer.

Fig. 19
figure 19

Simulated surface salinity maps considering the counterclockwise circulation during the spring tide on July 14 at a high tide (3:45), b lower low tide (11:00), c high tide (18:30), and d higher low tide (23:30). The 31.0 psu contour is represented by the outermost line. The contour interval is 1.0 psu

Full size image
Fig. 20
figure 20

Simulated surface salinity maps considering the counterclockwise circulation at lower low tides during the a middle, b spring, c middle, and d neap tide on July 11, 14, 17, and 20, respectively. The 31.0 psu contour is represented by the outermost line. The contour interval is 1.0 psu

Full size image

4 Discussion

The time series of observed water levels (Fig. 6) showed spring and neap tides, as well as the daily occurrence of high and low tides approximately twice a day. These temporal variations in water level were also showed in the numerically computed ones, indicating the practical utility of this model.

The decreases in salinity within the surface 1.0-m layer induced by river plume intrusion were captured moderately well by the model, in which the low-salinity water motion was driven primarily by the horizontal advection term in Eq. 7. However, the decrease in salinity at sites B2 and H4 computed by the model was too small (Fig. 7). These discrepancies might have occurred because of the treatment of salinities at the open boundaries, where the salinities were assumed to be constant. If the model covers a larger area, the impact of the open-boundary process can be reduced. Achieving a wide range of coverage would require a reduction in the computational load.

The increase in water temperature toward the surface was well simulated by the model, demonstrating adequate computation of the heat budget in the surface layer. However, the distribution of multiple thermoclines, that is, stepwise vertical temperature profiles, could not be simulated successfully. The error metrics were small for the temperatures measured on June 9 and July 14. Considering that it was cloudy during those 2 days, the seawater in the surface layer was heated for only a few hours or less. The large error metrics of the temperature on June 28 resulted from an insufficient simulation of the observed stepwise temperature profiles. On this day and a few days earlier, the actual seawater had been continuously heated in the daytime. Over these heating periods, the vertical mixing of seawater might has been suppressed under real conditions. The weakening of the mixing and the resulting formation of thermoclines could not be simulated. This problem may be resolved by improving the method for computing the vertical diffusion coefficients and using finer vertical grid resolutions in layers deeper than 10 m.

The river plume length was correlated with the water level. This correlation may be ascribed to the cyclic variations in the tidal flows. The river plume length was most closely correlated with the river discharge rate 25 h earlier in this case. Freshwater transport from the river mouth occurred mainly around LLTs (Fig. 11b). The freshwater discharged after an LLT remained near the mouth of the river (Fig. 11a, c, and d) and traveled several kilometers away around the succeeding LLT.

The difference between the river plume lengths calculated through the numerical simulations with and without the counterclockwise circulation may arise from the horizontal mixing of the low=salinity water mass with the surrounding offshore water mass induced by the inflow. This inflow conveyed the leading edges of river plume by longer distances along the coast (Fig. 19), and the resulting water mixing made the river plume to have higher salinities (Fig. 20).

A newly introduced explanatory variable (e.g., inflow velocity) may be hopeful for constructing a regression line that yields more accurate estimations of river plume dispersion. In this simulation, the inflow boundary condition was simply treated. A more advanced treatment of the inflow may be required to better simulate the flow pattern in the local area around the river mouth.

This study addressed the river plume dynamics in the summer season under weak wind condition. In the winter season, strong winds occur more frequently. Estimations of river plume length under strong wind condition may necessitate an inclusion of wind as an explanatory variable in addition to water levels and river discharge rate [36].

5 Conclusion

In this study, a numerical simulation model was developed to elucidate the dynamics of effluents from the Sakawa River in Sagami Bay and to realize a computer-aided marine environmental risk assessment. The model applied fractional time-step approach to elicit the advantages of numerical computation schemes for the respective terms in the governing equations. A fine vertical resolution was implemented to capture the dynamics of thermohaline fields in the surface layer. In situ observations were carried out to measure salinity and temperature at vertical intervals of 0.2 m in the surface 30-m layer, at eighty-three sites off the mouth of the Sakawa River. The results produced by the model were validated by comparison with observed data. Decreases in salinities within the surface 1.0-m layer were reproduced by numerical simulation. Regression analyses using the simulation results revealed that the river plume length correlated more strongly with the river discharge rate 25 h earlier than with the simultaneous river discharge rate. Applying regression analyses using such a correlation, the Sakawa river plume length under weak wind condition was well estimated. This result can contribute to the realization of a method for estimating river plume length using easily accessible data such as water level and river discharge rate.