Sensitivity Analysis on Lake Biwa under the A1B SRES climate change scenario using Biwa-3D Integrated Assessment Model:

Sensitivity analysis on lake temperature under the A1B SRES climate change scenario has been evaluated by coupling MRI-GCM with 20 km outputs and the Biwa-3D integrated assessment model. The non-hydrostatic 3D hydrodynamic model featuring Very Large Eddy Simulation (VLES) has been verified and compared with field observations in 2002 (the present-day year, PDY). A significant temperature increase in surface water, recording more than 34 degrees in the East Coast of the North Basin of the Lake Biwa, has been projected for the simulation for the very hot year (VHY), determined from comparison with MRI-GCM output for the year 2099, which may induce a catastrophic impact on lake water quality during the period. Weak stratification is predicted to start from March in VHY, compared with observations of stratification that doesn’t commence until after April in PDY. The thickness of the epilimnion which is around 15–20 m in PDY increases to 25–30 m in August in VHY due to a much higher atmospheric temperature. According to the model outputs, the vertical mixing may not always decrease due to the accumulated heat in the hypolimnion which has a higher temperature than the temperature during the winter. Accumulated heat in the hypolimnion may induce catastrophic degradation in the lake ecosystem.


INTRODUCTION
Several researchers discuss the issue of the impact of global warming on lake ecosystems; Matzinger et al. (2007) indicate through the survey of Lake Ohrid that global warming may cause an increase of nutrient input into lakes leading to eutrophication. Stich and Brinker (2009) suggest that the phosphorus concentration in deep and stratified lakes may decrease -whereas the temperature increases due to global warming. Lake Biwa is located in the middle of Japan and is an ancient lake which has been utilized as one of the most popular lakes providing drinking water for a population of 14 million people. The lake is under severe pressure arising from degradation of its ecosystem induced by several negative factors -increase of nutrient input from tributaries, decrease in dissolved oxygen in the bottom of the north basin, introduction of exotic species, and introduction of micropollutants in the lake water through rice herbicides .
The physical characteristics of Lake Biwa have been well studied and modeled by several researchers. Yamashiki (1999Yamashiki ( , 2000 and Yamashiki et al. (2003aYamashiki et al. ( , 2003b introduced VLES in predicting gyres formation using differential heating concepts. Shimizu et al. (2007) reproduced gyres and characteristics of internal waves using only meteorological conditions, particularly focusing on wind through the use of ELCOM (Estuary and Lake Computational Model). Akitomo et al. (2009aAkitomo et al. ( , 2009b) performed a complete survey in predicting temperature and current distribution using RANS and discussed the potential causes of development of counterclockwise (during the summer season) and clockwise (during the winter season) gyre formation associated with the development of a thermocline.
However, no research has previously been carried out which has coupled GCM output with a model of the Lake Basin. It should be cited that the GCM outputs initially contain uncertainty induced by (a) employed SRES scenario, (b) spatial resolution of the model, and (c) modeling (physical and numerical) scheme. However at AR4 the outputs of GCM are of a suitable quality to use as a subject of scientific discussion even though they still contain uncertainty. The difficulty in applying GCM output into the Lake Basin might be considered as being due to the following reasons: (i) Spatial scale: the resolution of GCM output is not sufficiently fine compared with the target basin scale, i.e most of lake catchment areas are too small compared with GCM output, and (ii) Anthropogenic impact in catchment: Impacts of human activities and catchment conditions can not be predicted precisely.
At AR4 the spatial resolution of GCM output becomes fine enough to overcome issue (i) as some of the modeling output is available with a grid spacing of less than 50 km for the whole globe. Issue (ii) is still a major issue and should not be ignored; however it does not prevent us from making predictions.
In this study, the primary goal is the sensitivity analysis of global warming on a lake environment, particularly on lake water temperature, using Biwa-3D featuring MRI-GCM with 20km grid spacing (Kitoh et al., 2009). Using this MRI-GCM one can predict the impacts of global warming at a resolution that is fine enough to run the model.

Governing Equations
The grid-filtered, incompressible Navier-Stokes and mass transport equations in a rotating environment can be expressed as follows: (1) where x i is eastward (i = 1), northward (i = 2) and upward (i = 3) directions, f is the coriolis parameter, u i (i = 1, 2, 3) is water velocity for each direction, β is the scalar, μ is the dynamic viscosity, κ = k/ρ 0 where k is the molecular diffusivity of the water and ρ 0 is a reference density of the water, P is water pressure, ρ is water density, F is gravity force, and Ps is the source term of the scalar.
The unresolved subgrid-scale (SGS) dynamic stress λ ij and flux χ j are modeled using the Mixed Scaling Formulation Model given below (Yamashiki, 1999(Yamashiki, , 2000Yamashiki et al., 2003aYamashiki et al., , 2003b: (4) in which the model coefficient C µ can be obtained from the following equations: in which is the grid-filter operator and is the testscale filtering operator.
The description of all variables in equations (4) and (5) are shown in Document S1 in Supplement S2.
Governing equations were discretized using a Semi-Lagrangian scheme (SLS) (Kalnay, 2003). A staggered scheme is generally used for the most part of the program, except for variables related to Dynamic SGS Modeling that were determined in the middle of each cell. The Adams-Bashforth method was employed for time integration. All three dimensional variables are transferred to one dimensional list vectors in the coding to reduce memory consumption.

Effect of Aspect Ratio in VLES
Vertical convection which can be explicitly calculated according to the above mixing equation, may fail when the aspect ratio is large. In this case study, we employed an adjustment coefficient to circumvent insufficient vertical mixing during the mixing procedure as (6) in which A R : aspect ratio determined by horizontal/vertical grid spacing; n: mixing adjustment constant, determined through calibration with observation results. The introduction of the aspect ratio into the vertical mixing coefficient may facilitate a gradual change from a realistic vertical mixing coefficient when the aspect ratio is small to a situation where the subgrid-scale model over-estimates the vertical mixing coefficient when the aspect ratio is high.

Verification of the model for present-day year
In this survey, we verified the model performance through a comparison with field observation data obtained by Lake Biwa Environmental Research Institute using a Fine Scale Profiler (F-Probe). Ten observation points were stationed along a north-south cross-section of the lake, shown in Figure S1 in Supplement 1. The observed data are displayed using a three-dimensional color contour map created by three-dimensional interpolation based on the inverse-distance method which takes the stratification process into consideration in order to illustrate horizontal and vertical distribution of the lake temperature.
We determined the year 2002 as the present-day year (hereafter referred to as PDY) used for the model validation purpose. The average temperature and precipitation for Lake Biwa basin during the period of 1985-2004 is 15.4 ± 0.5 degrees Celsius and 3.9 ± 0.6 mm/day, respectively. The average temperature of the year 2002, employed in this study as PDY, is higher than average by only 0.2 degrees Celsius and thus can be considered as representative year of the current temperature trend. The precipitation of the year 2002 shows 4.5 mm/day, which is the same as the mean plus the standard deviation, can be considered as slightly high precipitation compared with the 1985-2004 average.
For the surface boundary condition generation, 6 stations of the Automated Meteorological Data Acquisition Systems (AMeDAS) are selected as forcing dataset for PDY in order to generate the necessary vorticity through wind and to cover the whole lake area using inner interpolation. We employed (a) a 500 m horizontal and 2 m vertical grid spacing corresponding to 83 × 128 × 51 cells in total and (b) a 250m horizontal and 2m vertical grid spacing corresponding to 166 × 256 × 51 cells in total. Figure 1 illustrates a comparison of simulated monthly lake temperature with monthly field observation data provided by Lake Biwa Environmental Research Institute calculated using a 500 m horizontal grid spacing. The same comparison using a 250 m horizontal grid spacing is also shown in Figure S3 in Supplement 1. The tendency of the temperature distribution obtained from the field observation can be illustrated using numerical simulation. At the same time, the sharp boundary of stratification during the summer time in the lake margin can not be reproduced well and weak stratification may be found in the simulation. However, this might be caused by the interpolation procedure using the observation data obtained at the center of the lake, not at the margin of the lake.

Determination of the forcing dataset-VHY
For understanding the environment of Lake Biwa under warm climate condition in this study, we attempted to obtain very hot year (hereafter referred to as VHY) dataset as projected global warming signals, determined from comparison with 20-km grid GCM outputs produced by MRI (hereafter referred to as MRI-GCM) for the year 2099 through the Kakushin program promoted by the Ministry of Education and Culture (Kitoh et al., 2009). The SRES A1B scenario has been employed in this survey as this scenario is being considered to be representing the current trend of emissions. Figure S4 (a) in Supplement 1 illustrates the comparison between spatial averaged temperature among the 6 AMeDAS stations surrounding Lake Biwa and among 6 corresponding MRI-GCM grid-points in the calculation for 2001-2002. The annual-mean temperature difference between those results was −1.79 degrees Celsius, showing that the MRI-GCM under-estimates temperature for 2001-2002. Figure  S4 (b) in Supplement 1 shows the comparison of spatial averaged temperature among the 6 corresponding MRI-GCM grid-points between 2002 and 2099. Although the peak temperature in August does not rise too much in 2099, we observe that the duration of high temperature in 2099 increases compared with that in 2002 for both the AMeDAS and MRI-GCM outputs. The calculated annual-mean temperature rise between the target years is 3.3 degrees Celsius.
The average temperature created by the forcing data using above described procedure is 18.9 degrees Celsius, 3.5 degrees higher than the average temperature of the current climate. We classified this year as the very hot year (VHY), as the deviation is seven times the standard deviation of the 1985-2004 average temperature.
We made this comparison for several years and considered this difference as the important global warming signal obtained from this GCM output. We therefore imposed those increases into AMeDAS station data for the estimation of the lake status under VHY and obtained forcing dataset for VHY. Similar analysis has been made for the precipitation. GCM output values are then converted to estimated future AMeDAS station values using spatial and temporal interpolation.

Sensitivity analysis for the very hot year
Biwa-3D was applied for understanding lake environmental changes under the VHY. Two-years spin-up calculation has been made before applying numerical simulation for the corresponding year. Inflowing river conditions for daily discharge from 42 streams have been provided by modeling output using the Storage Function method by Biwako Public Works Office of Ministry of Land, Infrastructure, Transport and Tourism (MLIT), government of Japan. The river discharges for those 42 streams are adopted in association with the future precipitation calculated by the MRI-GCM.
3D temperature distribution is calculated for PDY and VHY by using the forcing datasets and shown in Figure S7 in Supplement 1. High lake surface water temperature is simulated in the summer season for the north basin of Lake Biwa. In the PDY simulation, the average lake surface water temperature in the north basin is around 29 degrees Celsius and only the east coast exceeds 30 degree Celsius. In VHY the whole lake surface water temperature during July-September exceeds 30 degrees Celsius and, according to the temperature distribution calculated using MRI-GCM outputs the maximum temperature exceeds 34 degrees Celsius.
The Lake surface temperature drops below that of the hypolimnion in the January in VHY prediction due to the accumulated heat during the spin-up and one-year calculation, which may induce vertical mixing if the lake surface water becomes cooler than that of the lake body, whereas the lake surface water becomes slightly cooler than the hypolimnion in PDY. This may change due to the increase in air temperature during winter and accumulated heat during summer. We assumed that the vertical lake water mixing always decreases due to VHY condition: however it may vary depending upon the projected winter air temperatures. On the other hand, the higher temperature in the hypolimnion during winter may induce significant degradation in the ecosystem which may lead to nutrient release, decrease of oxygen concentration, and increase in harmful phytoplankton in the epilimnion. These factors should be addressed in Part II of this series. Figure 2 and 3 illustrate calculated monthly averaged lake surface water temperature in PDY and VHY in the overall lake and vertical section, respectively. By taking a look at the July and August vertical temperature distribution in VHY in Figure 3, the accumulated heat during the summer season raises the lake water temperature to almost 30 degrees Celsius for the overall epilimnion (0-25 m depth). Also it should be noted that the thickness of the epilimnion which is around 15-20 m in PDY increases to 25-30 m in August of VHY. Entire vertical mixing in the lake occurs during January in PDY, whereas in 2002 due to the accumulated heat in the hypolimnion it is simulated to occur after December when surface water in the epilimnion becomes cooler than that of the hypolimnion. In VHY we may observe weak thermocline formation in March which is only observed after April in PDY.
However, the vertical mixing status may vary according to the air temperature in the winter period depending upon the choice of the projected climate. It should also be noted that the model shows over-cooling in the beginning of winter according to the comparison with observation data, and thus the continuation of the winter-season's vertical mixing, which is essential for all aquatic life, is still uncertain.

CONCLUSION
The sensitivity analysis of A1B SRES climate change scenario on Lake Biwa has been performed using the Biwa-3D integrated assessment model. Initially, lake water temperature has been projected for VHY, showing a significant temperature increase during summer. In winter due to accumulation of heat the temperature of the hypolimnion does not decrease significantly. Vertical  mixing, which supplies cooler, oxygen-rich water into the hypolimnion, still occurs due to surface cooling -with warmer atmospheric temperature in the winter compared with PDY -as the hypolimnion temperature become higher than that of PDY.
Also, it should be noted that the stratification in VHY becomes stronger than that of PDY with the increase of the thickness of the epilimnion, which is around 15-20 m in PDY increased to 25-30 m in August in VHY. The more pronounced stratification may cause severe water quality degradation during the summer time especially for dissolved oxygen at the bottom.
It should also be noted that the modeling results may vary as both MRI-GCM output and Biwa-3D contain certain biases and uncertainty.
The sensitivity analysis of lake water quality focusing on the dissolved oxygen concentration is scheduled for the next paper.

ACKNOWLEDGEMENT
The authors wish to thank Dr. Toshiyuki Nakaegawa at the Meteorological Research Institute for his valued advice, and Mr. Kenichi Tatsumi at DPRI Kyoto University for extracting climate data from MRI-GCM outputs.

SUPPLEMENTS
Supplement 1 Section 1 FIGURES Figure S1 Observation points (white circles) and 3D projection of vertical profile measurements using Fine Scale Profiler (F-Probe) in Lake Biwa. Figure