Hydrological and Chemical Budgets of Okama Crater Lake in Active Zao Volcano, Japan

: The Okama Crater Lake is located in the highly active Zao Volcano on the boundary of Miyagi and Yamagata Prefectures, Japan. At present, the lake stays relatively calm with neither bubbling, steaming nor gas smell at a pH of 3.2–3.4, though the lake did change color with steaming from the water surface in 1939 as the result of one of Zao’s volcanic activities. In order to clarify the geothermal effect on Okama and the groundwater ﬂow system below or around Okama, ﬁeld observations were performed in 2019 and 2020. Groundwater inﬂow and outﬂow in Okama were separately evaluated by estimating the hydrological and chemical budgets of the lake, based on the hydrometeorology, water temperature and river inﬂow measured in the ﬁeld. The average groundwater inﬂow and outﬂow were estimated at 0.012 m 3 /s and 0.039 m 3 /s during the non-rainfall periods of 2020, respectively. A surplus of groundwater outﬂow makes the lake level consistently decrease during non-rainfall periods or the completely ice-covered season. In the completely ice-covered periods, the water temperature consistently increased at 0–15 m above the lake bottom, which is probably due to thermal leakage from a hydrothermal reservoir below the lake bottom. The heat ﬂuxes averaged over December 2019–April 2020 and December 2020–March 2021 were calculated at 2.5 and 2.9 W/m 2 , respectively. A coupling between the estimated groundwater inﬂow and the calculated geothermal heat ﬂux was used to evaluate the temperature of inﬂowing groundwater.


Introduction
When a volcano has high potential to produce a disaster, its volcanic activity is monitored at present by video cameras, seismometers, gravimeters, tiltmeters and a GNSS (Global Navigation Satellite System) survey (US Geological Survey; URL https://www.usgs.gov/ media/images/volcanic-monitoring-types-and-methods-employed-usgs-vhp (accessed on 4 December 2021)). However, the prediction of volcanic eruptions is sometimes very difficult and, consequently, has produced many serious disasters, even in the 21st century (e.g., the 2010 eruption of Merapi Volcano, Indonesia [1,2], and the 2014 eruption of Mt. Ontake, Japan [3]). In order to increase the prediction accuracy and improve understanding, thermal and chemical studies are effective additions to the above geophysical approaches, although continuous thermal and chemical monitoring is rarely used among volcanoes because the sensors suffer heavy damage due to the high temperatures and/or corrosive environments.
If a lake or pond exists in an active volcano, the volcano carries a high risk of producing phreatic/phreatomagmatic explosions and lahars at the time of eruptive activity (e.g., Taal Volcano, the Philippines, with Taal Lake [4], Kusatsu-Shirane Volcano, Japan, with Yugama headwater region of the Nigori River, where part of Okama's water is likely to leak [14]. However, the location and altitude of the springing points are not specified in Konno [14]. The dot-dashed line shows a water divide of Okama, and a dashed line represents the uppermost part of a water divide of the Nigori River. The (c) photo shows the location of site M (black circle) and site MD (red circle). Water sampling was performed at site MD (lake water), site N (Nigori River water) and the three inflowing streams ( Figure 2).
Okama has three inflowing streams: the Goshiki Stream, the Western Stream and the Southern Stream. Of all the streams, only the Goshiki Stream is perennial with streaming water at any time; the other streams are ephemeral with streaming water only during a rainy period or in the snowmelt season. Referring to Konno [14] and Yoshimura [15], the lake was iron-rich and acidic at pH of 2.6 and 37 m deep at maximum in 1936. In 2020, the lake was still iron-rich, and the pH and maximum depth were 3.1-3.3 and ca. 25 m, respectively. Thus, the lake type of Okama belongs to a siderotrophic lake within inorganic acidotrophic lakes [15]. On the topographic map with a 1/25,000 scale, in 2014, the highest lake level is at an altitude of ca. 1550 m asl, which is the same as the level in 1936 [14]. Hence, the decrease in maximum depth from 37 m to 25 m indicates that sedimentation is in progress, probably due to surrounding slope collapse or sediment runoff from the inflowing streams. The dot-dashed line shows a water divide of Okama, and a dashed line represents the uppermost part of a water divide of the Nigori River. The (c) photo shows the location of site M (black circle) and site MD (red circle). Water sampling was performed at site MD (lake water), site N (Nigori River water) and the three inflowing streams ( Figure 2).
Okama is located at an altitude of ca. 1550 m above sea level (asl) with a depth of ca. 25 m at maximum and water surface area of 8.208 × 10 4 m 2 ( Figure 2). The drainage area, including the lake area, is 6.312 × 10 5 m 2 ( Figure 1). The drainage basin is adjacent to the headwater region of the Nigori River, where part of Okama's water is likely to leak [14]. However, the location and altitude of the springing points are not specified in Konno [14].
Okama has three inflowing streams: the Goshiki Stream, the Western Stream and the Southern Stream. Of all the streams, only the Goshiki Stream is perennial with streaming water at any time; the other streams are ephemeral with streaming water only during a rainy period or in the snowmelt season. Referring to Konno [14] and Yoshimura [15], the lake was iron-rich and acidic at pH of 2.6 and 37 m deep at maximum in 1936. In 2020, the lake was still iron-rich, and the pH and maximum depth were 3.1-3.3 and ca. 25 m, respectively. Thus, the lake type of Okama belongs to a siderotrophic lake within inorganic acidotrophic lakes [15]. On the topographic map with a 1/25,000 scale, in 2014, the highest lake level is at an altitude of ca. 1550 m asl, which is the same as the level in 1936 [14]. Hence, the decrease in maximum depth from 37 m to 25 m indicates that sedimentation is in progress, probably due to surrounding slope collapse or sediment runoff from the inflowing streams.
Hydrology 2022, 9, x FOR PEER REVIEW Figure 2. Bathymetric map of Okama Crater Lake with inflowing streams. The thick solid grey lines and dotted lines in the map show isopleths of water depth at 5 m, 1 m and 0.2 m respectively (modified after Yamasaki and Yagi [16] as well as Tsuchiya and Hirano [17])

Estimates of Hydrological and Chemical Budgets
When a crater lake is located in an active volcano, it is important to know groundwater in or around the lake is connected to the underground hydrother system and magma. Thus, in order to separately evaluate the groundwater inf outflow in Okama, the hydrological and chemical budgets of the lake were estima hydrological budget equation for a closed lake with no outflowing river, such as is given as follow: ΔV/Δt = (P-E)A0 + Rin + Gin-Gout where V is the water volume (m 3 ), A0 is the water surface area (m 2 ), P is the prec (m/s) onto the lake surface, E is the evaporation (m/s) at the lake surface, Rin is inflow (m 3 /s), Gin and Gout are the groundwater inflow and outflow (m 3 /s), resp and t is the time. The left side of Equation (1) indicates the water volume ch budget period Δt. Evaporation E was calculated by using the following bulk method: where QE is the heat flux (W/m 2 ) by evaporation, λ is the latent heat (J/kg) for evap β is the ratio of water vapor density to dry air density (=0.622), aE is the dimension transfer coefficient for latent heat, uz is the wind speed (m/s) at z (m) above the face, p is the air pressure (Pa) at z, ez is the vapor pressure (Pa) at z, e0 is the s vapor pressure (Pa) at lake surface temperature Ts (K), and ρa and ρw are the air a densities (kg/m 3 ), respectively. Here, z= 4.0 m was adopted in Equation (2) as th  [16] as well as Tsuchiya and Hirano [17]).

Estimates of Hydrological and Chemical Budgets
When a crater lake is located in an active volcano, it is important to know how the groundwater in or around the lake is connected to the underground hydrothermal flow system and magma. Thus, in order to separately evaluate the groundwater inflow and outflow in Okama, the hydrological and chemical budgets of the lake were estimated. The hydrological budget equation for a closed lake with no outflowing river, such as Okama, is given as follow: where V is the water volume (m 3 ), A 0 is the water surface area (m 2 ), P is the precipitation (m/s) onto the lake surface, E is the evaporation (m/s) at the lake surface, R in is the river inflow (m 3 /s), G in and G out are the groundwater inflow and outflow (m 3 /s), respectively, and t is the time. The left side of Equation (1) indicates the water volume change per budget period ∆t. Evaporation E was calculated by using the following bulk transfer method: where Q E is the heat flux (W/m 2 ) by evaporation, λ is the latent heat (J/kg) for evaporation, β is the ratio of water vapor density to dry air density (=0.622), a E is the dimensionless bulk transfer coefficient for latent heat, u z is the wind speed (m/s) at z (m) above the lake surface, p is the air pressure (Pa) at z, e z is the vapor pressure (Pa) at z, e 0 is the saturated vapor pressure (Pa) at lake surface temperature T s (K), and ρ a and ρ w are the air and water densities (kg/m 3 ), respectively. Here, z= 4.0 m was adopted in Equation (2) as the height of the wind speed sensor as well as of the air temperature and relative humidity logger above the lake water surface. The dimensionless bulk transfer coefficient, a E , for latent heat Hydrology 2022, 9, 28 5 of 18 was given at a constant of 0.0013 for z = 4.0 m [18]. Under the condition of non-rainfall, Equation (1) is as follows: Here, the G value, or the net groundwater inflow, is provided via calculation of the water volume change and surface water area change from the measurement of lake level (Figure 2), the measurement of river inflow and the calculation of evaporation E by applying the hydrometeorological data at site M to Equations (2) and (3).
The chemical budget equation for the closed lake is as follows: where C L is the mean ionic concentration (g/L) of the lake; C Rin , C F and C Gin are the ionic concentrations (g/L) of inflowing rivers, precipitation and inflowing groundwater, respectively; and S is the net depositional flux (kg/s) related to a chemical reaction of the ion. The mean ionic concentration, C L , of the lake is given as that of groundwater outflow, since the depth at which the lake water leaks out as groundwater is unknown. The magnitude of S is inferred by considering the lake water chemistry based on the ionic analysis. If a non-rainfall period is adopted as the budget period, the second term on the right side of Equation (5) is zero. Then, the simultaneous equations from Equations (4) and (5) lead to: where B =∆(C L V)/∆t − C Rin R in + S. Here, the B values are given by a coupling with the analyzed chemistry of the stream and lake waters. Then, an ion that is not influenced by the chemical reaction should be selected as the ion specified in the ionic concentration in Equation (5), since quantifying the flux S is sometimes difficult. However, the flux S, if any, could be estimated by setting sediment traps on the lake bottom and then by identifying the mineralogy or chemical compounds of the trapped deposit with the X-ray diffraction (XRD) method.

Field Observations
In order to evaluate the temporal change in the water volume of Okama on the left side of Equation (1), lake level was measured at 1 h intervals by fixing two pressure loggers (type HOBO U20, Onset Computer Inc., Bourne, MA, USA) at site L (for water pressure) and site M (for air pressure) ( Figure 1). The difference between the water pressure and the air pressure provided the water depth at site L as the lake level. The pressure logger at site L was fixed throughout the year.
In order to calculate the first term on the right side of Equation (1) and evaporation E in Equation (3), a weather station (URL https://www.onsetcomp.com/products/kits/ u30-nrc-sys-c/ (accessed on 17 November 2021)) was set at site M to record air pressure, air temperature, relative humidity, wind velocity and rainfall at 1 h intervals. The weather station was removed from site M before snowfall began in October, because the snow accumulation becomes more than 1 m deep in winter. Meanwhile, in order to measure the hourly air temperature, relative humidity and air pressure throughout the year, two loggers were fixed at the Daikokuten observation hut (altitude, 1446 m asl), belonging to the Research Center for Prediction of Earthquakes and Volcanic Eruptions, Graduate School of Science, Tohoku University, which is 1.55 km southeast of site MD. Therefore, the daily mean air temperature and air pressure at Okama after the removal of the weather station at site M were acquired using the high correlations (R 2 = 0.991, p < 0.01) for the daily mean air temperature and (R 2 = 0.982, p < 0.01) for the daily mean air pressure between site MD and the observation hut. The annual mean air temperature at site MD was also evaluated using annual data from the observation hut.
On the lake shore near site M, the discharge of the three inflowing streams ( Figure 2) was measured on non-rainfall days to evaluate the river inflow R in in Equation (1). The stream water was also manually sampled to analyze the inorganic chemistry, related to C Rin in Equation (5). The EC25 (electric conductivity at 25 • C) and pH of the stream water were then measured by a thermocouple thermometer (Type TX1001, YOKOGAWA, Co., Ltd., Tokyo, Japan) and a portable pH and EC meter (Type D-74, HORIBA, Co., Ltd., Kyoto, Japan). The Nigori River water was sampled at site N (Figure 1), its pH and EC25 were similarly measured, and the chemistry was analyzed. The river water at site N is unlikely to contain water leaked from Okama because its altitude (1557 m asl) is higher than that of Okama's water surface (1550 m asl).
In order to explore the effect of the geothermal heat flux on the thermal condition of Okama, the water temperature was measured every hour at nine depths of the deepest point (site MD) using a mooring system of temperature loggers (type TidBit v2 with accuracy, ±0.2 • C, Onset Computer Inc., Bourne, MA, USA; URL https://www.onsetcomp.com/ products/data-loggers/utbi-001/ (accessed on 17 November 2021)) ( Figures 1 and 3). In order to ascertain the temporal change in lake water chemistry related to the left side of Equation (5), an electric conductivity (EC) logger for fresh water (type HOBO U24-001 with accuracies of ±0.5 mS/m and ±0.1 • C for electric conductivity and temperature, respectively, Onset Computer Inc., Bourne, MA, USA) was added in June 2020 to a position of 2 m above the lake bottom. The mooring system was set in the lake throughout the year. were then measured by a thermocouple thermometer (Type TX1001, YOKOGAWA, Co. Ltd., Tokyo, Japan) and a portable pH and EC meter (Type D-74, HORIBA, Co. Ltd., Kyoto, Japan). The Nigori River water was sampled at site N ( Figure 1), its pH and EC25 were similarly measured, and the chemistry was analyzed. The river water at site N is unlikely to contain water leaked from Okama because its altitude (1557 m asl) is higher than that of Okama's water surface (1550 m asl). In order to explore the effect of the geothermal heat flux on the thermal condition of Okama, the water temperature was measured every hour at nine depths of the deepest point (site MD) using a mooring system of temperature loggers (type TidBit v2 with accuracy, ±0.2 °C, Onset Computer, Inc., USA; URL https://www.onsetcomp.com/products/data-loggers/utbi-001/ (accessed on 17 November 2021)) (Figures 1 and 3). In order to ascertain the temporal change in lake water chemistry related to the left side of Equation (5), an electric conductivity (EC) logger for fresh water (type HOBO U24-001 with accuracies of ±0.5 mS/m and ±0.1 °C for electric conductivity and temperature, respectively, Onset Computer Inc.) was added in June 2020 to a position of 2 m above the lake bottom. The mooring system was set in the lake throughout the year.
Meanwhile, by lowering a profiler on boat (type ASTD102 with accuracies of ±0.1 mS/m and ±0.01 °C for EC and water temperature, respectively, JFE-Advantech, Inc., Japan; URL https://www.jfe-advantech.co.jp/products/ocean-rinko.html (accessed on 16 November 2021)), the water temperature and EC25 were vertically measured at 0.1 m depth intervals in June-October of 2019 and 2020. The EC and water temperature were recorded only during the downward movement of the profiler. Thus, the profiler's record is not influenced by the disturbance of suspended matters from its bottom touch. The water temperature and EC25 from the loggers were calibrated by those from the profiler with relatively high accuracies. Lake water and bottom sediment at site MD were sampled on boat by a Van Dorn sampler (URL https://www.kc-denmark.dk/products/water-sampler/van-dorn-watersampler.aspx (accessed on 6 September 2021)) and an Ekman-Burge grab sampler [19], respectively. The water temperature, EC25 and pH of the lake water were then measured by a thermocouple thermometer (Type TX1001, YOKOGAWA, Co. Ltd., Tokyo, Japan), a compact EC meter (Type LAQUAtwin EC-11, HORIBA, Co. Ltd., Kyoto, Japan) and a compact pH meter (Type LAQUAtwin pH-11B, HORIBA, Co. Ltd., Kyoto, Japan; URL https://www.horiba.com/en_en/water-quality/pocket-meters/ (accessed on 18 November 2021)). Meanwhile, by lowering a profiler on boat (type ASTD102 with accuracies of ±0.1 mS/m and ±0.01 • C for EC and water temperature, respectively, JFE-Advantech, Inc., Japan; URL https://www.jfe-advantech.co.jp/products/ocean-rinko.html (accessed on 16 November 2021)), the water temperature and EC25 were vertically measured at 0.1 m depth intervals in June-October of 2019 and 2020. The EC and water temperature were recorded only during the downward movement of the profiler. Thus, the profiler's record is not influenced by the disturbance of suspended matters from its bottom touch. The water temperature and EC25 from the loggers were calibrated by those from the profiler with relatively high accuracies.

Laboratory Experiments
In order to obtain the concentration C Gin of inflowing groundwater in Equation (5), pore water in the bottom sediment of site MD was sampled by an extractor kit with 60 mL syringe and a porous cup (URL https://www.soilmoisture.com/1900K2/ (accessed on 24 November 2021)), and its EC25 and pH were measured by a compact EC meter (Type LAQUAtwin EC-11, HORIBA, Co., Ltd., Kyoto, Japan) and a compact pH meter (Type LAQUAtwin pH-11B, HORIBA, Co., Ltd., Kyoto, Japan), respectively.
Concentrations of major ions (K + , Na + , Mg 2+ , Ca 2+ , Cl − , SO 4 2− ) for the sampled river water, lake water and pore water were measured by ion chromatography (type Dionex ICS-1600 (cation) and ICS-2100 (anion), Thermo Fisher Scientific Inc., Tokyo, Japan; URL https://www.thermofisher.com/ (accessed on 24 November 2021)). Only the bicarbonate ion (HCO 3 − ) concentration was measured using sulfuric acid titration. However, the HCO 3 − concentration was almost zero for a few samples with pH of 2.6-4.3. Thus, for the samples at pH of 4.3 or less, the titration was omitted. In order to obtain the concentration CGin of inflowing groundwater in Equation (5), pore water in the bottom sediment of site MD was sampled by an extractor kit with 60 mL syringe and a porous cup (URL https://www.soilmoisture.com/1900K2/ (accessed on 24 November 2021)), and its EC25 and pH were measured by a compact EC meter (Type LAQUAtwin EC-11, HORIBA, Co. Ltd., Kyoto, Japan) and a compact pH meter (Type LAQUAtwin pH-11B, HORIBA, Co. Ltd., Kyoto, Japan), respectively.

Water Chemistry
Concentrations of major ions (K + , Na + , Mg 2+ , Ca 2+ , Cl -, SO4 2-) for the sampled river water, lake water and pore water were measured by ion chromatography (type Dionex ICS-1600 (cation) and ICS-2100 (anion), Thermo Fisher Scientific Inc., Tokyo, Japan; URL https://www.thermofisher.com/ (accessed on 24 November 2021)). Only the bicarbonate ion (HCO3 -) concentration was measured using sulfuric acid titration. However, the HCO3 -concentration was almost zero for a few samples with pH of 2.6-4.3. Thus, for the samples at pH of 4.3 or less, the titration was omitted. The Nigori River water (Figure 4e) was sampled in the catchment neighboring the Okama catchment ( Figure 1) and, thus, appears to be unaffected by the surface water and groundwater in the Okama catchment. However, its high acidity (pH = 3.2) suggests that the Nigori River catchment is also affected by sulfate substances as part of Zao Volcano.  The Nigori River water (Figure 4e) was sampled in the catchment neighboring the Okama catchment ( Figure 1) and, thus, appears to be unaffected by the surface water and Hydrology 2022, 9, 28 8 of 18 groundwater in the Okama catchment. However, its high acidity (pH = 3.2) suggests that the Nigori River catchment is also affected by sulfate substances as part of Zao Volcano. Figure 5 shows relations between the EC25 and (a) SO 4 2concentration, (b) Ca 2+ concentration, (c) Ca 2+ + SO 4 2− concentration or (d) the pH for sampled lake and stream waters. There are clear linear relationships with high correlations of R 2 = 0.746-0.914. Especially, the relation between the EC25 and SO 4 2concentration ( Figure 5a) is highly correlated at R 2 =0.914. These relations indicate that, if the EC25 is measured by the EC logger or the profiler, then the Ca 2+ and SO 4 2concentrations and pH are also numerically obtained with high confidence levels. Here, the pH appears to be controlled by the SO 4 2− concentration ( Figure 4).

Water Chemistry
Hydrology 2022, 9, x FOR PEER REVIEW 8 of 19 Figure 5 shows relations between the EC25 and (a) SO4 2-concentration, (b) Ca 2+ concentration, (c) Ca 2+ + SO4 2-concentration or (d) the pH for sampled lake and stream waters. There are clear linear relationships with high correlations of R 2 = 0.746-0.914. Especially, the relation between the EC25 and SO4 2-concentration ( Figure 5a) is highly correlated at R 2 =0.914. These relations indicate that, if the EC25 is measured by the EC logger or the profiler, then the Ca 2+ and SO4 2-concentrations and pH are also numerically obtained with high confidence levels. Here, the pH appears to be controlled by the SO4 2concentration ( Figure 4).  When the thermocline reached the bottom, the bottom temperature should have increased abruptly from 6.8 to 9.5 °C in 2019 or from 4.9 to 9.2 °C in 2020. This is because, during the cooling season, the thermocline descends to the bottom, while the upper layer keeps thermally uniform.  The EC25 values (and therefore the SO4 2-concentration) greatly increased downward at or below the thermocline. These great spatio-temporal increases in the EC25 appear to have been produced by the density underflow from the sediment runoffs of inflowing streams or the surface flow from rainfalls. In fact, such a sediment runoff was observed on the lake's shore after heavy rainfall (169 mm/day) on 28 July 2020 (Figure 8). The water  The EC25 values (and therefore the SO4 2-concentration) greatly increased downward at or below the thermocline. These great spatio-temporal increases in the EC25 appear to have been produced by the density underflow from the sediment runoffs of inflowing streams or the surface flow from rainfalls. In fact, such a sediment runoff was observed on the lake's shore after heavy rainfall (169 mm/day) on 28 July 2020 (Figure 8). The water When the thermocline reached the bottom, the bottom temperature should have increased abruptly from 6.8 to 9.5 • C in 2019 or from 4.9 to 9.2 • C in 2020. This is because, during the cooling season, the thermocline descends to the bottom, while the upper layer keeps thermally uniform.
The EC25 values (and therefore the SO 4 2− concentration) greatly increased downward at or below the thermocline. These great spatio-temporal increases in the EC25 appear to have been produced by the density underflow from the sediment runoffs of inflowing streams or the surface flow from rainfalls. In fact, such a sediment runoff was observed on the lake's shore after heavy rainfall (169 mm/day) on 28 July 2020 (Figure 8). The water temperature, EC25 and pH of the surface water flowing over the exposed silt layers were then measured at 10.5-12.7 • C, 388-495 mS/m and 2.3-2.5, respectively, thus showing very high EC25 and low pH. The bulk density of the surface water was then calculated at 1001.6-1001. 9 kg/m 3 , since the EC25 at 388-495 mS/m corresponds to total dissolved solid (TDS) at 1.9-2.5 g/L, using the conversion coefficient 5.00 (mg/L)/(mS/m) of the EC25 to TDS for the portable EC meter of HORIBA, Co., Ltd. Meanwhile, the Okama bottom water on 31 July 2020 (Figure 7) had the bulk density of 1000.4 kg/m 3 for water temperature at 5 • C, and the EC25 was 87 mS/m (equivalent to TDS at 435 mg/L). Thus, the surface water is consistently heavier than the whole lake water. This means that, after flowing into Okama, the surface water flows down the bottom slope as a density underflow [20]. temperature, EC25 and pH of the surface water flowing over the exposed silt layers were then measured at 10.5-12.7 °C, 388-495 mS/m and 2.3-2.5, respectively, thus showing very high EC25 and low pH. The bulk density of the surface water was then calculated at 1001.6-1001. 9 kg/m 3 , since the EC25 at 388-495 mS/m corresponds to total dissolved solid (TDS) at 1.9-2.5 g/L, using the conversion coefficient 5.00 (mg/L)/(mS/m) of the EC25 to TDS for the portable EC meter of HORIBA, Co., Ltd. Meanwhile, the Okama bottom water on 31 July 2020 (Figure 7) had the bulk density of 1000.4 kg/m 3 for water temperature at 5 °C, and the EC25 was 87 mS/m (equivalent to TDS at 435 mg/L). Thus, the surface water is consistently heavier than the whole lake water. This means that, after flowing into Okama, the surface water flows down the bottom slope as a density underflow [20]. Then, the suspended sediment concentration (SSC) at 50-100 mg/L in the surface water could contribute slightly to the magnitude of the bulk density. However, during and just after a heavy rainfall of more than 100 mm/day, the SSC could increase greatly and produce a turbid layer at the bottom of the lake, contributing to the sedimentation. During development of the thermocline, the EC25 increased even more greatly at 0.2 m or less above the lake bottom. This suggests that the deposition of suspended sediment is slow enough to produce a turbid layer, such as the nepheloid layer often observed near the ocean bottom [21]. Figure 9 shows a time series of the hourly water temperature, EC25 and pH from the EC logger at 2 m above the bottom of site MD ( Figure 3) and hourly rainfall at site M for the time period from 15:00, 31 July to 12:00, 19 October 2020. The continual rainfall gradually increased the EC25 and decreased the pH. This appears to reflect the integrated stagnation of the high EC25 and low pH water from the intrusion of density underflow into the lower layer (Figure 7b,d). During this intrusion of the underflow, mixing with lake water could occur. Thus, the EC25 at the logger may be lower than that of the surface water. Meanwhile, the water temperature gradually increased with slight decreases of 0.03-0.05 °C, being highly responsive to each rainfall. Then, the suspended sediment concentration (SSC) at 50-100 mg/L in the surface water could contribute slightly to the magnitude of the bulk density. However, during and just after a heavy rainfall of more than 100 mm/day, the SSC could increase greatly and produce a turbid layer at the bottom of the lake, contributing to the sedimentation. During development of the thermocline, the EC25 increased even more greatly at 0.2 m or less above the lake bottom. This suggests that the deposition of suspended sediment is slow enough to produce a turbid layer, such as the nepheloid layer often observed near the ocean bottom [21]. Figure 9 shows a time series of the hourly water temperature, EC25 and pH from the EC logger at 2 m above the bottom of site MD ( Figure 3) and hourly rainfall at site M for the time period from 15:00, 31 July to 12:00, 19 October 2020. The continual rainfall gradually increased the EC25 and decreased the pH. This appears to reflect the integrated stagnation of the high EC25 and low pH water from the intrusion of density underflow into the lower layer (Figure 7b,d). During this intrusion of the underflow, mixing with lake water could occur. Thus, the EC25 at the logger may be lower than that of the surface water. Meanwhile, the water temperature gradually increased with slight decreases of 0.03-0.05 • C, being highly responsive to each rainfall.   Figure 10d is an enlargement of the red arrow period in Figure 10c. The five black arrows in Figure 10a-c present the measurement times of the vertical profiles in Figures 6 and 7. In mid-October 2019 and late October 2020, the water temperature increased abruptly at the bottom and 2 m above the bottom, which was preceded by an increase at 10 m above the bottom at or around the end of September (Figure 10a). This indicates the descent and bottom arrival of the thermocline during the cooling season (Figures 6 and 7). At the bottom and 2 m above the bottom, the temperature appears to slightly increase throughout the heating and cooling seasons before this abrupt temperature increase ( Figure 9).

Thermal Variations in Okama
After the abrupt temperature increase, Okama was thermally uniform, as shown by the 18 October profile in Figure 6a. Thereafter, the lake surface cooled to 0 °C and remained completely ice-covered until early May 2020 or early April 2021, as shown by the almost 0 °C temperature at the lake surface (Figure 10c,d). In early January of 2020, the surface temperature became less than 0 °C. This indicates that, following the downward ice growth, the temperature logger just below the surface buoy was enclosed by lake ice. For the completely ice-covered period of 20 December 2019-1 May 2020, a small oscillation of the water temperature was observed at the bottom and 2 m above the bottom. This is probably due to the small vertical movement of the covered ice responding to air pressure variation, since both the surface buoy and the lower buoy in the mooring system ( Figure  3) were captured, possibly by the covered ice grown.
Such an ice capture of the two buoys could occur via the consistent decrease in lake water level; in early November, the lake level may have started to decrease due to a lack of rainwater supply ( Figure 11). In the completely ice-covered period of December 2020-March 2021, the water temperature at 5 m depth was 0.03-0.09 °C larger than that at 15 m above the bottom (Figure 10b). This means that due to the consistent decrease in the lake level, the logger at 5 m depth had come down to the level as low as the logger at 15 m   Figure 10d is an enlargement of the red arrow period in Figure 10c. The five black arrows in Figure 10a-c present the measurement times of the vertical profiles in Figures 6 and 7. In mid-October 2019 and late October 2020, the water temperature increased abruptly at the bottom and 2 m above the bottom, which was preceded by an increase at 10 m above the bottom at or around the end of September (Figure 10a). This indicates the descent and bottom arrival of the thermocline during the cooling season (Figures 6 and 7). At the bottom and 2 m above the bottom, the temperature appears to slightly increase throughout the heating and cooling seasons before this abrupt temperature increase (Figure 9).

Thermal Variations in Okama
After the abrupt temperature increase, Okama was thermally uniform, as shown by the 18 October profile in Figure 6a. Thereafter, the lake surface cooled to 0 • C and remained completely ice-covered until early May 2020 or early April 2021, as shown by the almost 0 • C temperature at the lake surface (Figure 10c,d). In early January of 2020, the surface temperature became less than 0 • C. This indicates that, following the downward ice growth, the temperature logger just below the surface buoy was enclosed by lake ice. For the completely ice-covered period of 20 December 2019-1 May 2020, a small oscillation of the water temperature was observed at the bottom and 2 m above the bottom. This is probably due to the small vertical movement of the covered ice responding to air pressure variation, since both the surface buoy and the lower buoy in the mooring system ( Figure 3) were captured, possibly by the covered ice grown.
Such an ice capture of the two buoys could occur via the consistent decrease in lake water level; in early November, the lake level may have started to decrease due to a lack of rainwater supply ( Figure 11). In the completely ice-covered period of December 2020-March 2021, the water temperature at 5 m depth was 0.03-0.09 • C larger than that at 15 m above the bottom (Figure 10b). This means that due to the consistent decrease in the lake level, the logger at 5 m depth had come down to the level as low as the logger at 15 m above the lake bottom ( Figure 3). In this situation, the lower buoy could float at the same level as the surface buoy and be captured by the covered ice. On 26 June 2020 and 18 May 2021, when the first survey in the open water season of Okama was performed, the mooring system was found to have moved to a shallower point 66.5 m northwest and 81.7 m southeast of site MD, respectively. This suggests that the mooring system moved due to a wind-driven drift of the covered ice during May (Figure 10c). The covered-ice drift could have started with the partial ice loss at the lake surface and continued until the abrupt increase in surface water temperature (Figure 10c,d). From 29 March to 26 April 2021, the water temperature at the surface and 1 m depth increased from 0.12 to 0.66 • C and from 0.19 to 1.70 • C at 1 m depth, respectively (Figure 10d). This probably reflects the passage of solar radiation onto the lake water through ice cracks during the season of snowmelt and ice melt.
Hydrology 2022, 9, x FOR PEER REVIEW 12 of 19 above the lake bottom ( Figure 3). In this situation, the lower buoy could float at the same level as the surface buoy and be captured by the covered ice. On 26 June 2020 and 18 May 2021, when the first survey in the open water season of Okama was performed, the mooring system was found to have moved to a shallower point 66.5 m northwest and 81.7 m southeast of site MD, respectively. This suggests that the mooring system moved due to a wind-driven drift of the covered ice during May (Figure 10c). The covered-ice drift could have started with the partial ice loss at the lake surface and continued until the abrupt increase in surface water temperature (Figure 10c,d). From 29 March to 26 April 2021, the water temperature at the surface and 1 m depth increased from 0.12 to 0.66 °C and from 0.19 to 1.70 °C at 1 m depth, respectively (Figure 10d). This probably reflects the passage of solar radiation onto the lake water through ice cracks during the season of snowmelt and ice melt. During the complete ice cover, the temperature gradually increased between the lake bottom and 15 m above the bottom (Figure 10a,b). At the surface and 1 m depth, such a gradual increase was not observed (Figure 10d), and the temperature at 15 m above the bottom was almost equal to that at 5 m depth (Figure 10b). Hence, the upper limit of the gradual increase could be located at ca. 17 m above the lake bottom, equivalent to ca. 3 m depth. The lake water is at almost 0 °C at the interface of water and ice, and the downward ice growth tends to cool the lower water [22]. Hence, the temperature increase at 0-15 m above the bottom is probably due to the geothermal heat supplied at the lake bottom, which could originate from the hydrothermal reservoir below the Okama bottom basin [23]. Figure 11 shows temporal variations of the daily mean air temperature, lake surface water temperature and lake level for 1 August-30 November 2020; daily mean relative humidity and wind speed for 1 August-18 October 2020; and the diurnal rainfall for 1 July-18 October 2020. The air pressure and air temperature during 19 October-30 November 2020 at site M, after the removal of the weather station, were evaluated by the regression lines with high correlations (R 2 = 0.982 and 0.991, respectively) between site M and the Daikokuten observation hut ( Figure 12). Figure 11. Temporal variations of the daily mean air temperature, lake level and lake surface water temperature for 1 August-30 November 2020; daily mean relative humidity and wind speed for 1 August-18 October 2020; and the diurnal rainfall for 1 July-18 October 2020. The daily mean air temperature and air pressure (at lake level) for 19 October-30 November were evaluated using the high correlations for the data between site M and the Daikokuten observation hut (Figure 12). Figure 11. Temporal variations of the daily mean air temperature, lake level and lake surface water temperature for 1 August-30 November 2020; daily mean relative humidity and wind speed for 1 August-18 October 2020; and the diurnal rainfall for 1 July-18 October 2020. The daily mean air temperature and air pressure (at lake level) for 19 October-30 November were evaluated using the high correlations for the data between site M and the Daikokuten observation hut (Figure 12). Figure 11. Temporal variations of the daily mean air temperature, lake level and lake surface water temperature for 1 August-30 November 2020; daily mean relative humidity and wind speed for 1 August-18 October 2020; and the diurnal rainfall for 1 July-18 October 2020. The daily mean air temperature and air pressure (at lake level) for 19 October-30 November were evaluated using the high correlations for the data between site M and the Daikokuten observation hut (Figure 12). The lake level declined during a non-rainfall period from 18 to 28 August 2020, as well as after the removal of the weather station at site M. Referring to Equation (1), this suggests that during a relatively lengthy non-rainfall period, the groundwater outflow from the lake is larger than the groundwater inflow plus the river inflow minus the evaporation over the lake surface, i.e., Gout > Gin + Rin-EA0.

Meteorology at Okama
In the period from 18 to 28 August 2020, the lake level decreased at an average of 33 mm/day. Meanwhile, using Equations (2) and (3), the mean evaporation at the lake water During the complete ice cover, the temperature gradually increased between the lake bottom and 15 m above the bottom (Figure 10a,b). At the surface and 1 m depth, such a gradual increase was not observed (Figure 10d), and the temperature at 15 m above the bottom was almost equal to that at 5 m depth (Figure 10b). Hence, the upper limit of the gradual increase could be located at ca. 17 m above the lake bottom, equivalent to ca. 3 m depth. The lake water is at almost 0 • C at the interface of water and ice, and the downward ice growth tends to cool the lower water [22]. Hence, the temperature increase at 0-15 m above the bottom is probably due to the geothermal heat supplied at the lake bottom, which could originate from the hydrothermal reservoir below the Okama bottom basin [23]. Figure 11 shows temporal variations of the daily mean air temperature, lake surface water temperature and lake level for 1 August-30 November 2020; daily mean relative humidity and wind speed for 1 August-18 October 2020; and the diurnal rainfall for 1 July-18 October 2020. The air pressure and air temperature during 19 October-30 November 2020 at site M, after the removal of the weather station, were evaluated by the regression lines with high correlations (R 2 = 0.982 and 0.991, respectively) between site M and the Daikokuten observation hut ( Figure 12).

Meteorology at Okama
The lake level declined during a non-rainfall period from 18 to 28 August 2020, as well as after the removal of the weather station at site M. Referring to Equation (1), this suggests that during a relatively lengthy non-rainfall period, the groundwater outflow from the lake is larger than the groundwater inflow plus the river inflow minus the evaporation over the lake surface, i.e., G out > G in + R in − EA 0 .
In the period from 18 to 28 August 2020, the lake level decreased at an average of 33 mm/day. Meanwhile, using Equations (2) and (3), the mean evaporation at the lake water surface was calculated to be 0.51 mm/day, while the river inflow was 0.0033 m 3 /s, equivalent to 3.5 mm/day per lake surface area. Hence, the net groundwater inflow decreases the lake level at −36 mm/day. On 24-25 September 2020, a large rainfall occurred totaling 231.6 mm. Thereby, the lake level rose 0.47 m. A decline of the lake level was not observed during rainfalls of more than 30 mm/day. This means that, due to negligibly small evaporation during rainfalls, the absolute value of the net groundwater inflow is equivalent to the 30 mm/day rainfall over lake surface plus the correspondent river inflow. After the removal of the weather station at site M, the lake level consistently decreased at 27 mm/day, reflecting, small rainfalls, if any. As a side note, the monthly rainfall at a weather station 20.3 km east of Okama was only 2.0 mm in November 2020.
The data of meteorology and lake level in Figure 11 can be applied to the hydrological budget estimate for the lake, shown by Equation (1).
The surface water inflow with high EC25, low pH and milky color in Figure 8 appears to have been produced continuously by continual rainfalls of 14-48 mm/day during 1-27 July 2020, in addition to the large rainfall of 169 mm/day on 28 July 2020.

Geothermal Heat Flux in Okama
In the completely ice-covered periods of December-April, an increase in the lake water temperature was observed between the lake bottom and 15 m above the bottom at the deepest point ( Figure 10). This means that some geothermal heat is supplied at the lake bottom. Hence, using the bathymetric map in Figure 2 and a time series of the daily mean water temperature from Figure 10, a time series of the daily mean heat flux and its 10-day moving average were acquired ( Figure 13). Then, the periods of 17 December 2019-27 April 2020 and 15 December 2020-29 March 2021 in Figure 10 were selected for the calculation, since the water temperature at the surface increased after 27 April 2020 and 29 March 2021, indicating the incidence of solar radiation via a partial loss of covered ice from spring ice melt.

Estimating Groundwater Inflow and Outflow in Okama
Using Equations (1)- (7), the groundwater inflow and outflow in Okama were calculated in the daily database, where one day was adopted as the fundamental budget time period Δt in Equations (1) and (5). Here, the SO4 2-concentration was taken as the representative ionic concentration for the chemical budget of Okama because the SO4 2-concentration highly correlated to the EC25 and pH ( Figure 5). Then, the net depositional flux S in Equation (5) was assumed to be neglected for SO4 2-. The EC25 record of the EC logger at 2 m above the bottom of site MD was converted to a daily series of the SO4 2-concentration, where the daily series of the volume-averaged SO4 2-concentration as CL in Equation (5) was obtained using a relation between the EC25 values in the logger and the three vertical profiles in Figure 7b (Figure 14). Then, assuming a horizontal multi-layered structure of the chemistry in the lake, the volume-averaged EC25 values from the three profiles (vertical axis in Figure 14) were calculated using the bathymetric map in Figure 2. Additionally, the water surface area A0 was temporally changed following the lake level change in Figure 11 and the bathymetry in Figure 2. The total stream inflow measured on the nonrainfall days was given as the river inflow Rin in Equations (1) and (5). The mean EC25 and SO4 2-concentration values (177 mS/m and 823 mg/L, respectively) of pore water in the bottom sediment at site MD were adopted as those of the groundwater inflow. First, the area, A(z), at each logger's height above the lake bottom, acquired by the bathymetric map, was multiplied by the following factor: where ρ w is the water density (kg/m 3 ), c p is the specific heat (J/Kg/K) of the water, and T(z) is the water temperature (K) at the logger's height z(m) above the bottom. The total heat storage (J) of the lake on a certain day was obtained by summing the following: at 0-17 m above the lake bottom, where ∆h is each height difference between the loggers. Then, a horizontal multi-layer thermal structure in the lake was assumed, i.e., horizontally taking the same temperature at the same water depth. Finally, the heat storage change per day per unit area (W/m 2 ) was calculated by dividing a daily change (J/day) in the total heat storage of the lake by the area (6.627 × 10 4 m 2 ) at 17 m above the bottom (Figure 2).
Focusing on the temporal variations of the 10-day moving average, the heat flux varied at a range of −0.5-5.6 W/m 2 in December 2019-April 2020 and 1.5-4.7 W/m 2 in December 2020-March 2021. The heat flux averages over these periods were 2.5 W/m 2 in December 2019-April 2020 and 2.9 W/m 2 in December 2020-March 2021. These values are similar to the heat flux (2.8-10.7 W/m 2 ) at Lake Kuttara belonging to the active Kuttara Volcano [9] but larger than that (0.29 W/m 2 ) at Lake Shikotsu, which is surrounded by three active volcanos [24]. By using the audio-frequency magnetotelluric (AMT) method, Ichiki et al. [23] revealed that a low resistivity zone, such as a hydrothermal reservoir, exists with the center at ca. 250 m below the lake bottom. The calculated geothermal heat flux in Figure 13 appears to reflect thermal leakage from the hydrothermal reservoir.

Estimating Groundwater Inflow and Outflow in Okama
Using Equations (1)- (7), the groundwater inflow and outflow in Okama were calculated in the daily database, where one day was adopted as the fundamental budget time period ∆t in Equations (1) and (5). Here, the SO 4 2− concentration was taken as the representative ionic concentration for the chemical budget of Okama because the SO 4 2− concentration highly correlated to the EC25 and pH ( Figure 5). Then, the net depositional flux S in Equation (5) was assumed to be neglected for SO 4 2− . The EC25 record of the EC logger at 2 m above the bottom of site MD was converted to a daily series of the SO 4 2− concentration, where the daily series of the volume-averaged SO 4 2− concentration as C L in Equation (5) was obtained using a relation between the EC25 values in the logger and the three vertical profiles in Figure 7b (Figure 14). Then, assuming a horizontal multi-layered structure of the chemistry in the lake, the volume-averaged EC25 values from the three profiles (vertical axis in Figure 14) were calculated using the bathymetric map in Figure 2. Additionally, the water surface area A 0 was temporally changed following the lake level change in Figure 11 and the bathymetry in Figure 2. The total stream inflow measured on the non-rainfall days was given as the river inflow R in in Equations (1) and (5). The mean EC25 and SO 4 2− concentration values (177 mS/m and 823 mg/L, respectively) of pore water in the bottom sediment at site MD were adopted as those of the groundwater inflow.  Figure 14. Relation between the volume-averaged EC25 from the three profiles in Figure 7b and the EC25 from the EC logger at 2 m above the bottom of site MD in 2020. Table 1 shows the net groundwater inflow (Gin-Gout), as well as the separated groundwater inflow Gin and groundwater outflow Gout, calculated for the six non-rainfall periods in Figure 11. The correspondent river flow observed is also shown. For the non-rainfall periods, the first day after a rainfall was omitted because the water supplied by the rainfall appears to raise the lake level one day later ( Figure 11). All the calculated results of the net groundwater inflow are negative. Thus, the groundwater outflow is larger than the groundwater inflow. Of all the calculated results, Periods (1), (2) and (6) include more than three non-rainfall days. The calculated results for these three periods could thus be more reliable. Gin-Gout, Gin and Gout averaged for the three periods were-0.028 m 3 /s, 0.012 m 3 /s and 0.039 m 3 /s, respectively. This negative net groundwater inflow decreases the lake level at-0.029 m/day, and thus is comparable in magnitude to the observed declining rates Figure 14. Relation between the volume-averaged EC25 from the three profiles in Figure 7b and the EC25 from the EC logger at 2 m above the bottom of site MD in 2020. Table 1 shows the net groundwater inflow (G in -G out ), as well as the separated groundwater inflow G in and groundwater outflow G out , calculated for the six non-rainfall periods in Figure 11. The correspondent river flow observed is also shown. For the non-rainfall periods, the first day after a rainfall was omitted because the water supplied by the rainfall appears to raise the lake level one day later ( Figure 11). All the calculated results of the net groundwater inflow are negative. Thus, the groundwater outflow is larger than the groundwater inflow. Of all the calculated results, Periods (1), (2) and (6) include more than three non-rainfall days. The calculated results for these three periods could thus be more reliable. G in − G out , G in and G out averaged for the three periods were −0.028 m 3 /s, 0.012 m 3 /s and 0.039 m 3 /s, respectively. This negative net groundwater inflow decreases the lake level at −0.029 m/day, and thus is comparable in magnitude to the observed declining rates of −0.033 m/day and −0.027 m/day ( Figure 11). Table 1. Measured river inflow R in ; calculated net groundwater inflow (G in − G out ); as well as separated groundwater inflow G in and outflow G out for the six non-rainfall periods in Figure 11.

No.
Period R in (m 3 /s) G in − G out (m 3 /s) G in (m 3 /s) G out (m 3 /s) If the heat flux H Gin (W/m 2 ) by the groundwater inflow G in at 0.012 m 3 /s is equal to the geothermal heat flux, 2.5 W/m 2 in 2019-2020 and 2.9 W/m 2 in 2020-2021 at the lake bottom, the inflowing groundwater temperature T Gin is estimated using the following equation: H Gin = ρ w c p G in (T Gin − T L ) (10) where ρ w is the groundwater density, c p is the specific heat (J/kg/K) of water and T L is the representative lake water temperature (K) (273.15 K = 0 • C + 273. 15). Here, if T L is supposed to be the bottom water temperature 275.45 K (2.3 • C) and 275.25 K (2.1 • C) averaged over the periods of 2019-2020 and 2020-2021, respectively, T Gin is given at 278.75 K (5.6 • C) and 279.05 K (5.9 • C), respectively. These values are larger than the annual mean air temperature 277.95 K (4.8 • C) at Okama inferred from those at Daikokuten. In general, the shallow groundwater in the catchment without any geothermal heat source tends to take the water temperature near to the local, annual mean air temperature [25,26]. The geothermal heat flux in Okama is thus significantly larger than that in the non-geothermal catchment. This suggests that, though Okama stays relatively calm, thermal leakage from the hydrothermal reservoir below the lake basin exists. In the completely ice-covered periods, the infiltration of snowmelt water at the snow-ground interface should be much smaller than that of rainwater in the rainfall season [27,28]. This means that the groundwater inflow to Okama in the ice-covered period is smaller than that in the rainfall season. For example, if G in = 0.010 m 3 /s at the calculated H Gin = 2.5 and 2.9 W/m 2 , Equation (10) gives T Gin = 279.45 K (6.3 • C) and 279.85 K (6.7 • C), respectively. Hence, the calculated temperatures 5.6 • C and 5.9 • C appear to correspond to the lowest values in the ice covered period.

Conclusions
In order to explore the geothermal effect on Okama Crater Lake in the active Zao Volcano, as well as the hydrological and chemical cycles in or around the lake related to the underground geothermal heat source, the lake water temperature was monitored at eight depths throughout the year and the hydrological and chemical budgets of the lake were estimated. The temperature monitoring revealed that, during the completely ice-covered period, the water temperature consistently increases between the lake bottom and 15 m above the bottom. This temperature increase is probably due to geothermal leakage from the hydrothermal reservoir below the lake bottom. The geothermal heat flux was calculated at 2.5 W/m 2 and 2.9 W/m 2 during the completely ice-covered seasons of 2019-2020 and 2020-2021, respectively.
Meanwhile, estimates of the hydrological and chemical budgets for Okama provided the groundwater inflow at 0.012 m 3 /s and the groundwater outflow at 0.039 m 3 /s, averaged over the three non-rainfall periods of 2020. In the chemical budget estimate, the SO 4 2− concentration was selected because SO 4 2− was dominant in the water chemistry and highly correlated with the electric conductivity at 25 • C (EC25) as well as the pH of lake water and stream water. The negative groundwater inflow at −0.028 m 3 /s can produce a decline in the lake level of −0.029 m/day during non-rainfall periods.
A coupling of the geothermal heat flux and groundwater inflow, which was obtained by two independent methods, evaluated the temperature of inflowing groundwater. Then, it was assumed that the geothermal heat flux calculated at 2.5 W/m 2 or 2.9 W/m 2 is supplied by the groundwater inflow estimated at 0.012 m 3 /s. The groundwater temperature was calculated at 5.6 • C or 5.9 • C. These values are significantly larger than the annual mean air temperature (4.8 • C) at Okama. In the catchment without a geothermal heat source, the shallow groundwater is supposed to take a temperature near to the local annual mean air temperature. Hence, geothermal leakage from the underground hydrothermal reservoir is judged to be significant in the lake, and, therefore, it is necessary to continue monitoring Okama.
As a next step, an estimate of the thermal budget in the lake, including the heat flux at lake surface, is needed to determine whether the geothermal heat flux obtained by this study is reasonable or not in magnitude, as well as to understand how the inflowing and outflowing groundwaters are related to the hydrothermal flow system below Okama.