Modeling runoff into a Region Of Freshwater Influence for improved ocean prediction : Application to Funka Bay

Regions Of Freshwater Influence (ROFI) existing between oceans and estuaries often yield significant fishery resources. To inform fishermen of the real-time ocean state, an operational prediction system focusing on a specific coastal region would normally require accurate and timely runoff data for all rivers. In this paper, hydrometeorological procedures providing the required runoff data on a daily basis were coupled with an Ocean General Circulation Model (OGCM) to evaluate the impacts of runoff processes on ocean simulations within a ROFI. The procedures we adopted employ a distributed tank model based on water mass and heat budgets derived from predicted meteorological datasets. An exponential relationship between runoff rates and watershed areas was used to determine model parameters in order to estimate the runoff from many other small rivers. The coupled model reproduced a surface salinity field in the bay that was in good agreement with observations, and simulated the expected clockwise circulation generated by the high net total discharge associated with snowmelt. Our results underline the fact that implementation of hydrological processes into ocean simulations is essential for a better understanding of water circulation driven by runoff into semi-enclosed bays over interannual timescales.


INTRODUCTION
The Regions Of Freshwater Influence (ROFI), conceptualized simply as a coastal basin into which numerous rivers discharge, occupies an oceanographic setting between the outer sea and the estuary (Simpson, 1997).Freshwater runoff from rivers provides buoyancy and momentum input together with coastal currents, and generates a complex physical regime and circulation system with properties that are entirely different to those of the outer ocean.Owing to the deposition and re-suspension of land-derived materials from riverine sources, the low salinity water spreading from the estuary should retain high nutrient concentrations that may often be associated with eutrophication.Therefore, an accurate determination of riverine influences on buoyancy and nutrients is essential for ROFI studies, and should lead to a better understanding of ocean circulation and ecosystems relevant to coastal fisheries, aquaculture, and marine harvesting.
Funka Bay, a typical ROFI in southwestern Hokkaido of northern Japan, exhibits a bowl-shaped bottom topography with a diameter of about 50 km that interfaces with the northwestern Pacific Ocean (Figure 1).More than 90% of the terrain surrounding Funka Bay is mountainous with a steep gradient (1/25) (Figure S1 in Supplementary Material, SM1) and is generally covered with beech forests (Fagus crenata) (Figure S2).The geology of the watersheds is largely characterized by volcanic Neogene or Paleogene sediments (andesite), implying a high permeability (Figure S2).The rivers flowing into Funka Bay have small watersheds (up to approximately 500 km 2 ).There are no runoff values for most of the watersheds and, while some data are derived, it is rarely enough to provide timely input data for OGCM simulations.
Rapid runoff carrying abundant nutrients and sediments can stay in the bay for over one month due to the weak water exchange (Miyake et al., 1988) and thereby greatly influence horizontal distributions of chlorophyll a (Figure S3) associated with fishery productivity.For example, scallops cultivated in coastal areas represent an important fishery resource (Shumway, 1991) that could be influenced by the freshwater runoff.Runoff provided by snowmelt forms nutrient-rich, low salinity water and its buoyancy leads to a summertime clockwise eddy (Figure S4) in the surface layer of Funka Bay (Satoh et al., 2003).This eddy influences the concentrations of micronutrients, such as iron and cadmium, associated with primary production (Kudo and Matsunaga, 1998).The clockwise eddy can also flush out the oxygen-deficient bottom water that is locally generated by an intermittent increase in the oxygen consumption rates in the bottom layer which can severely harm flatfish and shrimp production.
In this paper, hydrometeorological procedures capable of evaluating accurate, daily runoff for operational ocean simulation of the ROFI are proposed based on the heat and water balances within watersheds around Funka Bay and which make use of the Grid Point Value datasets of the Meso-Scale Model (GPV-MSM) provided by the Japan Meteorological Agency.These procedures are coupled with the OGCM to evaluate the impacts of implementing runoff processes into simulations of the Funka Bay, Northern Japan.

METHODOLOGY
The hydrometeorological procedures used to evaluate runoff were composed simply of two parts in order to minimize the calculation cost for ocean forecasts, though there are other more realistic runoff models or estimation procedures.First, water mass and heat budgets associated with radiation were estimated using meteorological data (Figure 2a), and then secondly, runoff was evaluated employing the distributed tank model (Figure 2b).These procedures are referred to as the Hydrometeorological and multi-Runoff Utility Model (HaRUM) for descriptive purposes.

Data sets
Re-analysis and predicted meteorological datasets GPV-MSM as inputs to the HaRUM were used for mean hourly air temperature, T a [°C], precipitation, P [mm], cloud cover, Cl, relative humidity, Rh, and wind speed, W [m s −1 ].Four time series of daily river runoff were provided by the Hokkaido Civil Engineering Office, covering the Yurappu, Oshamanbe, Nukibetsu, and Osaru rivers.These data were used for the validation of estimated runoff and to fix the parameters in the model.The analysis period for model calibration was from 2006 to 2007 because the runoff data covered 1996-2007, while the GPV-MSM data covered 2006-2011.The GPV/MSM data, with horizontal resolution of 0.125 × 0.1°, were interpolated onto watershed grids at an interval of 500 × 500 m using the spline method (Figure 1).The watersheds were defined using the drainage direction matrix derived by combining stream networks and topography data using Gaussian-interpolation.The datasets used for topography were the 50 m meshed Digital Elevation Model (Geospatial Information Authority of Japan: GSI), GTOPO30 (U.S. Geological Survey), JEGG500 (Japan Oceanographic Data Center), and ETOPO1 (National Oceanic and Atmospheric Administration).The stream networks were derived from the GSI and by manual tracing

Water and heat budgets
An analysis of water mass and heat budgets at the ground surface using the GPV/MSM data are based on Kondo and Watanabe (1992).The recharge, R, used as an input to each tank model is expressed by the water budget, R = P − E + M, where P, E, and M are precipitation, actual evapotranspiration from land surfaces, and snowmelt [mm], respectively.The net radiation flux density at the surface, where f is the albedo at the surface, S↓ is the solar insolation and L a ↓ and L g ↑ are the downward and upward longwave radiation fluxes from the atmosphere and ground surface, respectively.If snow covers the ground surface (S > 0) and the snow surface temperature T s < T a , snowmelt (M) can occurr based on R n and the heat balance (Kondo and Watanabe, 1992).Detailed calculation procedures are described in SM5.

Distributed tank-model
The distributed tank model (e.g., Huang et al., 2007) is employed here, which is an idealized "lumped" hydrological model (e.g.Kim et al., 2005) for analyzing the water balance in a river basin.Assuming homogeneous landuse, the distributed tank model consists of three serial tanks represented as n = 1 − 3 (Figure 2b).
The first tank represents the conceptual discharge system at the land surface.The upper (lower) side outlet of the first tank produces a preferential (overland) flow near the surface, when the water level S 1 exceeds the specified depth threshold, Z 1 (Z 2 ) in mm.The bottom outlet of the first tank represents infiltration into subsurface soil layers.The side outlets of the second and third tanks generate interflow and base-flow, respectively.The bottom outlets of the second and third tanks depict shallow groundwater infiltration and 'leakage' via deep percolation, respectively.A large amount of 'leakage' should ultimately generate submarine groundwater discharge because all watersheds face the bay.Evapotranspiration in the case of no water storage in the first tank would then come from the lower tanks.
The distributed tank model using a spatial resolution of 500 m was applied to each grid.Each tank outlet in the Nth (N ≤ 105) watershed (Figure S4) has an associated parameter to calculate outflow using P, E, and M at time intervals of Δt = 1 [hour] based on a linear proportionality: where, Q m (N) and L n (N) are outflows from the mth side and the bottom outlets of the nth tank, respectively; A m (N) and B n (N) are the scaling coefficients of outflows for the mth side and nth bottom outlet, respectively; S n (N) is the water level of the nth tank (positive in the downward direction).
Outflow from the mth outlet only occurs when S n (N) > Z m (N).Using Q m (N) derived by Equation (1) at time t [hour], the total outflow from the Nth watershed is calculated by Each standard tank model requires 13 parameters be determined for each grid point.However, this represents an impossible task owing to the large number of parameters that require calibration.To solve this problem, we assumed that parameters are expressed as a function F(A(N)) of the Nth watershed area under conditions of homogeneous landuse and geology.A(N) is area of the Nth watershed.After trying several functions, F(A(N)) was matched with an exponential form F(A(N)) = α exp(−βA(N)): This idea is based on general hydrological speculations, implying that smaller watersheds exhibit higher discharge rates.All the parameters and the F(A(N)) dependencies used in the tank model need to be manually calibrated by a trialand-error method using four observed runoff data values and two objective criteria: (1) to produce the highest correlations between observed runoff (Q obs ) and estimated (Q sim ) and ( 2) to maintain errors from the summed runoffs over the total period, |ΣQ sim /ΣQ obs | at ≤ 3% for each river.These two criteria should lead to both quantitative and qualitative reproducibility of the temporal variations of runoff and of the budgets of discharge volume.The initial conditions of the snow field and water level given to each tank model were derived by conducting a cyclic simulation for the year 2006.After the model calibration, to determine the parameters, more cyclic simulations were conducted to derive better initial conditions.This routine was repeated twenty times to satisfy the above two objective criteria.Each tank model adopted the last derived values for initial conditions and ran from the beginning of March, 2006 to the end of December, 2009 without the initializations after each flood event.

RESULTS AND DISCUSSIONS
Figure 3a shows the consistency between the observed and estimated runoffs within the four watersheds during the entire period and indicates correlations at significant levels.Both correlation values and errors at Yurappu (0.63 and 1%), Osaru (0.87 and 2%), Oshamanbe (0.71 and 3%), and Nukibetu (0.88 and 0.1%) lend support to our approach.The Nash-Sutcliffe efficiency coefficient (Nash and Sutcliffe, 1970) ranged between 0.6-0.7 except for Yurappu (0.4), which falls into the range (0.36-0.75) of "satisfactory" performance (Motovilov et al., 1999).The estimated and observed seasonal maxima around May indicate large runoffs arising from snowmelt.As a result of the model calibration, the water equivalent of snow on the mountain areas must be evaluated as three times that of precipitation because the local precipitation given by GPV-MSM would be an underestimate owing to its spatially coarse resolution for complicated mountainous topography.The observed low water runoff associated with the base-flow were reproduced well along with the local peaks generated by intense rainfall after the snow melting period, as shown between July 2006 to February 2007.The temporal tendency after each peak agrees well with the observation and is mainly formed by the interflow discharged from the 3rd outlet.These facts confirm that HaRUM satisfactorily reproduces the daily series of river runoffs required for more realistic ocean circulation models.
All of the coefficients were determined from Equation (3) using F(A(N)), the optimized exponential-type relationship between watershed area and F-value (see SM6); the time series of all runoffs were estimated and summed to derive the total discharge from all watersheds, Q all rivers = Σ Q total (N).The robust time series of the total discharge (Figure 3b) shows a seasonal variation that is similar to that for each runoff individually (Figure 3a).The maximum freshwater loading to Funka Bay occurs with a discharge of approximately 300 m 3 s −1 in late April, which is comparable to the high-water runoffs of large rivers (e.g., Yodo river, 278 m 3 s −1 ).The annual mean of the total discharge was 106 m 3 s −1 .Horizontal heterogeneity in each runoff was evident (Figure 3c) for 26 main rivers (Figure S4).The combined annual mean runoff of the Yurappu (14% of the total input of all rivers), Oshamanbe (7%), Nukibetsu (8%), and Osaru (20%) rivers amounts to approximately half of the total discharge.
The impacts of using the HaRUM on the salinity and flow field in the bay were evaluated using a coupled model consisting of both HaRUM and a three-dimensional OGCM developed at Kyoto University (e.g., Ishikawa et al., 2009; see SM7 for detailed settings).The analysis period of 2008-2009 was selected since a significant contrast in summertime surface salinity was found in this period (Figure S8) which should be influenced by the runoff into the bay; most noticeably the summer-time salinity in 2009 was much less than in 2008.Two numerical experiments were conducted: one case used the data sets of the monthly mean climatological runoffs (MMCR) provided by Dai and Trenberth (2002), and the other used runoffs estimated by the HaRUM (Figure S9) over the analysis period.
The results given in Figure 4a show that the summertime clockwise circulation was simulated by the coupled model.The difference of the surface salinity maps between both cases of HaRUM and MMCR (Figure 4c) shows that the salinity in the case of HaRUM was less than in the case of MMCR and was thus closer to the observation (Figure 4b).The correlation coefficient between simulated and observed salinities in the case of HaRUM (0.72) was higher than that of MMCR (0.63, Figure S10).Thus, the surface salinity simulation was improved by interfacing the HaRUM with the OGCM. Figure 4c shows that the surface current velocity of the clockwise circulation in the case of HaRUM was approximately twice as large as that in the case of MMCR.This difference suggests that the daily runoff calculated by the HaRUM can increase the modeled velocity of the summertime clockwise circulation in the bay.A more realistic simulation of the clockwise circulation can be expected by employing the HaRUM since the circulation should be a geostrophic current driven by the buoyancy of the runoff (Satoh et al., 2003;Takahashi et al., 2010).
Furthermore, the runoff reproduced by the HaRUM induced the interannual variation of the clockwise circulation in the bay.The difference of the summertime surface velocity field between 2008 and 2009 (Figure 4d) suggests that the circulation in 2009 intensified along the northern coast.The difference of the surface salinity field indicates that the low salinity water along the coast in 2009 was less saline than in 2008.This can be validated through a comparison of chlorophyll a distributions estimated from MODIS/Aqua images (Figure S3) gathered as of June 23, 2008 andJune 25, 2009 because chlorophyll a concentrations can be a reasonable index of nutrient rich, low salinity water along the coast.The averaged summertime runoff (Figure S9) from April to August in 2009 (236 m 3 s −1 ) was much more than that in 2008 (172 m 3 s −1 ), which induced stronger clockwise circulation in 2009 and low salinity water along the coastal region.The discharge generated by snowmelt in the mountain areas made up the majority (74-77%) of the total runoff (Figure S9), while the discharge in 2009 (177 m 3 s −1 ) was much more than in 2008 (133 m 3 s −1 ).The interannual variation in summertime runoff was attributed to the difference in the water equivalent of snow in the mountain areas.The water equivalent of snow in 2009 was much more than in 2008 (Figure S11), as confirmed by snow accumulation observed by the Automated Meteorological Data Acquisition System (AMeDAS).The AMeDAS station,

SUMMARY
Hydrometeorological procedures (HaRUM) for evaluating daily runoff flowing into Funka Bay were proposed in a synergistic approach aimed toward improved real-time ocean simulations.The HaRUM employed a distributed tank model on the basis of the water mass and heat budgets derived by utilizing the GPV-MSM datasets.The simulated daily runoff reproduced the four main observed runoff well, thereby enabling us to establish an exponential-type relationship between the runoff rates and watershed areas.This estimated relationship was used in the derivation of all parameters for the distributed tank models employed subsequently to estimate the runoff in many other small rivers.The coupled model of the OGCM and the HaRUM successfully reproduced the surface salinity and flow fields in the bay associated with the high net total snowmelt runoff.We pointed out that the advantage of incorporating the HaRUM into an ocean simulation is that it reproduces more realistic salinity fields over the interannual timescale.However, the reproducibility of water circulation affected by runoff should be validated by direct observation, both from a fisheries perspective and from the viewpoint of oceanography.The riverine fluxes from small watersheds have been neglected so far in many ocean models because of their coarse spatiotemporal resolution.As a result, the spatial effects of these small rivers have not so far been discussed, even though they could be important for improving local predictions associated with fishery harvesting.However, recent research in coastal ocean modeling reaches very high resolutions of up to ~1 km and permits eddy-resolving simulation.This paper underlines the fact that implementation of hydrometeorological models into high-resolution ocean prediction systems is essential, not only to simulate more realistic states in semi-enclosed bays, but also to provide the key link between hydrological and oceanographic processes.

Figure 1 .
Figure 1.Map of the study area showing watersheds around Funka Bay.The two upper-left panels show the geographical location of the bay in Japan.The white lines and colored areas indicate river paths and watersheds, respectively.Open circles represent river observation stations with emphasis on four major watersheds (solid line).Gray lines indicate isobaths in meters.The dashed line defines the mouth of the bay.

Figure 2 .
Figure 2. (a) Flowchart for the estimation of river runoff and (b) structure of the tank model with three serial tanks.The meanings of the quantities depicted in the diagram are given in the text.

Figure 3 .
Figure 3. (a) Temporal variation of the observed (circles with solid lines) and simulated (bars) daily runoff for each river from 2 March 2006 to 31 December 2007, (b) estimated volume transport of total discharge of all rivers, (c) annualmean river runoff of 26 major rivers (Figure S4).

Figure 4 .
Figure 4. (a) Snapshots of the simulated salinity and flow fields at the surface (−4m) in the case of HaRUM in 23 June, 2008.(b) Horizontal distribution of the observed salinity field in 23 June, 2008.Black dots indicate observation stations conducted by R/V Ushio-Maru (Hokkaido University).(c) Differences of simulated salinity and flow fields in the case of HaRUM minus MMCR in 23 June, 2008.(d) Differences of summer-time salinity and flow fields in 2009 minus 2008 averaged over three months (15 May-15 August).Colors indicate salinity fields (a and b) and differences of salinity (c and d).