Vertical mixing and oxygen flux caused by daily sea breezes in a shallow stratified lake

Vertical transport caused by mixing is essential for understanding physical processes in lakes. However, mixing processes in shallow lakes are not well understood because of the lack of turbulence measurements. This study presents observations of vertical mixing and oxygen flux in a shallow lake, Lake Kitaura of the Lake Kasumigaura continuous lake system, which is located along the central eastern coast of the Japanese mainland. Mooring and microstructure surveys were conducted in August 2020. The vertical eddy diffusivity was estimated from the Ellison scale using high-frequency sampled temperature data from the mooring location, and the estimations were consistent with those observed from the microstructure profiler. The estimated eddy diffusivity revealed a mixing structure and oxygen flux in the lake during the study period. The daily cycle of stratification and mixing was caused by daily heating due to solar radiation and by winds from daily sea breezes, respectively. The daily stratification maximum occurred around noon, which suppressed vertical mixing. Vertical mixing was intensified due to sea breezes in the afternoon, which led to movement of oxygen from the surface layer to the bottom layer. The maximum vertical mixing was observed at 18:00. The oxygen concentration did not increase during nighttime when the negative surface heat flux was observed; thus, nighttime cooling may not primarily contribute to the vertical oxygen supply.


Introduction
Vertical mixing plays a crucial role in physical processes of aquatic ecosystems, such as oceans (Munk and Wunsch 1998), estuaries and lakes (Wüest and Lorke 2003). In particular, the vertical flux and distribution of water constituents, such as nutrients and dissolved oxygen, are important for understanding and managing water quality. Vertical fluxes are primarily controlled by vertical eddy diffusivity and mixing; thus, numerous studies have focused on mixing-induced vertical flux in lakes (e.g., Rao et al. 2008;Yamazaki et al. 2009).
One of the major water quality issues in lakes, reservoirs, and estuaries is hypoxia induced by suppressed vertical oxygen flux from the surface layer toward the bottom layer (Breitburg 2002). Therefore, water quality management that considers vertical mixing and the oxygen supply is crucial for the sustainable usage of freshwater. Hypoxic water appears in lakes and shallow estuaries worldwide (Diaz, 2001) and is often defined as water with a dissolved oxygen (DO) concentration less than 2-3 mg L −1 (Lee and Lwiza 2008;Masunaga and Komuro 2020). A recent study by Woolway et al. (2021) showed that enhanced stratification due to global warming accelerates deoxygenation in lakes. Hypoxia triggers the release of nutrients and metals from sediments, and this has been recognized as one of the major issues in freshwater management. For example, nitrogen is released from sediments into lake water when the DO concentration is less than 1 mg L −1 (Zhang et al. 2014), phosphorus is released from lake sediments under hypoxic conditions (Wang et al. 2008), and heavy metals (Pb, Zn, and Fe) (Kang et al. 2019). However, details of oxygen supply mechanisms and associated processes have not been well-studied because of the lack of mixing observations in shallow waters. The mixing intensity of water bodies is often analyzed using microstructure shear measurements from a free-fall microstructure profiler, and this is the most popular and reliable method. In particular, estimates of the turbulent kinetic energy dissipation rate and vertical eddy diffusivity using a free-fall microstructure profiler are useful to evaluate turbulent transport because they directly consider the microstructure shearing caused by turbulent eddies. Numerous studies have employed free-falling microstructure profilers to observe vertical mixing in oceans and lakes (e.g., Auger et al. 2013;Masunaga et al. 2016;Itoh et al. 2021). However, it is difficult to mount free-falling microstructure profilers on mooring systems because profilers must be manually operated using a tether cable without vibrations; thus, this method is not suitable for continuous field observations. Time series of mixing parameters can also be observed by high-frequency acoustic current and temperature measurements using an Eulerian mooring observation system (Lorke and Wüest 2005;Ivey et al. 2018). Lorke and Wüest (2005) proposed a method for estimating the rate of turbulent kinetic energy dissipation from an acoustic Doppler current profiler (ADCP) based on the inertial dissipation method. Ivey et al. (2018) developed a new method for inferring the vertical eddy diffusion using high-frequency sampled temperature data from a mooring system. Although typical Eulerian acoustic current and temperature observations do not directly observe microstructure turbulence throughout the full viscosity range, the mixing parameters estimated from these approaches are reasonable (Ivey et al. 2018;Lorke and Wüest 2005). Thus, Eulerian mixing parameterization methods allow us to obtain information on time-series of mixing structures in shallow lakes.

Study area
Lake Kasumigaura, a continuous lake system, is located along the central coast of the Japanese mainland (Fig. 1). The system consists of two lakes, Lake Nishiura (West Lake) and Lake Kitaura (North Lake), the latter of which was the focus of this study. The major axis length and width of the lake are approximately 20 and 4 km, respectively. The averaged depth is approximately 6 m. The lake is connected to the Pacific Ocean through the Tone River ( Fig. 1), but a river barrage located at the southern end of the lake prevents saltwater intrusion from the ocean. Thus, the lake water is completely fresh, and its density is solely determined by temperature.
The summer season is characterized by intense stratification with suppressed vertical mixing in Lake Kasumigaura, resulting in hypoxia in the bottom layer (e.g., Komatsu et al. 2010;Masunaga and Komuro. 2020). A previous study by Masunaga and Komuro (2020) showed that vertical mixing in the lake is explained by the Monin-Obukhov length scale, which states that vertical mixing is primarily controlled by the balance between wind stress and surface heat flux. They also reported that daily-scale DO variation in the bottom boundary layer is controlled by vertical mixing. Although Masunaga and Komuro (2020) reported that mixing and DO are explained by surface heat and momentum fluxes, they did not measure mixing intensity and associate DO transport. To further investigate interior vertical mixing and associated DO transport, we observed the time-series of mixing structures of a shallow lake in the Lake Kasumigaura system using Eulerian ADCP and high-frequency temperature measurements. Fig. 1 a-b The map of the study area in Lake Kasumigaura of Japan. The JMA weather stations in Kashima and Tsukuba are represented by yellow dots in b. c The black star represents the location of the mooring site and the black dots indicate water and sediment sampling locations 1 3

Field observations
Field observations were conducted in Lake Kitaura in summer and fall of 2020. A mooring system was mounted in the middle of the lake (longitude: 140°33.96′ N; latitude: 36°0.043′ E). The location of the mooring is shown in Fig. 1c. As this study focused on the stratified period during summer, only the observed results from August 1 to 31 are presented (31 days in total). The mooring consisted of a 2 kHz Aquadopp current profiler (ADCP; Nortek), TidbiT v2 temperature loggers (Onset Computer Corporation), and an SBE-56 high-resolution temperature logger (Sea-Bird Scientific). The ADCP was mounted on the bottom landing frame and faced upward, and it resolved the water column with a bin size of 0.5 m. The first observed bin obtained from the ADCP was 0.9 m above the lakebed. The observation burst interval of the ADCP was 600 s (10 min), and each burst contained 15 pings and their average value, which were obtained with a sampling frequency of 1 Hz. The TidbiT v2 loggers were mounted along a thermistor chain that was placed from the surface to a depth of 6 m with an interval of 0.6 m. The TidbiT v2 logger was also mounted on the lakebed located at a depth of approximately 6.5 m. Temperature data were recorded at a time interval of 5 min. Four SBE-56 temperature loggers were also mounted along the thermistor chain at depths of 1.2, 2.4, 3.6, and 4.8 m, and their sampling frequencies were set to 1 Hz. The resolutions of the TidbiT v2 and SBE-56 temperature loggers were 0.02 and 1 × 10 -4 °C, respectively. The entire mooring period was between July 30 and November 2, 2020, except for the thermistor chain with TdbiT v2 sensors (June 14-November 2). Hourly DO and wind data were obtained from continuous measurements at the Kamaya Monitoring Station (KMS), which is operated by the Japan Water Agency (JWA) and is located near our mooring site (Fig. 1c). The DO data were measured from the surface to the lakebed with a vertical interval of approximately 0.5 m.
In addition to the mooring observations, we deployed a microstructure profiler, the VMP-250 (Rockland Scientific), to observe the microstructure turbulence of the lake. The VMP-250 carries two shear probes, an FP07 fast-response thermistor, optical fluorescence and turbidity sensors, and standard hydrographic parameter sensors for conductivity, temperature, and depth (CTD). It records data with a sampling frequency of 512 Hz. The conventional free-fall microstructure observation method cannot make measurements near the lake surface (0-2 m) because the stable free-fall mode cannot be established at the beginning of profiling (the stable free-fall is required for estimating correct mixing parameters). To observe the microstructure turbulence near the surface and lakebed, we employed the following two observation methods: (1) a free-fall mode and (2) an uprising mode. In the uprising mode, the VMP-250 with a buoyancy collar was released from the lakebed using a remotecontrolled weight releaser. The VMP-250 was deployed in the following four periods: morning (11:34-11:51) and afternoon (17:24-17:41) on August 5, and morning (10:24-10:39) and afternoon (15:39-15:55) on August 7 (Fig. 2a, d). Five or six free-falling and uprising vertical profiles were collected during each period. In total, we obtained 41 vertical profiles from the microstructure measurements. The average profiling speed was approximately 0.57 and 0.61 m s −1 for the free-falling and rising modes, respectively.
To estimate the heat flux, meteorological parameters were obtained from various weather stations operated by the Japan Meteorological Agency (JMA) near the study site (https:// www. data. jma. go. jp/ obd/ stats/ etrn/). Air temperature data were obtained from the JMA Kashima weather station. Relative humidity, pressure, and shortwave radiation were obtained from the JMA Tsukuba weather station. The cloud cover ratio was obtained from the JMA Tokyo weather station. The locations of the JMA weather stations are shown in Fig. 1b.
We also estimated the oxygen demand of the lake water and sediments using laboratory experiments. Water and Fig. 2 Time series of a temperature and b square of the buoyancy frequency within the entire period of thermistor chain observations (June 14-November 2).The duration between two vertical black dash dotted line show the analysis period of this study sediment samples were collected on July 29 and September 1, 2020, at the mooring site and three other locations in the lake, as shown in Fig. 1c. Water samples were obtained using an electric water pump at a height of 0.5 m above the lakebed. A sediment core sampler, KB Core Sampler (Rigo), was used to collect lake sediment samples. The height and inner diameter of the core sampler were 50 and 7 cm, respectively. Changes in DO due to the water and sediment oxygen demands (WOD and SOD) were measured using sealed acrylic cylinder tanks under dark conditions with a fiber optic oxygen meter (Fibox3, PreSens). For the SOD experiments, sediment samples in the cylinder tanks were covered with filtered lake water without suspended solids (a spool-type cartridge filter with a filter length of 1.0 mm was used to remove suspended solids). The initial DO was set at the saturation level (approximately 8.0 mg L −1 ) and temperature was maintained at 25 °C during the laboratory experiments to reproduce a typical condition in Lake Kitaura in the summer season. The oxygen demand was measured for two water samples and two sediment samples at each location on each date; therefore, we conducted a total of 16 experiments each for WOD and SOD. For more details on the laboratory experiments, refer to Ishii et al. (2009).

Vertical eddy diffusivity estimations
The vertical eddy diffusivity was estimated using three independent methods: (1) a conventional direct microstructure observation method using the VMP-250, (2) a recently developed method proposed by Ivey et al. (2018) using the Ellison scale (L E ) and the vertical velocity gradient (S v ), and (3) a method based on the relationship between L E and the Ozmidov length scale (L O ).
The first method is a commonly used estimation method proposed by Osborn (1980) as follows: where K ρ is the vertical eddy diffusivity, Γ is the mixing efficiency, ε is the rate of turbulent kinetic energy dissipation, and Ν is the buoyancy frequency given as follows: where g is the gravitational acceleration, ρ is the water density, ρ 0 is the reference density of 1000 kg m −3 and z is the vertical coordinate. The mixing efficiency was assumed to be Γ = 0.2, a commonly used constant (Osborn, 1980). The rate of turbulent kinetic energy dissipation is given by the following isotropic formula: where u' is turbulent velocity fluctuations observed by the shear probe and ν is the kinematic viscosity. ε was estimated by fitting the universal theoretical Nasmyth spectrum (Nasmyth 1970) with a total of 512 data points (bin size, 1 s).
The second method is to use an approach recently developed by Ivey et al. (2018), as described by the following equations: where T' is the turbulent fluctuation in temperature, T � 1∕2 represents the root-mean square of T', dT∕dz is the background vertical temperature gradient, and u and v are the zonal and meridional velocity components, respectively. The subscript of I18 in Eq. (4) indicates the method proposed by Ivey et al. (2018). Note that temperature, T, in Eq. (5) has to be density in oceans influenced by both of salinity and temperature. However, our study site was a freshwater area, so it was possible to use only temperature to compute L E . L E was obtained from the moored thermistor chain, and this method allowed us to continuously observe the eddy diffusivity during the entire mooring observation period. In this study, L E was computed from the high-frequency sampled temperature data from the SBE-56 and the background vertical temperature gradient from the TidbiT temperature loggers. Background low-frequency oscillations in temperature data must be removed before calculating L E (Cimatoribus et al. 2014). In our observations, the buoyancy frequency (N/2π) reached an order of O (10 -2 ) cps. In order to remove background oscillations with frequencies lower than the buoyancy frequency, the temperature fluctuations, T', were high-pass filtered observed temperatures with a cutoff frequency of 6.67 × 10 -2 cps (cycle of 15 s). The temperature gradient was obtained from temperature data at 0.6 m above and below the SBE 56 (i.e., at dz = 1.2 m). Then, L E was computed within a time interval of 5 min. As the resolution of the TidbiT temperature logger is 0.02 °C, the noise level is assumed to be 1.67 × 10 -2 °C m −1 . Any temperature gradient data below this noise level were replaced with the noise , level value. The L E was also measured using the VMP-250 microstructure profiler. We used high-pass filtered temperature data obtained from the FP07 thermistor with a cutoff frequency of 2 cps (approximately 3.3 cpm) to compute L E from the VMP-250 data. To remove noise-like high-frequency fluctuations in the velocity data, the velocity data obtained from the ADCP were averaged with a 4-h running mean.
The third method is based on the relationships between L E , L O , and the Thorpe scale (L T ). L O is defined as the largest vertical length of a vertical eddy unaffected by buoyancy, as described by Ozmidov (1965)as follows: where ε is the kinetic energy dissipation rate and N is the buoyancy frequency. L T is the root mean square of the Thorpe displacement profile (Thorpe 1977), and it is known to be equivalent to L E (Mater et al. 2013). L T is positively correlated with L O , and the L T -L O relationship has been frequently used to infer ε from a conventional CTD profile, as follows: where c is an empirical constant equal to L O /L T . Several values of c have been reported in previous studies, e.g., c = 0.79 (Dillon 1982), 0.66 (Crawford 1986), 1.50 (Peters 1997), and 0.90 (Ferron et al. 1998). Wesson and Gregg (1994) showed that L O /L T ranges from approximately 0.25-4. In this study, we assumed that L O is a linear function of L E on a logarithmic scale; that is, where a and b are empirical constants measured from microstructure observations. Considering this, the vertical eddy diffusivity for the third method can be written as follows: where Γ is the mixing efficiency (0.2), and Ν is the buoyancy frequency. Although this method requires constants a and b to be estimated from a microstructure profiler, velocity measurements using ADCP are not required.

Background conditions
Time series of temperature and N 2 during the entire mooring period are shown in Fig. 2. In our study area, August was the warmest and the most stratified month in the year. The temperature gradually increased throughout the study period (August 2020) ranging from 24.8 to 32.8 °C ( Fig. 2a and 3a).
The black chain-dotted lines in Fig. 3a, b show the mixed layered depth, which is defined as the depth at which the temperature differs from the lake surface by 0.5 °C (Levitus 1982). Strongly stratified periods (N 2 > 1 × 10 -3 rad 2 s −2 ) intermittently appeared, and the N 2 reached a maximum of 1.3 × 10 -2 rad 2 s −2 on August 16 (Fig. 3b). To evaluate the water column stability, we estimated the Richardson number Ri, which is determined by the following: It is known that the condition of Ri < 0.25 is necessary for instability of the background flow. The estimated Ri was lower than 0.25 during mixing periods (N 2 ~ O(10 -4 ) rad 2 s −2 ) and higher than 0.25 during stratified periods (N 2 > O(10 -3 ) rad 2 s −2 ) (Fig. 3c, d). Therefore, Ri provides a good representation of the vertical mixing conditions in the study area.
Atmospheric forcing, surface heat flux, wind speed, wind direction, and the daily averaged nondimensionalized Monin-Obukhov length scale (I MO ) are presented in Fig. 4. The parameter I MO indicates the stability of the water column, which is explained by the ratio between the surface wind stress and the buoyancy flux, as follows: where h is the height of the water column (depth), L MO is the Monin-Obukhov length scale, u w* is the wind frictional velocity, expressed as u w* = (ρ a C Dw U w 2 ) 1/2 ; where ρ a is the air density, 1.2 kg m −3 ; C Dw is the wind drag coefficient, 1.5 × 10 -3 ; and U w is the observed wind speed at 10 m above the lake surface), k is the von Karman constant (0.41), and J is the surface buoyancy flux. The surface buoyancy flux was calculated as follows: where e is the thermal expansion coefficient, g is the gravitational acceleration, H surf is the lake surface heat flux (the sum of shortwave and longwave radiation, and sensible and latent heat fluxes), ρ 0 is the reference density of 1000 kg m −3 , and C P is the specific heat of water (4184 J kg −1 K −1 ). The shortwave radiation was obtained from the Tsukuba JMA weather station. The long-wave radiation and sensible and latent heat fluxes were estimated using conventional parameterizations (Rayson et al. 2015;Masunaga and Komuro 2020). When the I MO is higher than unity, vertical mixing throughout the water column is expected. The surface heat flux showed a clear daily intensified heat flux toward the lake due to daytime shortwave radiation (Fig. 4a). The wind ranged from approximately 0 to 10 m s −1 in variable directions (Fig. 4b,  c). The daily averaged I MO , or buoyancy flux, was less than zero for seven out of the 31 days in August 2020. When I MO was positive, it ranged from 3.7 × 10 -1 to 2.7 × 10 2 . In approximately one-third of the study period (11 days) was occupied by 0 < I MO < 1. Therefore, both stable stratified and unstable mixing periods were observed in our observational data.

Vertical eddy diffusivity
The observational periods of the microstructure survey using the VMP-250 were categorized into the following two regimes: (1) the stratified period on August 5 (Ri = O(1) and I MO < 1), and (2) the mixed period on August 7 (Ri = O(0.1) and I MO ~ 1) (Fig. 3a-d). The mixing on August 7 was induced by intensified winds and weakened surface heat flux (Fig. 4a, b). In the first VMP-250 observation (morning of August 5), ε reached approximately 2.0 × 10 -6 W kg −1 , and the dissipation was much weaker in the deeper layers than in the surface layer (Fig. 5a). The vertical eddy diffusivity, K ρ , was weak in the deep layer during this period, and it was similar to the molecular diffusivity of heat, O(10 -7 ) m 2 s −1  . 5i). In the afternoon of August 5, vertical mixing was intensified by wind forcing (U w ~ 6 m s −1 ), and K ρ was higher than 1.0 × 10 -5 m 2 s −1 in the whole water column (Fig. 5j). The stratification near the surface was drastically reduced from the morning to the afternoon on August 5, and N 2 was reduced by a factor of 100 (Fig. 5e, f). In the observations on August 7, ε and K ρ appeared to be similar in the morning and afternoon (Fig. 5c, d, k, l). Similar to the change in stratification observed on August 5, N 2 was reduced near the surface in the afternoon (Fig. 4e-h). Meanwhile, the stratification was intensified near the bottom on both August 5 and 7 ( Fig. 5e-h, depth ~ 5-6 m), which was caused by a deepened thermocline in the afternoon. As shown in the results above, mixing conditions drastically changed within several hours in the study area, especially between the morning and afternoon on August 5. The observed ε was compared with the law-of-the-wall parameterization given as follows: where D is the depth (the distance from the lake surface). ε low is plotted as a blue solid line in Fig. 4a-d. The law-ofthe-wall estimation was close to the observed ε on August 5, but it overestimated the dissipation on August 7. The correlation coefficient was estimated at 0.69; however, the averaged ratio between ε low and ε was approximately 34, viz., ε low was roughly O(10) times higher than e. This overestimation may have been caused by buoyancy effects that are not considered in the simple law-of-the-wall parameterization. Therefore, this simple law does not work as a quantitative analysis method for inferring ε and K ρ .
The Ozmidov length scale (L O ) and Ellison length scale (L E ) obtained from the VMP-250 are compared in Fig. 6a. Coefficients a and b in Eq. (9) were computed from this comparison, and they were estimated at 5.2 and 1.5, respectively (i.e., L O = 5.2L E 1.5 ). It has been reported that the L O /L E ratio becomes higher than unity under well-mixed conditions, and the turbulent stage becomes "mature" (Smyth et al. 2001;Garanaik and Venayagamoorthy 2019;Masunaga et al. 2022). The average L O /L E was 2.2; therefore, well-mixed turbulent mixing would be expected to occur during our microstructure observations. Figure 6b compares the vertical eddy diffusivity computed from Eq. (1) (K ρ ) and Eq. (9) (K ρLE ). The correlation coefficient (r) between K ρ and K ρLE was 0.83 (p < 0.01), and the standard deviation between them was less than O(1) (Fig. 6c). Thus, the vertical eddy diffusivity can be inferred from L E with a standard deviation error less than O(1), without the need for microstructure shear data.
We compared the Ellison length scale obtained from the mooring and microstructure observations (Fig. 6d). The hourly averaged L E from the SBE-56 at each moored thermistor depth (1.2, 2.4, 3.6, and 4.8 m) was compared with the vertically averaged L E estimated from the VMP-250 within bins of the moored-thermistor depths ± 0.6 m. The correlation coefficient between L E from the mooring observations and that from the microstructure observations was 0.54 (p < 0.05). K ρLE and K ρI18 from the mooring observations were also compared with K ρ from the microstructure surveys (Fig. 6d, e). The correlation coefficients between K ρLE -K ρ and between K ρI18 -K ρ were both 0.44 (p < 0.1). K ρLE and K ρI18 were consistent with each other (Figs. 2e and 7); there was a strong correlation between them (r = 0.93, Fig. 6 a-c Comparisions of the mixing parameters estimated from the VMP-250 microstructure profiler; a Ozmidov length scale (L O ) as a function of the Ellison scale (L E ); b, c comparison of the vertical eddy diffusivity from two methods, K ρ and K ρLE . d-f Comparisons of the mixing parameters determined from the microstructure and mooring observations. The red line in a indicates the regression of L O against L E . Vertical black dashed lines and red curve in c show the standard deviation and fitted first mode normal distribution, respectively p < 0.01), and the standard deviation of log 10 (K ρLE /K ρI18 ) was 0.38 (Fig. 7).
Although K ρLE and K ρI18 show reasonable agreements with K ρ from the direct microstructure measurements with significant consistency (Fig. 6e, f) (r = 0.44, p < 0.1), they do not fit K ρ very well. There are two possible reasons for generating errors between K ρLE /K ρI18 and K ρ . First, the mooring and VMP-250 observations did not measure exactly the same waters. Due to technical difficulties, the horizontal distance between these two surveys was approximately 10-30 m. Intermittent turbulent patches are generally exhibited with horizontal scales of several tens of meters (e.g., Rosenblum and Mamorino 1990;Masunaga and Yamazaki 2014). Therefore, the distance between the two observations may cause errors.
Second, the definition of the constant mixing efficiency, Γ = 0.2 in Eq. (1) and (9), is supposed to be causing errors. It is well known that Γ varies O(10 -3 -1) with mixing conditions, and previous observational and numerical studies have shown that Γ varies O(10 -3 -1) with mixing conditions and Γ increases as the buoyancy Reynolds number decreases (e.g., Shih et al. 2005;Garanaik and Venayagamoorthy 2019;Masunaga et al. 2022). The buoyancy Reynolds number is represented as follows: Probability density functions (PDFs) of Re b for the four VMP-250 period are shown in Fig. 5m-p. The distribution of Re b was generally observed on the lower side and 78% of the total Re b was lower than 10 in the morning of August 5 (Fig. 5m). On other hand, during other periods, Re b distributed on the higher side and 90, 60, and 83% of the total Re b was higher than 100 in the afternoon of August 5, and the morning and afternoon of August 7, respectively. The lowest and highest values of Re b were 3.3 × 10 -2 and 2.0 × 10 6 , respectively. Re b widely ranged within O(10 -2 -10 6 ), which might result in errors between the eddy diffusivities from the three methods with the constant Γ. As described above, although several issues were involved in the eddy diffusivity estimations, the L E -based vertical eddy diffusivity estimation method provided reasonable vertical eddy diffusivity values with errors roughly less than O(1) (Fig. 6c). We assumed that K ρLE and K ρI18 indicate appropriate eddy diffusivity values. This study used the vertical eddy diffusivity, K ρLE estimated using Eq. (9), in the following discussion. K ρLE near the surface (1.2 m depth) ranged from 1.3 × 10 -8 to 8.9 × 10 -3 m 2 s −1 and fluctuated daily (Fig. 3e).

Daily vertical mixing structure
The observed results from the microstructure profiler showed intensified mixing during the afternoon (Fig. 5).
Daily composite plots of the observed parameters were constructed to investigate the daily cycle mixing (Fig. 8).
Here, the daily composite represents the average value at the same time on each day of the study period. The surface temperature was high due to the daytime heat radiation (Fig. 8a, f). The peak surface temperature appeared at approximately 12:00, and heat was gradually transported downward in the afternoon (12:00-18:00) (Fig. 8a). The thermocline also moved deeper in the early afternoon, and the maximum N 2 layer was located approximately depth of 4.5 m during the night (Fig. 8c). The north-westward wind was intensified with an amplitude of 5.5 m s −1 in the afternoon. The wind velocity spectra showed a clear peak of the meridional wind component at 24 h (Fig. 9a). As the wind speed increased between 8:00 and 18:00, the vertically averaged vertical diffusivity increased from 2.5 × 10 -6 to 3.6 × 10 -4 m 2 s −1 (Fig. 8g). The daily composite Ri was also consistent with the observed daily mixing structure (Fig. 8e). I MO decreased in the morning (7:00-11:00) as the stratification increased, and it increased in the afternoon (12:00-16:00) as vertical mixing intensified (Fig. 8i). The dark gray shaded areas in Fig. 8i indicate I MO < 0 (i.e., negative buoyancy flux appeared during these periods). This implies that convective mixing likely occurred at night. One of the possible reasons for the generation of intensified wind in the afternoon is the sea breeze caused by differential heating of the land and ocean. Sea breezes have been studied for over a century (Davis et al. 1890), and they occur on coastlines worldwide (Gille et al. 2003). Zhang et al. (2009) andSobarzo et al. (2010) reported that resonant oceanic motions occur due to sea breezes along continental shelves during summer. Woodson et al. (2007) showed diurnal upwelling caused by sea breezes in a bay. Furthermore, sea breezes are recognized as important contributors to enhanced turbulent mixing and gas exchange in estuaries (Orton et al. 2010). As the Japan mainland faces Fig. 7 Comparisons of the vertical eddy diffusivity estimated by the two methods from the mooring observations, K ρLE and K ρI18 . The red line in a indicates the regression of L O against L E . The thick and thin red lines in a show the mean and standard deviation, respectively. Vertical black dashed lines and red curve in b show the standard deviation and fitted first mode normal distribution, respectively the Pacific Ocean, intensified sea breezes are expected in summer (Yoshikado 1981). The daily averaged meridional wind speed from 12:00 to 18:00 was highly correlated with the maximum air temperature (r = 0.79, p < 0.01, Fig. 9b). Therefore, the intensified wind in the afternoon was attributed to sea breezes caused by differential heating in the vicinity of the study area and the Pacific Ocean.
The daily averaged stability parameter Ri was explained by the I MO (r = − 0.50, p < 0.05) (Fig. 9c). It is well known that a mixing parameter can be represented by another mixing parameter with an appropriate slope (e.g., Garanaik and Venayagamoorthy 2019; Masunaga et al. 2022).
Ri and I MO are represented as (B/S) 1 and (S/B) 3/2 , respectively, where B is the buoyancy term and S is the shear term. The relationship between Ri and I MO can be written as Ri ~ I MO −2/3 . The comparison of Ri and I MO shown in Fig. 9c follows the theoretical relationship of Ri ~ I MO −2/3 (the slope of -2/3 is shown as a red line in Fig. 9c).Thus, the observed mixing and stability can be appropriately explained by atmospheric forcing. To further investigate mixing in response to atmospheric forcing, two other I MO values were computed using long timescale (low-frequency) wind (I MOlow ) and short timescale (highfrequency) wind (I MOhigh ). The cut-off frequency of the low-pass and high-pass filters for these two wind components was set at 9.3 × 10 -6 s −1 (equivalent to a cycle of 30 h). I MOlow was significantly correlated with Ri (r = -0.51, p < 0.05); however, the stability parameter computed from high-frequency winds, I MOhigh , was not significantly correlated with Ri (r = − 0.23, p > 0.1) (Fig. 9d, e). Therefore, the short timescale sea breeze explains intensified mixing less than the daily timescale breeze; however, it does not explain the stability of the water column on a daily or longer timescale. For such timescales, long-timescale wind forcing, such as under synoptic scales, may control the mixing and stability of the water column.

Oxygen transport
The time series of the DO is plotted in Fig. 10a, b. The DO distribution was largely influenced by vertical mixing and low DO concentration (e.g., DO < 3 mg L −1 ) intermittently appeared below the mixed layer depth (Fig. 10a). The DO in the bottom layer decreased at a rate of approximately 3 mg L −1 day −1 (approximately 0.13 mg L −1 h −1 ) under strongly stratified conditions ( Fig. 10b; e.g., for August 12-16, a slope of -3 mg L −1 day −1 is shown as a thick gray dashed line). The decrease in the bottom-layer DO was likely due to biological oxygen consumption. To evaluate vertical DO Fig. 9 a Power spectra of zonal (u w ) and meridional (v w ) wind components, and b daily averaged meridional wind speed obtained from 12:00 to 18:00 as a function of the daily maximum air temperature. c-e The daily and depth-averaged Richardson numbers ( where C DO is the DO concentration. The DO flux was computed from the vertically averaged K ρLE obtained at the four layers, and DO gradient between the surface (~ 0.5 m) and bottom (~ 6.0 m) layers. The time series of the DO flux is shown in Fig. 10c. The positive DO flux is the flux toward the lakebed. Negative DO flux values occupied only 2.4% of the total estimated DO flux and were, therefore, ignored in the plot. The DO flux varied over a wide range between O(10 -9 ) and O(10 -3 ) mg L −1 m s −1 . The DO supply from the surface to the lakebed was estimated by the DO flux and height of the water column (h = 6.5 m). Under high DO flux conditions, F DO = 1.0 × 10 -3 mg L −1 m s −1 , and the hourly DO supply toward the lakebed was computed as 0.55 mg L −1 h −1 . This intensified DO flux was much greater than the decrease in DO under strongly stratified conditions (approximately 0.13 mg L −1 h −1 ). Therefore, the bottomlayer DO would increase under conditions of intensified DO flux. In contrast, under conditions of limited DO flux (F DO = 1.0 × 10 -8 mg L −1 m s −1 ), the hourly DO supply was 5.3 × 10 -6 mg L −1 h −1 , and the vertical DO supply was not expected.
The results from the WOD and SOD laboratory experiments are presented in Fig. 11. The change in DO (DDO) for the SOD experiments is represented as the ΔDO in milligrams per unit area (mg m −2 ; Fig. 11b). The ΔDO rate was computed using the linear least square method and was −2.7 mg L −1 day −1 (−0.11 mg L −1 h −1 ) and (16) −995.2 mg m −2 day −1 (−41.5 mg m −2 h −1 ) for the WOD and SOD experiments, respectively. The vertically averaged ΔDO rate (in mg L −1 t −1 , where t is the time unit) due to the SOD was estimated by dividing ΔDO by h × 10 3 (h is the height of the water column, 6.5 m) and was computed as 0.15 mg L −1 day −1 (−6.4 × 10 -3 mg L −1 h −1 ). According to the results of the WOD and SOD experiments, the oxygen demand in the water column was determined to be −2.8 mg L −1 day −1 (− 0.12 mg L −1 h −1 ). This estimated oxygen demand is consistent with the rate of change in DO in the bottom layer during the strongly stratified period from the field observations (approximately 3 mg L −1 day −1 ,  Fig. 10b). The extinction coefficient of visible light in the lake is high (about 1-2 m −1 ) due to eutrophication (Nakamura and Aizaki 2016). Under these conditions, over 99% of the incoming light from the lake surface is absorbed by the lake water before reaching the bottom layer (6.5 m depth). Therefore, phytoplankton photosynthesis near the lakebed might be prohibited because of light limitation, so the increase in DO near the bottom must have been due to the vertical DO supply from the upper layer. Daily composite plots of the DO and the rates of change in DO and DO flux are shown in Fig. 12. DO increased from sunrise to 18:00 (approximately one hour before sunset) (Fig. 12a, b), which was likely due to phytoplankton photosynthesis in the surface layer. The daily peak of surface DO occurred at 15:00, after which oxygen was gradually transported toward the bottom layer where photosynthesis was prohibited (Fig. 12a). Near-surface DO (0-2 m depth) increased from 9:00 to 12:00, and DO in the deeper layer (2-6 m) increased in the afternoon (12:00-18:00) (Fig. 12c). The vertical DO transport in the afternoon was consistent with intensified mixing due to the sea breeze shown in Fig. 8. The vertically averaged DO increased during the day and decreased during the night (Fig. 12b, d), which can be simply explained by phytoplankton photosynthesis, and biological and chemical oxygen consumption, respectively. Here, daytime and nighttime were considered as the durations between the sunrise and sunset, and vice versa. The vertically averaged rate of change in DO (ΔDO/Δt) was 0.13 and -0.15 mg L −1 h −1 in daytime and nighttime, respectively.
The DO flux increased from sunrise to sunset, and the maximum DO flux, 2.9 × 10 -4 mg L −1 m s −1 , was observed at 18:00 (Fig. 12e). This was the same time at which the highest eddy diffusivity was observed (Fig. 8g). The minimum DO flux, 1.7 × 10 -6 mg L −1 m s −1 , was observed in the early morning at 5:00. The maximum DO flux induced by sea breeze was O(100) times higher than the DO flux during the night. It is expected that convective mixing should occur due to nighttime surface cooling, as reported by Monismith et al. (1990). Indeed, our results showed that a negative heat flux (heat transport from lake water toward the atmosphere) was observed during nighttime (Fig. 8f, i). However, DO in the bottom layer did not increase during the night, except around 6:00 (Fig. 12c), and stratification appeared throughout the entire night on average (Fig. 8ac). We conclude that the vertical DO supply is largely determined by intensified mixing due to the afternoon see breeze and that nighttime convective mixing partly influence the vertical DO supply.

Conclusions
This study conducted mooring and microstructure observations in a shallow stratified lake. We found intensified mixing due to sea breezes leading to the vertical oxygen supply from the surface toward the bottom. The vertical eddy diffusivity was estimated using the following three independent methods: (1) direct microstructure shear measurements using a free-falling or rising microstructure profiler; (2) mooring measurements using the Ellison scale (L E ) and vertical background shear (S), as proposed by Ivey et al. (2018); and (3) mooring measurements of L E and buoyancy frequency (N) with relationships between L E and the Ozmidov scale (L O ). The eddy diffusivity computed from the three methods was compatible with each other. It is worth noting that the third method allows us to estimate vertical eddy diffusivity only using a thermistor chain if the relation between L E and L O holds. The observed eddy diffusivity showed intensified vertical mixing and DO flux in the afternoon due to increased sea breezes. The sea breeze was caused by differential heating between the land around the lake and the Pacific Ocean. The results indicate that strong daytime heat flux worked upon lake water mixing in two ways: (1) surface heating intensified near the surface layer stratification, suppressing vertical mixing; and (2) differential heating in the vicinity of the lake and Pacific Ocean generated sea breezes in the afternoon, resulting in enhanced vertical mixing. It should be noted that the effects of strong daytime heating both suppressed and enhanced vertical mixing.
Understanding the oxygen supply and associated transport processes is an important issue for ecosystem management in lakes and reservoirs. This study succeeded in observing the vertical eddy diffusivity and oxygen flux using a highfrequency sampling mooring thermistor and conventional DO surveys. We found short-term daily mixing and a large vertical oxygen flux due to sea breezes with a duration of several hours, which has not been found in previous studies. Our methodology can be used not only for our study site, but it can also be used for other lakes and reservoirs. Conducting such observations of mixing and oxygen flux using a mooring system with high-frequency thermistors would contribute to improving the understanding of mixing and oxygen supply in lakes. operating the observation systems at the Kamaya Monitoring Station (Japan Water Agency).We also thank the editor and two anonymous reviewers for their thoughtful suggestions to improve the manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.