Estimation of long-term external nutrient loading from watersheds to Lake Biwa by a combined rainfall-runoff model and loading-discharge curve approach

External nutrient loadings to Lake Biwa were estimated using a combined tank model and loading-discharge curve approach. The model was applied to collective drainage basins of the lake’s Imazu (northwest), Hikone (northeast), and Otsu (south) areas. The hourly model was conducted using particular discharges from Kita (Ado) river, Takatoki (Ane) river, and Yasu River to obtain loading curves for phosphate (PO4) and silica (SiO2) by assimilating measured concentrations (2002–2003). The tank model was updated by adding an evapotranspiration routine and direct paths of groundwater discharges to the lake floor. The daily model was calibrated through analysis of water budget among the basin, inflow, lake and outflow, and then validated. The model was established and combined into a loadingdischarge curve to determine the long-term external nutri‐ ent loadings entering the lake (1980–2017). Seasonal varia‐ tion in nutrient loadings increased during spring and summer and decreased during winter. Annual phosphatephosphorus (PO4-P) loading ranged from 217 to 296 tons yr–1 in the North Basin and 45 to 76 tons yr–1 in the South Basin, while SiO2 loading fluctuated from 16,027 to 32,655 tons yr–1 and 2,518 to 5,490 tons yr–1 in the North and South Basins, respectively.


INTRODUCTION
Lake Biwa is the largest freshwater lake in Japan, providing water resources to more than 14 million people in Western Japan (Kumagai, 2008). It is a body of water known worldwide for its massive sink of SiO 2 (Goto et al., 2007) and has been classified a phosphorus (P)-limited lake (Tezuka, 1986). The magnitude of the SiO 2 sink is related to the loading of P from the watershed of the lake (Goto et al., 2007). The external nutrient loading rate is one of the biogeochemical processes that affects nutrient-limited primary production (Howarth, 1988), and this may lead to an Correspondence to: Huu Le Tien, Graduate School of Environmental and Life Science, Okayama University, 3-1-1 Tsushima-Naka, Kita-ku, Okayama 700-8530, Japan. E-mail: letienhuu@hueuni.edu.vn expansion in SiO 2 and P retention causing subsequent alteration of the aquatic ecosystem. According to a previous study, the ratio of groundwater leakage rate from the groundwater system into Lake Biwa to the loading from rivers is estimated to be more than 0.1 (Taniguchi and Fukuo, 1996). Therefore, to understand the nutrient cycle in the lake and the impact on algal biomass and phytoplankton community composition, evaluation of the long-term variation in external nutrient loading to the lake needs to consider loading not only from rivers but also from groundwater (Ding et al., 2019;Longley et al., 2019). Previous work has reviewed the role of external loading to water quality or eutrophication status in Lake Biwa (Hsieh et al., 2011), analyzed the water quality of several rivers that flow into the lake (Kunimatsu and Kitamura, 1981;Taniguchi and Tase, 1999), and discussed the relationship between inflowing river water quality and the geological environment of upstream areas (Morii et al., 1993). However, data on nutrient levels in rivers and groundwater flowing into the lake are sparse and assessments still lack a quantitative assessment of the long-term variation in external loadings. Furthermore, estimating nutrient loading from all rivers surrounding Lake Biwa requires a huge amount of spatial information and understanding of complicated processes. A loading-discharge curve simply provides the correlation between nutrient load (L) and water discharge (Q), which has been widely used for the estimation of nutrient transport (Galat, 1990;Alexander et al., 2002;Salvia-Castellví et al., 2005;Amano and Kazama, 2012). A rainfall-runoff model, that is a hydrological tank model, is commonly used for long-term runoff analysis and can represent flow distributions over a given time for each layer of the watershed area (Arifjaya et al., 2011). However, simultaneous application of combined rainfall data and loading curves for estimating long-term external nutrient loading including groundwater discharge has rarely been mentioned in previous hydrological studies (Nakayama et al., 2011;Iwata et al., 2013). Therefore, this paper focuses on and discusses the applicability of a combination method, which includes the rainfall-runoff model and the L-Q curve approach to evaluate the long-term external nutrient loadings via surface and groundwater discharge into the Lake Biwa.

METHODS
Lake Biwa consists of a monomictic North Basin (area: 616 km 2 , average depth: 45.5 m), and a warm polymictic South Basin (area: 58 km 2 , average depth: 3.5 m) (Tezuka, 1992) that are surrounded by about 115 Class-A rivers of various sizes that flow into the lake (see Figure 1). Only the Seta River and the Kyoto Canal provide outflows from the lake. The collective basins were divided into three parts in this study, namely Otsu, Hikone and Imazu sub-basin. Discharges increase during snowmelt (March-April), the Eastern Asian summer monsoon (June-July), and the typhoon season (September-October).

Tank model
Hydrological tank models (Sugawara, 1979) were used to trace the local water cycle in and around Lake Biwa. The model includes four serial storage tanks, which is a simple and efficient method to convert rain into discharge and flux paths for nutrients. Each basin has been given different rainfall and evapotranspiration rates using meteorological data on precipitation and air temperature provided by the Ministry of Land, Infrastructure, Transport and Tourism (MLIT, 2020) at Otsu, Imazu, and Hikone, respectively. Because rainfall was considered uniform in each sub-basin, Figure 1. The river systems of three catchments in Lake Biwa. Sampling rivers were divided into three catchments, Imazu (I), Hikone (H) in North Basin, Otsu (O) in South Basin. The filled circles indicate the class-A sampling rivers. the local rain was given according to the areal fraction to the partial sum of sampling rivers in each sub-basin. The discharge from each tank is described as in Figure 2. Parameters of the models are presented in Table SI, SII. The fourth tank was diverted to connect to a certain depth of the lake floor. Evapotranspiration is incorporated via subtraction from the top tank. Runoff from the side outlets of a storage tank is proportional to water head over the outlet, and infiltration is proportional to the water depth. Water balance and its component can be described as: q 12 = α 12 ℎ 1 − z 12 ; p 1 = β 1 ℎ 1 (4) q 2 = α 2 ℎ 2 − z 2 ; p 2 = β 2 ℎ 2 (6) where Q: discharge; Q r : river discharge; Q g : groundwater discharge (m 3 s -1 ). V: lake volume (m 3 ); A, S: drainage and lake surface areas (km 2 ); H: water level of the lake (m); Δt: Structure of the tank model. The surface discharge (q 11 + q 12 ), intermediate discharge (q 2 ), sub-base discharge (q 3 ), and base discharge (q 4 ), p is infiltration to the lower reservoir, α and β are the coefficients, h is the water level, z is the outlet height time interval (hr and d); p i , q i : infiltration and discharge (mm Δt -1 ); P, E: precipitation and evapotranspiration (mm Δt -1 ); h i : water depth of the tank; z i : the height of outlet above the base of a tank; α i , β i : the runoff and infiltration coefficient.
Following Dalton's model (Korzukhin et al., 2011), evapotranspiration rate has been established as: where E: evapotranspiration rate (mm s -1 ); ρ: air density (1.2 kg m -3 ); C E : bulk coefficient; U: average wind speed (m s -1 ); e s , e a : saturated and actual water vapor pressure of the air (hPa); p: air pressure (hPa), respectively. Because of difficulty in including wind speed, the factor C E U/p was considered as a constant, which was adjusted to obtain the proper ratio of total evapotranspiration to total precipitation (r E ) in equation (14). Basic parameters are defined as: to check budgets in the fraction of river water, groundwater (r E , r G ) and residence time of water (R L ) in a year.

L-Q curve
In this study, we built the L-Q curve based on consistent hourly nutrient and discharge data. The relationship among instantaneous nutrient load L (g s -1 ), weight concentration C (g m -3 ) and instantaneous discharge Q (m 3 s -1 ) was explored following the suggestion of Amano and Kazama (2012) as: Multiplying discharge with the above, the nutrient load is where k, m, and p are the parameters for normalized discharge.
Equation (17) represents the dimensional intercept of the regime curve at Q = Q 1 , the unit discharge where Q 1 is taken to be 1 m 3 s -1 so dimensionless discharge retains dimensional value of one. In equation (16), the exponent p is close to unity and it is possible to separate the loading by groundwater discharge as: The gamma (γ) is used to express extra loads due to groundwater discharge at the floor of the lake. There are three particular cases: γ = 0: No extra loading other than rivers γ = 1: No difference in concentrations with river water γ > 1: Higher concentrations in groundwater Tsurumaki and Kobayashi (1989) reported that the PO 4 -P concentration in groundwater at the three sites in Lake Biwa was more than 0.1 mg L -1 . In this study, we assumed that γ = 3 for PO 4 -P estimation and γ = 1 for SiO 2 estimation.
Only 54 among the 115 rivers were selected as sampling sites during the 2002-2003 field surveys (Table SIII). A partial drainage area of sampled rivers in sub-basins was summed up to 2,002 km 2 , which represents around 63% of the total drainage area of 3,601 km 2 . PO 4 and SiO 2 were determined using a colorimeter for the three sub-basins bimonthly: July (preliminary), September (50 points), and November (53 points) of 2002; January (50 points), March (54 points), May (53 points), and July (55 points) of 2003 during a project run by the Ministry of Environment (Harashima et al., 2006). Discharge corresponding with sampling time in the rivers was calculated as Q i = (A i /A N )/ Q N using an hourly scale tank model where Q i , A i , A N and Q N are calculated discharge, sampling area, sub-basin area and initial discharge, respectively. Instead of 155 basinwise models given individual river discharges, we considered a virtual river of the sub-basin that was calculated for the partial discharge of the individual sampling river in which nutrient concentrations were measured.
The tank model (on a daily-basis) was set up to evaluate long-term inflow discharge  and temporal changes in water level. The difference between outflow or runoff in Equation (1) from the North Basin to the South Basin and runoff of Otsu sub-basin can be compared to the measured water level of the lake to determine water balance and the reliability of the model. Measured water level data were collected from 1991 to 2015 at Mihogasaki station (South part).

Evaluation criteria
Statistical parameters were used to quantify model performance: the Nash-Sutcliffe efficiency coefficient (NSE), the percentage of bias (PBIAS), and root mean squared error-observations standard deviation ratio (RSR) (Gupta et al., 2009;Moriasi et al., 2007). The respective equations were set as follows: where Q obs is the observed discharge, Q cal is the calculated discharge, Q obs is the mean observed discharge, and n is the total number of observations.

RESULTS AND DISCUSSION
The hourly model calculation was established using hourly streamflow data across 2 years, 2002 and 2003 ( Figure 3). The calibration period was from July 11 of 2002 until December 31 of 2002 and the validation period from January 1 until July 31 of 2003. Due to the limited availability of discharge data in the catchment surrounding Lake Biwa, observation data obtained from MLIT focused on the Yasu river for the Otsu sub-basin, Kita river of the Ado river for the Imazu sub-basin, and Takatoki river of the Ane river for the Hikone sub-basin. Precipitation at Otsu, Imazu, and Hikone were used in the model based on the standard precipitation R and is expected to represent the
The second approach focused on estimating the discharge from 2002 to 2003 using an hourly model to investigate the L-Q curve. The concentration of PO 4 and SiO 2 were plotted against instantaneous discharge. The external loading to Lake Biwa was assessed using the L-Q curve equation in Figure 4, which exemplifies a basic component of PO 4 and SiO 2 that was identified in the three rivers' systems around Lake Biwa. The relationship between the hourly river discharge and nutrient loading is shown in Figure 4. From this regression, the mean L was calculated from the mean Q. Furthermore, the result showed that the incline of the SiO 2 line was greater than that of PO 4 . These results reveal the equation for calculating loading based on discharge as: Otsu sub-basin: Total three sub-basins: where L 1 is PO 4 load, L 2 is SiO 2 load, and Q, Q 1 , Q 2 , Q 3 are the total discharge, partial discharges at Otsu, Imazu, and Hikone, respectively. After establishing the loading L-Q curve as equations (28) and (29) for a total of 3 sub-basins, the riverine nutrient loadings were estimated. The nutrient discharge from Imazu and Hikone sub-basins were used to integrate the external nutrient inflows into the North basin of Lake Biwa. Meanwhile, the nutrient loadings into the South basin were based on the calculation from the Otsu station. Ferguson (1986) mentioned that a loading method that uses a log-log Figure 4. The loading curve for PO 4 and SiO 2 which were regarded as unified for these nutrients. The correlation coeffecients between L and Q were R 2 = 0.99 (SiO 2 ) and R 2 = 0.97 (PO 4 ) loading curve can underestimate a river's loads from a statistical point of view. However, the log-log equation was found to minimize the deviation from the logarithmic data set, although a negligible scattering exists about the regression equation in the low region for PO 4 and SiO 2 as illustrated in Figure 4. Each straight line, which indicates a regression equation, fits the distribution well. The loadflow relationship for SiO 2 was higher than PO 4 indicating the sensitivity of PO 4 to water discharge.
The runoff rates (mm d -1 ) into the North and South basins were estimated and compared with the water level in Lake Biwa. The performance of the daily scale model with respect to the simulated water level was further examined using statistical criteria applied to the calibration and validation periods. Model calibration and validation statistics, which compared observed and simulated water levels for the daily time intervals, are displayed in Table SIV. Daily water level data across 13 years (1991 to 2003) were used to calibrate the remaining data from 2004 to 2015 and to validate the model performance ( Figure S1). The calibration and validation periods indicated that the simulated water level matched well with observations. A statistical comparison revealed NSE, RSR and PBIAS indices for the calibration (0.98, 0.15, and -0.01%) and validation (0.94, 0.25, and 0.25%) periods. The performance is in a 'good' category for both cases according to Moriasi et al. (2007).
The L-Q curve and the daily tank model were used to estimate external nutrients into Lake Biwa and establish the water budget. The annual river budget was 3.1 km 3 yr -1 , while the annual groundwater budget was 0.8 km 3 yr -1 (Table I). Fujino (1980) estimated the river and groundwater inflows to be 3.02 km 3 yr -1 and between 1 and 1.5 km 3 yr -1 , respectively. Figure 5a shows that the estimated annual PO 4 -P (converted from PO 4 ) loading increased from 1980, dropping in 1994 due to the drought in the region, and then fluctuated after 2000. Meanwhile, SiO 2 loading fluctuated from 1980, decreased in 1994 and 1998, and then increased afterwards. PO 4 -P loading into the North basin of Lake Biwa ranged from 217 to 296 tons yr -1 while loading in the south part averaged around 63 tons yr -1 , ranging from 45 to 76 tons yr -1 . PO 4 -P loading in the North basin is also considered as a source of nutrients for the South basin. According to Nagare et al. (2001), P mass measured from 1995 to 2000 was around 190 to 340 tons, and the output of P occupied 15% of the total budget (North and South Basin). In this study, the total PO 4 -P loading from the river in the North basin accounted for 81% of the total PO 4 -P loading. In addition, Goto et al. (2007) reported that the inflow discharge of silica was calculated to be 25,000 tons yr -1 , while the output was approximately  5,000 tons yr -1 . These findings are lower than our estimation of 30,244 tons yr -1 of average SiO 2 inflow, with 82% of the annual SiO 2 input to the watershed from loading of the North part where it ranged from 16,000 to 33,000 tons yr -1 . SiO 2 loading from the South ranged from 2,500 to 8,000 tons yr -1 . The estimation for the SiO 2 sink represented by Goto et al. (2013) showed that 20,000 ton SiO 2 yr -1 was retained in Lake Biwa, which accounts for 80% of the annual inflow. Figure 5b shows the estimated seasonal change in nutrient loading from rivers and groundwater. Loadings increased from March to June and reached a maximum in July, then decreased in winter. Agricultural practices such as the application of fertilizers and tillage might be one of the most important reasons for increasing nutrient loadings in rivers and groundwater during the spring and summer seasons. The rainy period from the beginning of June to late July, and also the typhoon periods in September and October contributed to the increased nutrient loading into Lake Biwa. We hypothesized that the PO 4 -P and SiO 2 discharges into Lake Biwa by rivers and groundwater flowing into the western area of the lake were smaller than those in the eastern areas and southeast areas. The difference in geological environments, land use, population density, and agricultural activities were the key criteria for our classification into 3 separate sub-basins. Nutrient loading in Hikone (western side) was the highest, while PO 4 -P loading in Imazu (eastern side) was lowest and SiO 2 loading was lowest in Otsu (southern side). We suggest that the urbanized and farmland areas give rise to higher loads of nutrients than forested area. Detailed analysis of interrelationships between the long-term variations in external nutrient loadings and nutrient cycling in the lake are tasks for the future.

CONCLUSIONS
In the present study, the long-term external nutrient loading into Lake Biwa was described by the combination of a hydrological tank model and L-Q curve in three steps. The hourly tank model was calibrated on a discharge basis by replacing central rain data with local rain. The local discharges of individual rivers were used to establish loading curves which were then used for estimating the total discharge of each sub-basin. The exponents of loading curves were close to unity, which suggests the load is proportional to the discharge (step 1). The daily tank model was calibrated on a water-level basis and the stage-discharge rating curve for the outflow (1991)(1992)(1993)(1994)(1995) showed that it was highly controlled and less varied. The second model was validated using water level data  under natural outflow conditions. Evapotranspiration and groundwater discharge were also considered (step 2). The daily-basis tank model was extended  and established to determine the long-term external nutrient loadings into the lake and were also used to analyze water budgets among the catchments. The results revealed that nutrient loadings increased during spring and summer and decreased during winter. Annual PO 4 -P loading ranged from 217 to 296 tons yr -1 in the North Basin and 45 to 76 tons yr -1 in the South Basin, while SiO 2 loading fluctuated from 16,027 to 32,655 tons yr -1 and 2,518 to 5,490 tons yr -1 in the North and South Basins, respectively.