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

. Wind stress, wind waves, and turbulence are essential variables and play a critical role in regulating a series of physical and biogeochemical processes in large shallow lakes. However, the parameterization of these variables and simulation of their interactions in large shallow lakes have not been strictly evaluated owing to a lack of field observations of lake hydrodynamic processes. 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 China. Then, a Wave and 15 Current Coupled Model (WCCM) is developed based on these observations and numerical experiments by rebuilding the wind drag coefficient expression, introducing wave-induced radiation stress, and adopting a simple turbulence scheme to simulate wind-driven currents in Lake Taihu. The results show that the WCCM can accurately simulate the upwelling process driven by wind-driven currents during the field observations. A comparison with a reference model indicates a 42.9% increase of the WCCM-simulated current speed, which is mainly attributed to the new wind drag coefficient expression. The 20 WCCM-simulated current direction and field are also improved owing to the introduction of wave-induced radiation stress. The use of the simple turbulent scheme in the WCCM improves the efficiency of the upwelling process simulation. The WCCM thus provides a sound basis for simulating shallow lake ecosystems.

Many efforts have been made on coupled current-wave model development, especially on the coupling of the Simulating WAves Nearshore model (SWAN; Booij and Holthuijsen, 1999) with existing three-dimensional current models (Chen et al., 2018;Liu et al., 2011;Warner et al., 2008;Wu et al., 2011). However, due to the difficulty in modifying existing model 100 codes (Chen et al., 2018), most coupling models have been developed using third-party software (e.g. Model Coupling Toolkit) rather than by directly merging the original codes. However, this is not yet an efficient way to modify the descriptions of some key variables in these models. Herein, a Wave and Current Coupled Model (WCCM) is developed by merging the codes of a three-dimensional lake current model (LCM) and SWAN.

Three-dimensional lake current model 105
Although most current models largely use same governing equations and solution methods, differences in the programming languages, operating environment, mesh, and description of key processes or parameters impede a full understanding of these models to allow further code modification. It is thus preferable to develop a new model for determining a suitable description of wind stress, wind waves, and turbulence. The LCM model with a concise and efficient programming is therefore developed to simulate water temperature, water level, and lake currents based on the classic method (Blumberg and 110 Mellor, 1987),

Governing equations
The governing equations of the LCM in the Cartesian coordinate system (Fig. 2) consist of the continuity equation, momentum equations, temperature equation, and density equation (Koue et al., 2018). The sigma (σ) coordinate system is introduced in the vertical direction to eliminate the influence of lakebed topography on the lake current simulations (Fig. 2). 115 Based on the derivation rule of a composite function, these equations in the Cartesian coordinate system (xʹ, yʹ, z, tʹ) are transformed into the σ coordinate system (x, y, σ, t) using Eqs. A1-1 through 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 = 2ωsinφ, where ω is the rotational angular velocity of the earth and φ is the geographic latitude; F x and F y are the wave-125 induced 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; T is the water temperature (°C ), K H and K V are horizontal and vertical turbulent diffusivity (m 2 s −1 ), respectively, S h and C P are the heat source term and heat capacity, (J m 3 s −1 , 4179.98 J kg −1 °C −1 ), respectively; and ε u , ε v , and ε T are the secondary terms introduced by the coordinate system transformation (Eqs. A2-1 through A2-3). 130 The key parameters and solutions of the continuity equation and momentum equations are demonstrated below, whereas the development and validation of the temperature and density simulations of the LCM will be reported in a separate paper.

Turbulence scheme
To improve the calculation efficiency, the value of the vertical eddy viscosity (A V ) is estimated using the Prandtl length l and Richardson number (R i ). 135 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 140 surface.

Boundary conditions
Wind stress at the lake surface: where ρ a is the air density, u w and v w are the wind speed components in the x-and y-directions 10 m above the lake surface 145 (m s −1 ), respectively, and C s is the wind drag coefficient.
The expression of C s for light winds differs from that for high winds, and a piecewise function is recommended to fit the changes of C s with wind speed (Large and Pond, 1981). A constant (C c ) is often used to represent C s below the critical wind speed (W cr ), while a proportional function is adopted for the increase of C s with wind speed over W cr . However, referring to Geernaert et al. (1987), C s approaches a constant of ~0.003 for wind speeds higher than 20 m s −1 . We therefore propose that 150 a logistic function is more reasonable to derive the expression of C s under high-wind conditions. The wind components in the x-and y-directions are used to calculate C s in the x-and y-directions, respectively.
where (| w |) and (| w |) are the logistic functions. 155 Friction at the lakebed: where C B is the bottom friction coefficient given by:

Wave-induced radiation stress 160
Wave-current interaction is a complicated process (Mellor, 2008) and remains poorly understood. Longuet-Higgins and Stewart (1964) first proposed the concept of wave-induced radiation stress, and Sun et al. (2006) derived the 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, 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, 170 1980) are used to discretize Eqs. (1)-(3) on the staggered grid (Figs. 2, 3). A detailed description of the solution of equations is provided in Appendix A3.

Simulating WAves Nearshore model
In view of the importance of wind waves in the hydrodynamic and ecological processes of shallow lakes, the SWAN model has been frequently used to simulate wind waves in Lake Taihu (Wang et al., 2016;Wu et al., 2019;Xu et al., 2013). The 175 governing equation for the SWAN is the wave action balance equation: 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, which represents the effects of generation, dissipation, and 180 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).
The action balance equation is solved in the Cartesian coordinate system using a 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
The SWAN and LCM were coupled to establish the WCCM model (Fig. 3). The current speeds u and v and water level ζ 185 computed by the LCM model are inputs for the SWAN model. The H S , T 0 , L, and θ m values computed using the SWAN model are used as inputs in the LCM model to compute the wave-induced radiation stresses F x and F y (Eqs. (13), (14)).

Configuration of the WCCM in Lake Taihu
The WCCM is used to simulate wind waves and lake currents in Lake Taihu during the process-based field observation periods. Referring to existing model studies in Lake Taihu (Hu et al., 2006;Mao et al., 2008;Liu et al., 2018), the horizontal 190 computation domain of Lake Taihu (Fig.1) for the LCM is divided into 72 × 72 = 5184 cells with a 1-km resolution to improve the computing efficiency. The water column is divided into five layers in the vertical direction, the time step is 30 s, and the α value 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 with the influence of wind stress (Li et al., 2011;Wu et al., 2018;Zhao et al., 2012). The 195 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, and wind speed and direction collected from the LHWS and TLLER (Fig.   1). The initial condition for the water level was determined via interpolating the water levels values 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 and YSI Sonde at the beginning of the model integration and the current speed was initialized to 0 m s −1 . 200 Ten parameters must be determined for the LCM simulation (Table 1). Among them, φ, g, κ, and ρ a are constants, while the A v and C B values can be calculated from the κ, z 0 and r s values. A H and z 0 values are the same as those used for the EFDC, and r s is set to 0.01 (Table 1).
The parameters in Eqs. (8) and (9) are determined as follows. A critical wind speed of 7.5 m s −1 is used to distinguish between light and high winds by equaling the wind speed for defining aerodynamically rough water surface (Wu, 1980). The 205 expression of the logistic function in Eq.
The SWAN model mesh is the same as the LCM horizontal mesh. Considering their randomness, the characteristic wind wave values 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 used to calculate the source term 215 (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).Our study also verifies that the SWAN in the WCCM can accurately simulated the change of significant wave height at LHWS during the 2018 field observations (Fig. B.1).
Considering the time of the wind peaks and cold start of the WCCM, the hydrodynamics time series of the latter half of the 220 2015 summer observation (from 0:00 on August 8 to 0:00 on August 12, 2015) were used to calibrate the WCCM, and those of the latter half of the 2018 winter observation (from 0:00 on December 26 to 0:00 on December 31, 2018) were used to evaluate the WCCM performance.
The WCCM can be used to simulate the changes of water temperature in Lake Taihu (Fig. B.2), which will be discussed in detail in a separate paper. Here, only the WCCM simulations of the lake currents are evaluated. 225

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

Comparison between the WCCM and EFDC
A comparison between different models is a useful method to study currents in large water bodies (Huang et al., 2010;Morey et al., 2020;Soulignac et al., 2017). 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 240 Lake Taihu (Li et al., 2011;Li et al., 2015;Wang et al., 2013). Here, the EFDC is used to evaluate the WCCM performance. 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 the EFDC to calculate the vertical eddy viscosity (Ji et al., 2001). The wind stress in the EFDC is calculated using the following equations (Hamrick, 1992;Li et al., 2015;Wu et al., 1980): 245 where s ′ and w s are wind drag coefficient and wind shelter coefficient in the EFDC, respectively (Fig. 4).
The mesh used for the EFDC simulation is the same as that in the LCM and WCCM. After consulting with the authors of the uncertainty and sensitivity analysis performed on the hydrodynamic parameters of the EFDC for Lake Taihu (Li et al., 2015), 250 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.

Numerical experiments
Four numerical experiments were designed to evaluate the accuracy of the WCCM and identify the relative importance of wind stress, wind waves, and turbulence in improving the simulation of the wind-driven currents.
 Experiment 1 (EFDC): numerical simulation of the lake currents using the EFDC. The Mellor-Yamada 255 turbulence scheme is used and the drag coefficient is given by Eqs. (20) and (21), but the wave-induced radiation stress is not considered (no coupling with SWAN).
 Experiment 2 (LCM_1): numerical simulation of the lake currents using the LCM with the same drag coefficient expression as in EFDC (Eqs. (20) and (21)), but a different turbulence scheme, as given in Eqs.
 Experiment 4 (WCCM): same experiment as LCM_2 but with considering wave-induced radiation stress to achieve the two-way coupling model.

Summer observation and model calibration in 2015
The average wind speed over Lake Taihu between 0:00 on August 8 and 0:00 on August 12, 2015 was 9.9 m s −1 (Fig. 5), with a maximum of 15.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 2015 summer observation.
The mean water level observed at the five stations was 3.64 ± 0.01 m with a maximum of 4.04 m recorded at the WL1 270 station at 12:00 on August 10 (Fig. 6), corresponding to 3.38 m recorded at the WL4 station. The mean measured surface, middle, and bottom current speeds at the LHWS (Fig. 7) were 5.0 ± 3.0, 5.5 ± 3.5, and 5.4 ± 3.6 cm s −1 , respectively. The contours of the water level simulated by the WCCM at 13:00 on August 10, corresponding to the time of the maximum wind speed, are similar to those of the EFDC simulation and show a decreasing trend from northwest to southeast (Fig. 8).
The surface current field simulated by these two models mainly flows from southeast to northwest, which is further demonstrated by the simultaneous streamtraces (Fig. B.3). The middle and bottom current fields of the southern part of the 280 lake are consistent 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 the significantly higher current speed simulated by the former (Fig. 8). There are vortexes produced by the WCCM in the upwind area, such as in Xukou Bay and northwest of Xishan Island (Fig. B.3). In contrast, the vortexes simulated by the EFDC tend to be located in the downwind 285 area, such as Zhushan Bay and Meiliang Bay.

Winter observation and model validation in 2018
The average wind speed over Lake Taihu is 9.2 m s −1 between 00:00 on December 26 and 00:00 on December 31, 2018 ( Fig.   9) with a maximum of 13.6 m s −1 at 22:00 on December 26, corresponding to a wind direction of 26.3°. Lake Taihu experienced a strong north-northeast wind event during the 2018 winter observation. 290 The mean water level over the five stations was 3.46 ± 0.01 m with a minimum of 3.23 m recorded at the WL5 station at 22:00 on December 26, corresponding to a secondary peak of 3.62 m recorded at the WL3 station (Fig. 10). The mean measured surface, middle, and bottom current speeds at the LHWS (Fig. 11) were 3.7 ± 2.0, 3.5 ± 2.0, and 4.2 ± 2.2 cm s −1 , respectively.
The EFDC, LCM_1, LCM_2, and WCCM-simulated water levels at each water level station significantly correlate with the 295 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 300 maximum wind speed, are similar with those of the EFDC and show a deceasing trend from southwest to northeast (Fig. 12).
The surface current fields simulated by these two models mainly flow from north to south, which can be further demonstrated by the simultaneous streamtraces (Fig. B.4). 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 305 former is significantly higher (Fig. 12). Clockwise vortexes form in Gonghu Bay in the surface, middle, and bottom current fields simulated by the EFDC (Fig. B.4), 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 2015 summer observation, a maximum water level difference 310 of 0.66 m occurred at 12:00 on August 10 between WL1 in the downwind lake area and WL4 station in the upwind lake area (Fig. 6). Prior to this maximum, all of the measured surface, middle, and bottom currents flowed along the wind direction and their speed significantly increased (Fig. 7). The strong southeast-east winds drive the entire water column at the LHWS to form wind-driven currents, thus resulting in a downwind upwelling (Wu et al., 2018). Similarly, generated by the strong north-northeast wind event during the 2018 winter observation, wind-driven currents also resulted in downwind upwelling 315 (Fig. 11). These upwelling processes provided an excellent opportunity to evaluate the performance of the WCCM in Lake Taihu.
The numerical solutions of the governing equations and most parameter values of the WCCM are similar to those of the EFDC. The main differences between the two models are the vertical eddy viscosity, wind drag coefficient, and waveinduced radiation stress. The numerical experiments show that the average correlation coefficient between the WCCM-320 simulated and measured current speeds increased by 36.4% compared with the LCM_1 results, or 42.9% compared with the EFDC results in 2018. The current direction and field simulated by the WCCM also improved, whereas the water level was simulated at a similar accuracy as by the EFDC. Compared with the reference model, the WCCM is more reliable to simulate wind-driven currents and subsequent downwind upwelling in Lake Taihu. The WCCM can also accurately simulate wind waves and water temperature in the lake (Fig. B.1 and B.2).

Wind drag coefficient
The wind drag coefficient 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 the current velocity in Lake Taihu (Li et al., 2015). Our numerical experiments also indicate that the correlation coefficients between the simulated and measured current speeds of LCM_2 and WCCM are significantly greater than those of EFDC and LCM_1 (Tables 3, 5). This implies 330 that the new expressions of C s (Eqs. (16) and (17)) mainly contributes to the enhanced correlation coefficients. Based on previous studies (Edson et al., 2013;Geernaert et al., 1987;Large and Pond, 1981;Xiao et al., 2013) and our field observations, these expressions were derived to describe the discontinuity of changing trend of C s with wind and directionality of the wind momentum transmission.
The magnitude of C s represents the transmission efficiency of the wind momentum to a waterbody and its change is 335 discontinuous. Surface waves can increase the roughness of a lake surface, and further influence the transmission efficiency (Xiao et al., 2013). The transmission efficiency on aerodynamically rough water surfaces is higher than that on aerodynamically smooth water surfaces (Lükő et al., 2020). Wu (1980) proposed that the atmospheric surface layer appears to be aerodynamically rough when the wind speed exceeds 7.5 m s −1 . This implies that there is a discontinuity of the C s curves at wind speeds of 7.5 m s −1 . Field observations in Lake Taihu (Xiao et al., 2013) indicate that the measured C s 340 initially decreased under light-wind conditions (< ~7.5 m s −1 ) and then increased under high-wind conditions (> ~7.5 m s −1 ).
The curves of C s plotted by the equation proposed by Edson et al. (2013) and Large and Pond (1981) also intersect at a wind speed of 7.5 m s −1 (Fig. 4). A wind speed of 7.5 m s −1 is therefore reasonable for defining the discontinuity of changing trend of C s with wind in Lake Taihu.
As shown in Eqs. (16) and (17), a logistic curve is used to describe the increase of C s for wind speeds > 7.5 m s −1 ; otherwise 345 C s is a constant. Under light-wind conditions, the mechanism of the C s change with wind speed remains incompletely understood and its mathematic description is non-deterministic (Fig. 4). According to a tremendous amount of measured C s values reported by Edson et al. (2013), the points between C s and wind speed evenly distribute on both sides of a constant under light-wind conditions. A constant is therefore suitable (Large and Pond, 1981). Under high-wind conditions, the proportional function is most frequently used to fit the C s change (Geernaert et al., 1987;Large and Pond, 1981;Wu et al., 350 1981;Zhou et al., 2009). However, the measured C s values indicate more rapid changes than described by the proportional function (Edson et al., 2013). Furthermore, Geernaert et al. (1987) concluded that C s increases to a constant (~0.003) by compiling all of the reported C s measurements. The logistic function is therefore used to fit the rapid increase that then tends toward a constant. It should also be noted that the curves of Eqs. (16) and (17) and w s × Eq. (21) used in this study are significantly lower than the other two curves (Fig. 4). The main cause is that the limited water depth and fetch in Lake Taihu 355 reduce the transmission efficiency of the wind momentum and restrict the development of wind-driven currents in the lake.
The directionality of wind momentum transmission is further addressed using different C s values in the x-and y-directions.
There have been numerous expressions designed to calculate the wind drag coefficient based on ocean environments without consideration of the directionality of wind momentum transmission (Geernaert et al., 1987;Large and Pond, 1981;Lükő et al., 2020;Wu, 1980;Zhou et al., 2009). However, the increase of transmission efficiency with wind speed (Lükő et al., 2020) 360 will result in a contradiction in these existing expressions that the same C s values are used in x-and y-direction while the components of wind speed in these directions are different. Moreover, wind waves and lake seiche also have directionality, which can affect the transmission efficiency of the wind momentum by changing the roughness and tilt of the lake surface.
Neglecting the directionality of wind momentum transmission can therefore over-or under-estimate the wind drag coefficient in any one direction in large shallow lakes. 365

Wave-induced radiation stress
Wave-induced radiation stress is first considered in simulating wind-driven currents in large shallow lakes. The results show that this consideration can improve the simulated current direction. The MAE UVD values of the LCM_2 (average MAE UVD of 56.3°; Table 3) in 2015 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 370 coefficients of LCM_2 in 2018 are slightly lower than those of the WCCM in 2018 (Table 5), which implies that waveinduced radiation stress can also contribute to the improvement of the WCCM-simulated current speed.
A comparison between the WCCM-and EFDC-simulated current fields further demonstrates the importance of waveinduced 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 375 EFDC exhibit counterclockwise vortexes in Zhushan Bay and Meiliang Bay (Fig. B.3), 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 turbulent owing to wave deformation resulting from the shallow water and lakeshore. The wave-induced radiation stress therefore reduces the likelihood that a vortex will form in this area. Conversely, the middle and bottom current fields simulated by the LCM_2 without wave-induced radiation stress 380 also show counterclockwise vortexes in Zhushan Bay and Meiliang Bay (Fig. B.5), which is similar to the EFDC result. It is very important for Lake Taihu that the absences of vortexes in the downwind area will reinforce the accumulation of buoyant cyanobacteria, and further promote cyanobacterial blooms within this area.

Vertical eddy viscosity
The vertical eddy viscosity play a less prominent role in the development of wind-driven currents than the other variables. In 385 this study, the 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 (Eqs. (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-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 efficient.

Challenges of the hydrodynamic model development for shallow lakes
Although the WCCM performance has been improved relative to the reference models of the EFDC, LCM_1, and LCM_2, the correlation between WCCM-simulated and ADP-measured current speed remains low, and the mean of the simulated current speed is lower than that of the measured current speed. Similar conclusions can be drawn from the model validation studies in other lakes (Huang et al., 2010;Jin et al., 2000;Ishikawa et al., 2021;Soulignac et al., 2017). There are three 395 possible explanations for this problem. (1) Based on the Doppler effect of sound waves, the ADP measures the threedimensional lake currents by detecting the movement of suspended particle matter (SPM) in water column. However, the spatiotemporal distributions of the concentration and physicochemical properties of the SPM are changeable in lakes (Zheng et al., 2017). This will undoubtedly influence the measurements of real currents in lakes.
(2) The spatiotemporal resolution of the numerical model input data can introduce errors into the lake current simulations, including mesh, underwater 400 topography, boundary conditions, and wind field. (3) The wind-induced hydrodynamics in large shallow lakes are not fully understood. For example, Eqs. (16) and (17) derived from the field observations are only effective when the wind speed is ≤ 15.5 m s −1 , which is the maximum of the field observations, meanwhile the contributions of the wind waves to the development of wind-driven currents are underestimated in Lake Taihu.

Conclusion 405
Strong summer or winter winds generate wind-driven currents in Lake Taihu, which subsequently results in downwind upwelling events. The WCCM is developed that reconsiders the expression of wind stress, wind waves, and turbulence based on these events and numerical experiments. This model can simulate the development of wind-driven currents with a 42.9% increase of simulated current speed compared with the EFDC results of 2018. The new expression for the wind drag coefficient is mainly responsible for increasing the correlation coefficient between the WCCM-simulated and measured 410 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 current speed simulation. The simple parameterized turbulence scheme is sufficient for simulating wind-driven currents in Lake Taihu. We emphasize that more process-based field observations using advanced instruments are required to fully understand the real hydrodynamic characteristics of large shallow lakes and further improve the performance of lake hydrodynamic models, especially for the interaction between wind waves and lake 415 currents.

Code and data availability
The source code of the EFDC model is freely available from https://doi.org/10.5281/zenodo.5602801 (Wu, 2021