Reconsideration of winds, wind waves, and turbulence in simulating wind-driven currents of shallow lakes in the Wave-current Coupled Model (WCCM) version 1.0

Winds, wind waves, and turbulence are essential variables and playing critical role in regulating a series of physical and biogeochemical processes in large shallow lakes. However, parameterizing winds, waves, currents and turbulence and simulating the interaction between them in large shallow lakes haven’t been evaluated strictly because of a lack of field observations of lake hydrodynamics process. To address this problem, two process-based field observations were conducted to record the development of summer and winter wind-driven currents in Lake Taihu, a large shallow lake in 15 China. Based on these observations and numerical experiments, a wave-current coupled model (WCCM) is developed by rebuilding expression of wind drag coefficient, introducing wave-induced radiation stress, and adopting a simple turbulence scheme, and then used to simulate wind-driven currents in Lake Taihu. The results show that, the WCCM can accurately simulate the upwelling process resulting from the wind-driven currents during the field observations. Comparing with other model, there is a 42.9% increase of WCCM-simulated current speed which is mainly attributed to the new expression of 20 wind drag coefficient. Meanwhile WCCM-simulated current direction and field are also improved due to the introduction of wave-induced radiation stress. Furthermore, the use of the simple turbulent scheme in the WCCM makes the simulation of the upwelling processes more efficient. The WCCM provides a sound basis for simulating shallow lake ecosystems.


Three-dimensional lake current model
The LCM is developed to simulate the water level and lake currents.

Governing equations 95
The governing equations for lake currents in the Cartesian coordinate system (Fig. 2) consist of the continuity equation and momentum equations (Koue et al., 2018). To eliminate the influence of lakebed topography on the lake current simulations, the sigma (σ) coordinate system is introduced in the vertical direction (Fig. 2).
Based on the rule of derivation of a composite function, the continuity equation and momentum equation in the Cartesian coordinate system (xʹ, yʹ, z, tʹ) are transformed into the σ coordinate system (x, y, σ, t) using Eqs. A1-1 to A1-5.
Where: u, v, and w are the components of the current velocity in the x-, y-, and σ-directions (m s −1 , m s −1 , s −1 ), respectively; h, ζ, and H are the lakebed elevation, water level, and water depth (m), respectively; f is the Coriolis force (s −1 ) defined by f = 105 2ωsinφ, where ω is the rotational angular velocity of the earth and φ is the geographic latitude; F x and F y are the waveinduced radiation stress in the x-and y-directions, respectively; ρ and ρ 0 are the water and reference density (kg m −3 ), respectively; g is the gravitational acceleration; A H and A V are the horizontal and vertical eddy viscosity (m 2 s −1 ), respectively; and ε U , and ε v , are the secondary terms introduced by the coordinate system transformation (Eqs. A2-1 to A2-6).

Turbulence scheme 110
To improve the calculation efficiency, the value of the vertical eddy viscosity (A V ) is estimated using the Prandtl length l and the Richardson number (R i ).
l and R i are given by: where κ is the von Ká rmá n constant, z 0 is the roughness height of the lakebed, and r s is the roughness height of the lake surface.

Boundary conditions
Wind stress at the lake surface: 120 V ( , ) = a s √ w 2 + w 2 ( w , w ) , where ρ a is the air density, u w and v w are the wind speed components in the x-and y-directions at 10 m above the lake surface (m s −1 ), respectively, and C s is the wind drag coefficient.
Here, according to the process-based observations and model calibration described in the following section 4.1, we define a new expression of C s that considers the discontinuity of changing trend and directionality of wind momentum transmission, 125 which differs from previously reported expressions of C s .
Friction at the lakebed: where C B is the bottom friction coefficient that is given by:

Wave-induced radiation stress
Wave-current interaction is a complicated process (Mellor, 2008). Up to now, it is not fully understood. Longuet-Higgins and Stewart (1964) firstly proposed the concept of wave-induced radiation stress and Sun et al. (2006) derived the 135 expressions of the stress for three-dimensional current numerical models: where H S is the significant wave height (m), T 0 is the wave period (s), L is the wavelength (m), θ m is the mean wave direction, 140 and θ 1 is the angle between the mean wave direction and geographical east direction.

Solution of equations
The splitting mode technique (Blumberg and Mellor, 1987) and alternation direction implicit difference scheme (Butler, 1980) are used to discretize the Equations (1 to 3) on the staggered grid ( Figs. 2 and 3). The detailed description of the solution of equations is indicated in Appendix A3. 145

Simulating WAves Nearshore model
In view of the importance of wind waves in the hydrodynamic and ecological processes of shallow lakes, the SWAN model, which has been proven suitable for simulating the wind waves in Lake Taihu (Wang et al., 2016;Wu et al., 2019;Xu et al., 2013) was used to simulate the spatiotemporal variation of wind waves in the lake. The governing equation for SWAN is the wave action balance equation: 150 where N is the action density spectrum, t, x, and y are the time and horizontal coordinate directions, respectively, σ 1 is the relative frequency, θ is the wave direction, c x , c y , 1 , and denote the wave propagation velocity in x, y, σ 1 , and θ space, respectively, and S is the source in terms of energy density representing the effects of generation, dissipation, and nonlinear wave-wave interactions. H S , T 0 , L and θ m are deduced from the value of N(x, y, t, σ 1 , θ) (Booij et al., 2004). 155 The action balance equation is solved in the Cartesian coordinate system using the first-order upwind scheme of the finite difference method (Booij and Holthuijsen, 1999;Booij et al., 2004).

Two-way coupling of the LCM with SWAN
SWAN and LCM were coupled together to establish the WCCM model (Fig. 3). The current speeds u and v, and the water level ζ that are computed by the LCM model are inputs of the SWAN model. On the other hand, H S , T 0 , L and θ m that are 160 computed by the SWAN model are used as inputs in the LCM model, for the computation of the wave induced radiation stresses F x and F y (Eqs. (13) and (14)).

Configuration of the WCCM in Lake Taihu
The WCCM is used to simulate the wind waves and lake currents of Lake Taihu during the period of process-based field observations. The computational domain of Lake Taihu (Fig.1) for the LCM is divided in 72 × 72 = 5184 cells at 1 km 165 resolution. The water column is divided into five layers in the vertical direction and the time step is 30 s. The value of α is 0.5. Lake Taihu is considered as a closed lake for the simulation because the influence of inflows and outflows on the current field is very small compared to the influence of the wind stress (Li et al., 2011;Wu et al., 2018;Zhao et al., 2012). The simulations therefore disregard the inflows and outflows. The model inputs at the air-water boundary include air temperature, surface air pressure, cloud cover, relative humidity, wind speed and direction provided by the LCWS and TLLER stations (Fig. 1). Among them, the measured wind speed at 5 m above the water surface was adjusted to 10 m (Wu et al., 2018) using the method suggested by the Coastal Engineering Research Center (1984). The initial condition for the water level was obtained by interpolation of the values of the water levels measured at stations WL1-WL5 at the beginning of the model integration. The initial water temperature was set to the measured values recorded by the ADP at the beginning of the model 175 integration and the current speed was initialized by 0 m s −1 .
Ten parameters need to be determined for the simulation of the LCM (Table 1). Among them, φ, g, κ, and ρ a are constants, while A v and C B values can be calculated from the values of κ, z 0 and r s . A H and z 0 values are the same values as the ones used for the EFDC (Environmental Fluid Dynamics Code) and r s is set to 0.01 (Table 1). Being described in the following section, the EFDC, which is a hydrodynamic numerical model, is used here to evaluate the WCCM's performance. The expression of 180 the wind drag coefficient is designed and calibrated using the process-based observation data of 2015.
The mesh of the SWAN model is the same as the horizontal mesh of the LCM. Considering the randomness of wind waves, the characteristic values of wind waves are typically represented by the statistical values of the high frequency pressure records over a 10-min period. The time increment of the SWAN model was therefore set to 600 s. The frequency band was set to 0.04-4 Hz and the wave direction ranged from 0° to 360° with an increment of 6°. The second generation mode was 185 used to calculate the source term (e.g., wind input, depth-induced wave breaking, bottom friction, triads). The parameter cdrag of the SWAN model was set to 0.00133 and the Collins bottom friction coefficient was set to 0.025. The calibration and validation of these parameters have been reported in previous studies (Xu et al., 2013;Wang et al., 2016). Because the influence of lake currents on the SWAN-simulated wind waves has already been analyzed in Lake Taihu (Li et al., 2007), only the WCCM simulation of lake currents was evaluated in this study.

Statistical analysis
To evaluate the WCCM model performance, the mean absolute error (MAE), the root mean square error (RMSE), and the correlation coefficient (r) between the measured and simulated values at both significance levels of p < 0.05 and p < 0.01 are considered (Koue et al., 2018). The magnitude of lake current speed is expressed as the mean ± standard deviation.
The mean absolute error of the horizontal current direction (MAE UVD ; Carvalho et al., 2012) is used to compare the simulated 200 and measured values: In addition, ArcGIS 10.2 (ESRI Inc., USA) was used to process the spatial data and Tecplot 360 (Tecplot Inc., USA) was used to draw contours of the water level, current field, and streamtraces. 205

Comparison between the WCCM and EFDC
Comparison between different models is a useful method to study currents in large water bodies (Huang et al., 2010;Morey et al., 2020). The EFDC is one of the most widely used models for shallow lakes worldwide (Chen et al., 2020) and offers a general-purpose modeling package to simulate three-dimensional flow, transport, and biogeochemical processes in surface water systems (Ji et al., 2001;Ji, 2008). The EFDC has been successfully applied in Lake Taihu modelling (Li et al., 2011;210 Li et al., 2015;Wang et al., 2013). The EFDC hydrodynamic model was developed by Hamrick (1992) and its governing equations are the same as Eqs. (1)-(3).
It uses the splitting mode technique to solve the continuity equation and momentum equation in the σ coordinate system. The Mellor-Yamada turbulence model is used in EFDC to calculate the vertical eddy viscosity (Ji et al., 2001). The wind stresses in the EFDC is calculated using the following equations (Hamrick, 1992;Li et al., 2015): The mesh used for the simulation with the EFDC is the same as the LCM and WCCM. After consulting the authors of the uncertainty and sensitivity analysis performed on the hydrodynamic parameters of the EFDC for Lake Taihu (Li et al., 2015), the optimal horizontal eddy viscosity was set to 1 m 2 s −1 , the roughness height to 0.005 m, and w s to 0.7. 220

Numerical Experiments
Four numerical experiments were designed to evaluate the accuracy of the WCCM and to identify the relative importance of winds, wind waves, and turbulence in improving simulation of the wind-driven currents as follows:  Experiment 1, denoted EFDC: numerical simulation of the lake currents using the EFDC. In this experiment, the Mellor-Yamada turbulence scheme is used and the drag coefficient is given by Eqs (18)-(19), but the 225 wave-induced radiation stress is not considered (no coupling with SWAN);  Experiment 2, denoted LCM_1: numerical simulation of the lake currents using the LCM, with the same expression of the drag coefficient as in EFDC (Eqs (18)-(19)), but a different turbulence scheme that is given by Eqs.(4)-(6), and still without consideration of wave-induced radiation stress;  Experiment 3, denoted LCM_2: same experiment as LCM_1 with a different expression of the drag 230 coefficient that is given by Eqs. (8)  Experiment 4, denoted WCCM: same experiment as LCM_2 but with consideration of wave-induced radiation stress to achieve the two-way coupling model.

Summer observation and model calibration in 2015 235
The average wind speed over Lake Taihu between 0:00 on August 8 and 0:00 on August 12, 2015 was 8.6 m s −1 (Fig. 4) with a maximum of 13.5 m s −1 at 13:00 on August 10, corresponding to a wind direction of 107.5°. Lake Taihu experienced a strong southeast-east wind event during the observation.  (Fig. 7). The surface current field simulated by these two models mainly flows from southeast to northwest, which can be further demonstrated by the simultaneous stream traces (Fig. B.1). The middle and bottom current fields of the southern part of the lake are consistent 250 with the surface current field, but those in the center and northern parts of the lake mainly flow from southwest to northeast.
A major difference between the WCCM-and EFDC-simulated current fields is in that the current speed simulated by the former is significantly higher (Fig. 7). There are vortexes produced by the WCCM in the upwind area, such as in Xukou Bay and northwest of Xishan Island (Fig. B.1). In contrast, the vortexes simulated by the EFDC tend to be located in the downwind area, such as Zhushan Bay and Meiliang Bay (Fig. B.1). 255
The EFDC, LCM_1, LCM_2, WCCM-simulated water levels at each water level station significantly correlate with the measured values (p < 0.01; Table 4 The water level contours simulated by the WCCM at 22:00 on December 26, 2018, corresponding to the time of the maximum wind speed, are similar with those by the EFDC. They show a decease trend from southwest to northeast (Fig. 11). 270 The surface current fields simulated by these two models mainly flow from north to south, which can be further demonstrated by the simultaneous stream traces (Fig. B.2). The middle and bottom current fields mainly flow from northwest to southeast.
The main difference between the WCCM-and EFDC-simulated current fields is that the current speed simulated by the former is significantly higher (Fig. 11). Clockwise vortexes formed in Gonghu Bay in the surface, middle, and bottom 275 current fields simulated by the EFDC (Fig. B.2), whereas this clockwise vortex is only located in the middle current field simulated by the WCCM.

Discussion
Influenced by the strong southeast-east wind event during the summer observation in 2015, the maximum water level difference at 12:00 on August 10 between WL1 located in the downwind lake area and WL4 station located in the upwind 280 lake area was 0.66 m (Fig. 5). Before this maximum, all of the measured surface, middle, and bottom currents flow along wind direction and their speed significantly increased (Fig. 6). It can be concluded that the strong southeast-east winds drive whole water column at the LCWS station to form monolayer wind-driven currents, and then result in a downwind upwelling (Wu et al., 2018). Similarly, generated by the strong north-northeast wind event during the winter observation in 2018, winddriven currents also result in a downwind upwelling, despite both of the wind-driven currents and upwelling of the winter 285 observation are weaker than those of the summer observation (Fig. 10). These summer and winter upwelling processes provided us excellent chance to evaluate the performance of the WCCM in Lake Taihu. whereas the water level is simulated at similar accuracy as by the EFDC. The WCCM can accurately simulate the winddriven currents, and subsequent the downwind upwelling in Lake Taihu.

Wind drag coefficient 295
This variable is a key parameter for hydrodynamic numerical models. The EFDC parameter sensitivity analysis shows that the wind drag coefficient is the most sensitive parameter for simulating current velocity in Lake Taihu (Li et al., 2015). The correlation coefficients between the simulated and measured current speeds of LCM_2 and WCCM, which use the new expression of C s , considering the discontinuity of changing trend and directionality of wind momentum transmission, are significantly greater than those of EFDC and LCM_1 (Tables 3 and 5). This implies that the use of the new C s mainly 300 contributes to the increase of the correlation coefficient.
A piecewise function is firstly proposed in this study to describe the discontinuity of the changing trend of the wind momentum transmission. The changing trend of the wind drag coefficient at the water surface is discontinuous (Wu, 1980). The atmospheric surface layer appears to be aerodynamically rough related to wind waves for wind speeds > 7.5 m s −1 and aerodynamically smooth for wind speeds < 3 m s −1 (Wu, 1980). The transmission efficiency of wind momentum to the water 305 under aerodynamically rough conditions is higher than that under aerodynamically smooth conditions (Lükő et al., 2020).
Field observations of Lake Taihu indicate that C s increases with wind speed (Xiao et al., 2013). A critical wind speed of 7.5 m s −1 is therefore adopted to describe the discontinuity. As shown in Eqs. (8) and (9), a logistic curve is used to describe the increase of C s for wind speeds > 7.5 m s −1 , otherwise C s is a constant value.
The directionality of wind momentum transmission is further addressed using different C s values in the x-and y-directions. 310 There have been numerous expressions designed to calculate wind drag coefficient based on ocean environments (Geernaert et al., 1987;Large and Pond, 1981;Lükő et al., 2020;Wu, 1980;Zhou et al., 2009). However, few expressions consider the directionality of wind momentum transmission. There is a contradiction that: C s increases with wind speed (Lükő et al., 2020;Xiao et al., 2013), while the tilt of water surface along wind direction in large shallow lakes with limited water depth and fetch due to the upwelling will decrease the wind momentum transmission efficiency. Because of this contradiction, the 315 same transmission efficiency of wind momentum is used in both of being perpendicular and parallel to wind directions, which may over-or under-estimate the wind drag coefficient in any one direction.

Wave-induced radiation stress
This is the first time for wave-induced radiation stress to be considered in simulating wind-driven currents in large shallow lakes. The results show that it can improve the simulated current direction. In 2015, the MAE UVD values of the LCM_2 320 (average MAE UVD of 56.3°; Table 3) are greater than those of the WCCM (average MAE UVD of 52.9°; Table 4). A similar result can be achieved by comparing the MAE UVD values between the LCM_2 and WCCM in 2018 (Table 5). Moreover, the correlation coefficients of LCM_2 in 2018 are slightly less than those of the WCCM in 2018 (Table 5), which implies that wave-induced radiation stress can also contribute to the improvement of the WCCM-simulated current speed. The comparison between the WCCM-and EFDC-simulated current fields further demonstrates the importance of wave-325 induced radiation stress. Although the current field simulated by the WCCM is similar to that by the EFDC, the vortex locations simulated by these models are quite different. In 2015, the middle and bottom current fields simulated by the EFDC exhibit counterclockwise vortexes in Zhushan Bay and Meiliang Bay (Fig. B.1), which are located in the downwind area, but the current fields simulated by the WCCM do not show the same phenomenon. This is because the interaction between wind waves and lake currents in the downwind area is violent owing to wave deformation resulting from the 330 shallow water and lakeshore. The wave-induced radiation stress makes the vortex less likely to form in this area. Conversely, the middle and bottom current fields simulated by the LCM_2 without wave-induced radiation stress also show counterclockwise vortexes in Zhushan Bay and Meiliang Bay (Fig. B.3), which is similar to the result of the EFDC.

Vertical eddy viscosity
Comparing with other variables, the vertical eddy viscosity play a less prominent part in the development of the wind-driven 335 currents. In our study, Mellor-Yamada level-2.5 turbulence closure model (Mellor and Yamada, 1982;Ji et al., 2001) is adopted in the EFDC and the other parameters are determined after parametric uncertainty and sensitivity analysis (Li et al., 2015), while a simple turbulence scheme (equation (4)-(6)) is adopted in the LCM_1. However, the accuracy of the LCM_1 is rather similar to that of the EFDC (Tables 2, 3, 4, and 5), which implies that the high-order turbulence scheme does not improve the lake current simulations (Koue et al. 2018), while the simple turbulence scheme makes the WCCM more 340 efficient.

Conclusion
The strong summer or winter winds generate wind-driven currents in Lake Taihu, and subsequent results in downwind upwelling. Based these processes and numerical experiments, a wave-current coupled model (WCCM) is developed by reconsidering the expressions of winds, wind waves, and turbulence. It can accurately simulate the development of wind-345 driven currents with a 42.9% increase of simulated current speed compared with the EFDC results of 2018. The new expression of the wind drag coefficient is mainly responsible to increase of the correlation coefficient between the WCCMsimulated and measured current speeds. The introduction of wave-induced radiation stress can contribute to the improvement of the simulated current direction and fields, and slightly improve the simulation of current speed. Moreover, the simple parameterized turbulence scheme is sufficient for the simulation of wind-driven currents in Lake Taihu. 350 It should be noted that despite the performance of the numerical models had been greatly improved, the correlation between simulated and measured current speed remains low, especially for the 2018 observation. Actually, few model studies are reported to conduct this correlation analysis because of lack data or worse results. Therefore, we urge that more processbased field observations are required to help us to fully understand the real hydrodynamic characteristics of large shallow lakes and further improve the performance of shallow lake current models.