Evaluation of the WRF lake module (v1.0) and its improvements at a deep reservoir

. Large lakes and reservoirs play important roles in modulating regional hydrological cycles and climate; however, 10 their representations in coupled models remain uncertain. The existing lake module in the Weather Research and Forecasting (WRF) system (hereafter WRF-Lake), although widely used, did not accurately predict temperature profiles in deep lakes mainly due to its poor lake surface property parameterizations and underestimation of heat transfer between lake layers. We therefore revised WRF-Lake by improving its (1) numerical discretization scheme; (2) surface property parameterization; (3) diffusivity parameterization for deep lakes; and (4) convection scheme, the outcome of which became WRF-rLake (i.e., revised 15 lake model). We evaluated the off-line WRF-rLake by comparing simulated and measured water temperature at the Nuozhadu Reservoir, a deep reservoir in southwestern China. WRF-rLake performs better than its predecessor by reducing the root-mean-square-error (RMSE) against observed lake surface temperatures (LSTs) from 1.4 °C to 1.1 °C and consistently improving simulated vertical temperature profiles. We also evaluated the sensitivity of simulated water temperature and surface energy fluxes to various modelled lake processes. We found (1) large changes in surface energy balance fluxes (up to 60 W m -2 ) 20 associated with the improved surface property parameterization and (2) that the simulated lake thermal structure depends strongly on the light extinction coefficient and vertical diffusivity. Although currently only evaluated at the Nuozhadu Reservoir, we expect that these model parameterization and structural improvements could be general and therefore recommend further testing at other deep lakes


Introduction
Inland waters such as lakes and reservoirs differ from their surrounding land with altered albedo, larger thermal conductance and heat capacity, and lower surface roughness, and therefore different radiative and thermal properties.These lake properties can exert significant influences on local and regional climate and are important in understanding lakeatmosphere interactions (Hostetler et al., 1994;Bonan, 1995;Lofgren, 1997;Krinner, 2003;Long et al., 2007;Samuelsson et al., 2010;Dutra et al., 2010;Subin et al., 2012;Deng et al., 2013;Xiao et al., 2016).
Lakes interact with atmosphere through energy, mass (mostly water), and momentum exchanges (Lerman et al., 2013).Generally, due to their larger thermal inertia and smaller roughness (Samuelsson et al., 2010;Xiao et al., 2016), lakes tend to attenuate surface diurnal temperature variation and enhance surface wind compared to the surrounding land.The influences of lakes on regional climate vary by season (Subin et al., 2012).In the early winter and spring, lakes warm and moisten overlying air masses, generating the so-called "lake effect" on precipitation manifested as heavier snowfall in downwind regions (Bates et al., 1993;Niziol et al., 1995;Scott and Huff, 1996;Zhao et al., 2012;Wright et al., 2013).In the early summer, reduced sensible heat fluxes into the atmosphere are observed in some temperate and high latitude lakes because of their lower surface temperature and smaller roughness than the surrounding land (Lofgren, 1997;Krinner, 2003;Dutra et al., 2010).Also, anticyclones (cyclones) tend to be intensified in summer (winter) through their interactions with the Great Lakes (Notaro et al., 2013;Xiao et al., 2016).During fall and early winter, when the lake surface is warmer than the overlying air, high-latitude lakes (e.g., the Great Bear Lake in Canada) release the heat collected during summer to the atmosphere, reducing snow accumulation in the surface areas around the lakes (Long et al., 2007).
Since lakes strongly affect lower atmospheric heat, water, and momentum states and fluxes, predicting weather and climate in lake basins or lake-rich regions requires realistic lake representations in numerical weather prediction models (NWPMs) and climate models.In the last decade, many efforts have been made focusing on the development and analysis of lake modules in such models (MacKay et al., 2009;Bonan, 1995;Mironov et al., 2010;Dutra et al., 2010;Stepanenko et al., 2010;Subin et al., 2012;Subin et al., 2013;Stepanenko et al., 2013).Although there exist sophisticated lake models that account for detailed spatiotemporal dynamics of various lake processes, it is not a common practice to fully couple atmospheric models with them because of their high computational cost and prohibitive complexity for coupling (MacKay et al., 2009).
In contrast, one-dimensional (1D) models have been widely used for coupling with atmospheric models, because 1D models are sufficiently fast to facilitate long-term coupled simulations and have performed well in simulating seasonal and interannual variations of lake water temperature (Peeters et al., 2002).Depending on how they parameterize eddy diffusivity, these 1D lake models mainly fall into two categories (MacKay et al., 2009;Perroud et al., 2009;Martynov et al., 2010;Stepanenko et al., 2010; Table 1): the Hostetler-type models (e.g., the Hostetler Model (Hostetler and Bartlein, 1990;Hostetler et al., 1994), Minlake (Fang and Stefan, 1996), SEEMOD (Zamboni et al., 1992), LIMNMOD (Karagounis et al., 1993), MASAS and CHEMSEE (Ulrich et al., 1995), CLM4-LISSS (Subin et al., 2012), and WRF-Lake (Gu et al., 2015)) and the more sophisticated turbulence models (e.g., the bulk model of Kraus and Turner (1967), DYRESM (Imberger et al., 1978), PROBE (Svensson, 1978), GOTM (Burchard et al., 1999), SIMSTRAT (Goudsmit et al., 2002), and LAKE (Stepanenko and Lykosov, 2005;Stepanenko et al., 2011)).The Hostetler-type models use parameterized eddy diffusivity to model vertical mixing in the lake body.The eddy diffusivity is dependent on surface wind speed and lake stratification and its parameterization follows that of Henderson-Sellers (1985), which is formulated under the assumptions of unstratified Ekman flow (Smith, 1979).Although the representativeness of this eddy diffusivity scheme for lakes has not been systematically evaluated, these models have been applied in numerous coupled simulations with regional climate models (Bates et al., 1995;Hostetler and Giorgi, 1995;Small et al., 1999) and general circulation models (Bonan, 1995;Krinner, 2003).The more complex turbulence models use the - turbulence scheme and parameterize eddy diffusivity based on the Kolmogorov-Prandtl relation.Thus, the turbulence models require two additional equations for the turbulent kinetic energy () and its dissipation rate ().Although models from both categories have been intensively used in climate models, they share the potential and common drawback that the eddy diffusivity representations may be inappropriate when temperature gradients are weak (MacKay et al., 2009).
In addition to the above two categories, bulk mixed layer models have also been developed, including FLake (a relatively simple 2-layer model based on the similarity theory (Mironov et al., 2010)) as a typical example.Gula et al. (2012) and Mallard et al. (2014) have coupled FLake to WRF in 1-way and 2-way model configurations, respectively.However, it is difficult for these oversimplified lake models to capture seasonal stratification and to accurately simulate water temperature in deep lakes.Thus, the performance of the bulk mixed layer models in climate-modelling studies is still limited.WRF-Lake is a 1-D lake model that has been widely applied to study how lakes affect local weather and regional climate (Gu et al., 2015;Mallard et al., 2015;Gu et al., 2016;Xiao et al., 2016).It is an eddy-diffusion model adapted from the Community Land Model (CLM) version 4.5 (Oleson et al. 2013) with further modifications by Gu et al. (2015).However, some previous studies and our evaluation here suggest that WRF-Lake is inaccurate when simulating deep lakes and reservoirs (Mallard et al., 2015;Xiao et al., 2016).To better represent thermodynamic processes of deep lakes and reservoirs in regional climate simulations, we here thus revised the WRF-Lake model.Our revisions fall into four categories: (1) including an improved spatial discretization scheme; (2) improving the surface property parameterization; (3) improving the diffusivity parameterization for deep lakes; and (4) improving the convection scheme to avoid unphysical mixing.We evaluated the improved lake model, WRF-rLake (i.e., revised lake model), at the Nuozhadu Reservoir, a deep reservoir (up to 200 m deep) in southwestern China and conducted sensitivity experiments to evaluate the effect of dominant processes and parameters on simulated lake water temperature and surface heat fluxes.

Model Name Main Features References
Hostetler-type Model 2 Model description and improvements

Description of the current WRF-Lake model
The WRF-Lake model is a mass and energy balance scheme which vertically divides the lake into 20-25 layers and solves the 1D heat diffusion equation.The lake includes 0-5 snow layers; 10 lake liquid water and ice layers (hereafter lake body); and 10 sediment layers (Subin et al., 2012).The lake body water content is assumed to be constant and sediment layers are fully saturated.An infinitely small interface is assumed between the first lake layer and the overlying atmosphere to calculate lake surface fluxes of heat, water mass, momentum, and radiation.After subtracting latent and sensible heat fluxes from the surface net all-wave radiation, the residual energy fluxes are set as the top boundary condition to force the heat diffusion equation.Mixing processes in the lake body include molecular diffusion, wind-driven eddy diffusion, and convective mixing (Hostetler and Bartlein, 1990).We describe below the basic model structure to indicate model processes and parameters that were analyzed and modified in this study.

Surface energy balance
In the WRF-Lake model, the energy balance at the lake surface is: where  is absorbed shortwave radiation,  is net longwave radiation, and  is downward heat flux into the lake,  is sensible heat, and  is latent heat.The units for all variables in equation ( 1) are W m -2 .In this equation,  is regulated by albedo (); , , and  are largely dependent on the thermal conductivity of the top lake layer (κ), aerodynamic resistance for heat ( ℎ ) and that for vapor (  ), respectively.A more thorough discussion of the relationship between lake surface fluxes and aerodynamic resistances is provided by section 2.1.8 in Subin et al. (2012).

Radiation transfer and absorption in the lake body
For unfrozen lakes, a fraction () of the incoming solar radiation is absorbed by the surface water within the first 0.6 m.The remaining solar radiation is absorbed by the lake body following Beer's law: where  (W m -2 ) is the distribution of solar radiation with depth  (m; positive downward),  is set to a constant 0.4 (Oleson et al., 2010),   is the base of the surface absorption layer (0.6 m), and  is the light extinction coefficient (m -1 ) which is calculated as a function of the lake depth according to Hakanson (1995): where  (m) is the lake depth.
Though there exist more sophisticated radiation schemes in other lake models (e.g., the 9-band scheme by Paulson and Simpson, 1981), we kept the current WRF-Lake radiation scheme since it includes the essential components in a waterbody physical radiation parameterization: an intensity-decaying formulation as a function of penetration depth following the Beer-Lambert law (Jerlov, 1976) and a scheme for absorption coefficients.Such an approach is also accepted by many other 1D lake models (Fang and Stefan, 1996;Stepanenko and Lykosov, 2005).To improve the model performance, we tentatively set the cutoff depth   to be 0 m in this version as 0.6 m is usually an overestimated value, especially for shallow lakes (Deng et al., 2013).Although adopting this   value (0 m) demonstrates acceptable performance in this work, a more lake-specific cutoff depth may be needed for better model performance.

The heat diffusion solution
After the surface energy balance and radiation absorption are calculated, the following 1D thermal diffusion equation is solved: where   is the volumetric heat capacity of water (J K -1 m -3 ),  is water temperature (K),  is time (s),  is depth (m),  is the thermal diffusivity (m 2 s -1 ), and  is the solar radiation as in equation ( 3).
The thermal diffusivity  in WRF-Lake comprises molecular diffusivity (  ) and wind-driven eddy diffusivity (  ).
is a function of thermal conductivity and volumetric heat capacity of water (1.433 × 10 −7 m 2 s -1 ).  follows the Henderson-Sellers parameterization (shown in Appendix) and is determined by wind speed at 2 m above the water surface, latitude dependent Ekman profile, and lake stratification dependent Brunt-Väisälä frequency (Henderson-Sellers et al., 1983;Henderson-Sellers, 1985).However, for deep lakes, e.g., Lake Michigan (Martynov et al., 2010), Subin et al. (2012) argued that increasing the eddy diffusivity by a factor of 10, and in some cases by a factor of 100, substantially improved simulated lake water temperatures and surface fluxes.Thus, in WRF-Lake, the eddy diffusivity for lakes deeper than 15 m was multiplied by a factor ranging from 10 2 to 10 5 depending on lake depth and surface temperature (Gu et al., 2015).

Convective mixing
The convective mixing module is executed after the heat diffusion equation solution and is identical to that of Hostetler's lake model, which assumes that no temperature instability can be sustained in the lake body.In particular, at the end of any numerical time step, if there is denser water overlying lighter water, then convective mixing happens rapidly, mixing the lake body from surface to the unstable layers (Hostetler and Bartlein, 1990).

Improvements by WRF-rLake
As described below, we modified the existing WRF-Lake model (and named the revision as WRF-rLake) by improving the spatial discretization scheme, parameterization for surface properties and vertical diffusivity, and the convection scheme (Figure 1).

Vertical discretization
The WRF-Lake model discretizes the water body into 10 layers with the top-most layer fixed to 0.1 m (Gu et al., 2015) and each of the other nine layers constituting 10% of the total depth.Although this discretization works well with shallow lakes (e.g., depth < 50 m), it may be problematic for deep lakes in three aspects.The first is insufficient resolution.
Discretizing lakes into only 10 layers is too coarse for the 200-m-deep Nuozhadu Reservoir, resulting in only several layers to cover the thermocline.The second is depth loss, because the 9 layers below the first layer only account for 90% of the lake body, which, for very deep lakes, may cause up to approximately 10% lake depth loss and will potentially lead to energy loss of the simulated lake system.The third is the large layer thickness transition between the first two layers.In the case of the Nuozhadu Reservoir (~200 m deep), the default discretization results in a sharp thickness increase from 0.1 m to 20 m (200 times) between the top two layers, which may be numerically problematic.
We therefore introduced a new discretization scheme for the lake body where both the vertical resolution and layer depth distribution are improved.For lakes deeper than 50 m, a 25-layer adaptive discretization scheme is used where the topmost layer is set to 0.1 m and the remaining layers have their thicknesses increasing exponentially by a fixed factor (FF) that depends on lake depth (Table 2).For the Nuozhadu Reservoir, FF is taken to be 1.29, which results in a series of layer depths (m) from the top to the bottom of: 0.1, 0.1, 0.2, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.3, 1.6, 2.1, 2.7, 3.5, 4.6, 5.9, 7.6, 9.8, 12.6, 16.3, 21.0, 27.1, 35.0, 44.7 m.Similar adaptive discretization techniques are widely used by various geoscientific models to improve model performance, e.g., the grid stretching in GFDL HiRAM (Harris et al., 2016) and nonuniform meshing in MPAS (Skamarock et al., 2018) etc.When a lake is less than 50 m deep, a new 10-layer discretization scheme is applied where the top layer is fixed at 0.1 m and the remaining depth is allocated evenly among the remaining nine layers.

Surface properties and parameters
As discussed in section 2.1.1,the aerodynamic resistances for heat ( ℎ ) and vapor (  ) heat fluxes are critical for surface energy balance predictions.The aerodynamic resistances are functions of momentum ( 0 ) and scalar roughness lengths ( 0ℎ for sensible heat and  0 for latent heat).In WRF-Lake,  0 is set to 1 mm, 5 mm and 2.5 mm for unfrozen lakes, frozen lakes without snow, and frozen lakes with snow, respectively; while  0ℎ and  0 are always kept equal to  0 .
However, for unfrozen lakes, the momentum roughness length of 1 mm is often too large.Open seas are believed to have larger roughness lengths than lakes due to better developed surface moving waves, while the Engineering Sciences Data Unit (ESDU) (1972) documentation suggested  0 = 1 mm for normal and 0.1 mm for calm seas.A number of lake studies also have shown 1 mm to be a maximum  0 value.For example, measurements on Lake Washington in the U.S. show  0 is generally below 1 mm and ranged between 0.01 -1 mm (Atakturk and Katsaros, 1999); measurements on Lake Ngoring, a high-altitude lake in the Tibetan Plateau, found  0 ranged between 0.001~1 mm and seldom went beyond 1 mm (Li et al, 2015).Thus, in order to produce more realistic roughness lengths for lakes, CLM4-LISSS parameterized  0 based on the forcing wind, friction velocity, fetch, and lake depth.Adopting these parameterizations has produced more accurate surface heat fluxes and LSTs over many natural lakes (Subin et al., 2012).
We therefore adopted the CLM4-LISSS parameterization of roughness lengths with some further modifications: •  0 is set to 2.4 mm for frozen lakes with snow and 1 mm for frozen lakes without explicit snow layers; and  0ℎ and  0 are computed as functions of  0 and friction velocity u * (m s -1 ).
• For unfrozen lakes,  0 is parameterized as: where γ is a dimensionless empirical constant (0.1),  is the dimensionless Charnock coefficient described below, and  is the acceleration of gravity (Fairall et al., 1996;Charnock, 1955;Smith, 1988): where  (m) is the lake fetch (assumed to be 25 times of lake depth when observations are unavailable),  (m s -1 ) is 2 m wind speed,  (m) is lake depth,   = 0.01,   = 0.11, and  = 1. and  account for influences from fetch and depth, respectively.
We also made further improvements in roughness lengths parameterizations for unfrozen lakes.In LISSS,   should be 100 but was tentatively set to 22, corresponding to the use of  instead of u * in equation ( 7).Subin et al. (2012) recommended that future lake models should relax this assumption by setting   = 100 and directly applying u * rather than , which we have done in WRF-rLake.Because u * depends on surface roughness lengths, which in turn depend on u * , we introduced a fixedpoint iteration for the equations relating u * and surface roughness lengths.It is worth noting that the parameterization of roughness lengths for frozen lakes could also be improved.However, as the Nuozhadu Reservoir is unfrozen throughout the year, we did not make any modifications to the representations for lake ice.Future work should investigate lakes with frozen periods to further improve the roughness length parameterization.

Improved diffusivity parameterization
The Henderson-Sellers parameterization (see Appendix) for wind-driven eddy diffusivity underestimates mixing in deep lakes and was therefore increased in WRF-Lake by a multiplicative factor (Gu et al., 2015).However, this treatment may trigger new problems.As   declines exponentially with depth (equation A1), it is more likely to be underestimated in deeper layers than in the topmost one.Thus, enlarging   by the same factor for the whole lake may introduce new problems in two ways.
First,   may be overestimated in lake surface layers.A number of empirical studies have estimated the effective vertical heat diffusivity in lakes and coastal oceans using heat flux measurements and tracer distribution measurements.These studies showed that vertical heat diffusivity in natural lakes seldom exceeds ~1 cm 2 s -1 and in oceans, where the forcing winds are usually stronger and moving surface waves are better developed, vertical heat diffusivity exceeds ~10 2 cm 2 s -1 (Hutchinson, 1957;Li, 1973;Kullenberg et al., 1973;Kullenberg, 1974;Jassby and Powell, 1975;Sarmiento et al, 1976;Quay et al., 1980).
We thus set a maximum of 10 2 cm 2 s -1 for wind-driven eddy diffusivity (corresponding to the maximum value for open seas) to avoid overestimation in surface layers.
Second, for deep layers,   may still be underestimated by WRF-Lake because   is forced to decline exponentially with depth.Additional diffusion terms should be included to account for turbulent mixing generated by other mechanisms in deeper layers where wind driven eddies cannot penetrate.So, in WRF-rLake, we adopted the enhanced diffusion term,   (m 2 s -1 ), from CLM4-LISSS (Ellis et al., 1991;Fang and Stefan, 1996;Subin et al., 2012), which was originally suggested by Fang and Stefan (1996) to compensate for unresolved turbulence (even below ice or at large depth) and is given as follows: where  (s -1 ) is the Brunt-Väisälä frequency (a minimum  2 is suggested to be 7.5 × 10 −5 s -2 , see Fang and Stefan (1996)).
However, as discussed by Subin et al. (2012), this suggested value for   is only several times larger than   and may still be too small for deep lakes.Therefore, for lakes deeper than 50 m, we imposed an increase in   by a factor of 100 for all layers and argue that more analyses are required to robustly represent unresolved turbulence.

Convective mixing
In WRF-Lake, density instability is the prerequisite for convective mixing.As we have mentioned above, density for each water layer is first calculated based on lake water temperature, then adjacent layers are compared for their densities to determine whether there should be convective mixing.However, since the water density is calculated as a function of temperature, when the temperature gradient between two adjacent layers is small (usually less than 10 -3 K m -1 , but still with a lighter layer overlying a heavier layer), small numerical errors may incorrectly trigger convective mixing.In WRF-rLake we therefore set a density gradient threshold of 10 -4 kg m -3 m -1 (which is equivalent to a temperature gradient threshold of about 10 -3 K m -1 ) to avoid this unphysical convective mixing.In our tests at Nuozhadu Reservoir, convective mixing can penetrate as deep as 5 meters below the surface in WRF-Lake, while it is almost always restricted within the top 1 m in WRF-rLake.

Study area
The Nuozhadu Reservoir is near the downstream end of a group of reservoirs along the Lancang River (22˚38' N, 100˚26' E), which is located in southwestern China and is called the Mekong River when it leaves China and enters Laos (Figure 2).The Nuozhadu Reservoir has a dam height of 262 m and a normal water level1 of 812 m above sea level.The water depth upstream of the dam in year 2015 is around 200 m.The water surface area has increased more than 10 times after construction of the reservoir.Owing to its particular location, great depth, and large surface area, the Nuozhadu Reservoir serves as a good example for research on the impacts of artificial inland waters on regional climate.

Forcing data
Our study period covers from January 1, 2015 to December 31, 2015, a year when the reservoir was under normal operation.We ran the lake module off-line, driven directly by forcing data acquired from local meteorological stations rather than WRF-simulated fields, in order to evaluate the lake module free from potential biases originating in WRF.The forcing data of the first day was repeated 7 times to form a one-week spin-up.The downward shortwave radiation and downward longwave radiation were obtained from the China Meteorological Forcing Dataset (Kun et al, 2010; Chen et al, 2011) with a temporal resolution of 3 hours and spatial resolution of 0.1°  0.1°.A linear interpolation was applied to these data to obtain hourly forcing for WRF-rLake.Although it probably underestimates peak radiation values, linear interpolation may still be 5 considered to be an acceptable approximation given no data of higher temporal resolution is available.
Other forcing data, including atmospheric temperature, atmospheric pressure, atmospheric specific humidity, atmospheric wind speed in the east and north directions, and precipitation, were acquired with a 1 hour temporal resolution from a meteorological observatory run by China Huaneng Group Co., Ltd., the construction unit of the Nuozhadu Reservoir (Figure 3).The station (22°40′N, 100°23′E) is located about 5 km upstream (or northwest) of the Nuozhadu Dam with an 10 observational height of 10 m.

Observations for model initialization and evaluation
Water temperature was measured hourly from 712 m to 804 m above sea level at an interval of 2 m near the Nuozhadu 5 dam during 2015 (Figure 4).Since the reservoir is in a tropical region, surface water temperature is higher than 20 ℃ throughout the year.In 2015, the water level was dropped from January to June in preparation for the rainy season and rose gradually thereafter from July to December.
Measured water temperatures on 1st January 2015 was used to initialize the first 90 m of the water body.The remaining 110 m of depth (for which observations were not made) were initialized using simulated results at the same site on the same day by Delft3D FLOW (Chen, 2017), a three-dimensional hydrodynamic model proven to be sufficiently accurate at the Nuozhadu Reservoir.

Numerical experiments
To examine the incremental improvements in the WRF-rLake simulations, we performed four sets of off-line numerical experiments analyzing four key parameterizations (Table 3): 1) Vertical discretization ("Lyr" set): the default WRF-Lake 10-layer and modified 25-layer settings were contrasted to assess the impacts of different vertical discretization schemes.
2) Vertical Diffusivity ("Diff" set): we examined the impact of vertical diffusivity by applying different diffusivity schemes: the original scheme by Hostetler and Bartlein (1990), the scheme by Gu et al. (2015), the scheme by Gu et al. (2015) with enhanced diffusion term, and the scheme as discussed in section 2.2.3 or called "modified" diffusivity as adopted by our new model WRF-rLake.
3) Roughness length ("Rou" set): momentum and scalar roughness lengths are set to 1 mm, 10 mm, calculated as in Subin et al. (2012), or calculated with further modification as in section 2.2.2, to examine the effects of different roughness lengths schemes.4) Light extinction coefficient ("Ext" set): through model tests, we conclude that in addition to the schemes we modified, the light extinction coefficient is also a key parameter for accurately modelling deep lakes (Hocking and Straskraba, 1999).Although the default parameterization of light extinction coefficient has been applied in previous WRF-Lake studies (e.g., Gu et al., 2015), we tested the impacts of different values of this coefficient.Given that no Secchi disk measurements were available in our study site, no empirical constrains for the light extinction coefficient could be directly developed.Thus, we tested a range of light extinction coefficient values: 0.13 m -1 (default), 0.30 m -1 , 1.00 m -1 , and 3.00 m -1 .Although measurements have reported larger variability of light extinction coefficient (e.g., 0.05 to 7.1 m -1 in Subin et al. ( 2012)), we found simulated temperature profiles were insensitive to values outside of the 0.13 to 3.0 m -1 range.We concluded that the best performance could be achieved by increasing the light extinction coefficient to ~1.00 m -1 , which thus is adopted in our baseline run (BL).
In addition, a control run (CTL) was configured with default roughness lengths, extinction coefficient, and vertical diffusivity (Gu et al. 2015) as in the default WRF-Lake; and a baseline run (BL, i.e., our proposed new model structure) was configured with all modifications described in section 2.2 and a tuned light extinction coefficient to demonstrate the effects of each improvement in WRF-rLake.#: as the reservoir was ice-free in year 2015, the constants given in this column refer to the roughness lengths for unfrozen lakes.

Comparison of simulated temperature fields between WRF-Lake and WRF-rLake
The simulation results near the dam, the same place where the observations were collected, were used to conduct the evaluation.In comparing WRF-Lake simulations (CTL) to the observations (Figure 5), the LSTs are underestimated most of the time with a Mean Bias Error (MBE) of -0.69 °C.The largest bias is found in the 1st quarter where it reached -3.37 °C.This bias mainly resulted from the overestimation of roughness lengths by prescribing them to 1 mm, which resulted in too large outgoing latent and sensible heat fluxes (section 4.2).The simulation by BL is better in terms of RMSE of simulated LSTs, which is reduced to 1.14 °C from 1.35 °C by CTL.Vertical temperature profiles were also improved, as the thermocline in the top 10 meters is reproduced in hot seasons, in contrast to the CTL and Diff_3 simulations.More detailed statistical metrics of the vertical temperature profile (Table 4) suggest that BL gives the best simulations among the three simulations compared.Applying all the modifications to vertical diffusivity (discussed in section 2.2.4), we obtained the best simulations by BL; however, we note that the original diffusivity parameterization of Henderson-Sellers (1985) may be inappropriate for deep lakes as this parameterization finds it hard to yield good diffusivity simulations for deep lakes without tuning.Thus, we suggest more thorough evaluation and modification to this parameterization should be carried out in future research.

Effects of vertical discretization
In these experiments, Lyr_1 adopted the default WRF-Lake 10-layer scheme and was identical to BL except for discretization.Lyr_2 (identical to BL) uses the modified 25-layer scheme.Overall, compared to Lyr_1, Lyr_2 reduced the RMSE against monthly observed lake temperatures profiles from 1.64 °C to 1.13 °C.Lyr_1 predicted too much vertical mixing in the top 10 m, failing to reproduce the thermocline in this zone in summer (Figure 7).In the first two months, the temperature difference in top 10 m between Lyr_2 and Lyr_1 was not obvious.But as the lake warms up, the thermocline in Lyr_2 (top 10 m) is intensified, hindering heat transfer to the underlying layers, which in turn further strengthened the thermocline.However, this process is not captured by the Lyr_1 model; rather, a well-mixed water column in the top 10 m (top two layers) is simulated throughout the year.From March to June, the overmixing by Lyr_1 also resulted in colder LSTs, which in turn produced lower sensible and latent heat fluxes by up to 10 W m -2 and 20 W m -2 , respectively (not shown).Lower outgoing heat fluxes kept more energy stored in the lake body, explaining why the whole lake body became increasingly warmer throughout the summer compared to Lyr_2.This underestimation of outgoing heat fluxes was reversed when the surface temperature in Lyr_1 finally exceeded that of Lyr_2 and produced more sensible and latent heat fluxes to the atmosphere from October to December.
In summary, the 10-layer scheme produced more uniform temperature fields across the water column and, due to its coarser spatial discretization, poorly predicted the thermocline where temperature changes rapidly.We further tested a 100layer scheme but its simulations did not differ significantly from the 25-layer scheme (not shown here), so we conclude that the 25-layer scheme is sufficient for deep lakes like the Nuozhadu Reservoir.For the monthly-averaged vertical temperature profile (Figure 8), the BL scheme best captures the decrease of water temperature with depth and LSTs.BL yields the smallest RMSE of 1.13 °C against monthly observed lake temperatures profiles, while Diff_1, Diff_2, and Diff_3 yield 1.62 °C, 1.47 °C, 1.47 °C, respectively.Diff_1 yields strong stratification 10 within 10 meters of the surface, indicating that the eddy diffusivity by Hostetler and Bartlein (1990) is too small and prevents heat from transferring from the surface to depth.This suppression reduces heat stored in the lake body and weakens the thermal inertia of the lake, leading to biased high LSTs in summer and biased low LSTs in winter.Diff_2 and Diff_3 produce identical vertical temperature profiles, indicating that the enhanced diffusion term is relatively small compared to eddy diffusivity and makes little difference.Diff_2 and Diff_3 produce lake layers of almost 5 uniform temperatures in the top 10 m (Figure 8).This pattern is further confirmed by the strong vertical mixing in the first 10 m inferred by the vertical profile of overall diffusivity that includes molecular, eddy, and enhanced diffusivity (Figure 9).The overall diffusivity by Diff_2 in the top 10 meters could be as large as 10 3 cm 2 s -1 , which is too large compared to estimates in the literature (section 2.2.3), supporting the necessity for limiting eddy diffusivity (section 2.2.4).Li (1973).
The overall vertical diffusivity below 20 m by Diff_2 or Diff_3 is of the same magnitude as molecular diffusivity (1.43 × 10 −3 cm 2 s -1 ).By increasing the enhanced diffusion term by 100 times for deep lakes, BL yielded more reasonable vertical diffusivity below 20 m (Figure 9) which is quite close to an estimation at a very similar lake: Li (1973) estimated the vertical diffusivity in Lake Zürich, a lake with similar topography and depth (~130 m deep) to the Nuozhadu Reservoir and concluded its value for vertical diffusivity mostly ranged between 0.1 ~ 1 cm 2 s -1 .
With the constrained eddy diffusivity and "enhanced diffusion term", BL produced the best temperature profiles compared to observations (Figure 8).Further, the total diffusivity affects both the shape of vertical temperature profiles and the quantity of energy transferred down from the surface, which further influences LSTs.

Effects of roughness lengths
We next compare BL, Rou_1, Rou_2, and Rou_3, experiments that were conducted with different roughness length parameterizations, where the latter three have roughness lengths of 1 mm, 10 mm, and the parameterization from the Subin et al. (2012) scheme.In BL, the value of roughness values varied but were almost always smaller than Rou_1 and Rou_2 (generally less than 0.5 mm).BL and Rou_3 produced almost identical LSTs (Figure 10) and latent and sensible heat fluxes (Figure 11), suggesting that the modification added in section 2.2.2 (i.e., setting   equal to 100 and applying u * rather than  in equation 7), although being more realistic, do not influence the overall performance of WRF-rLake significantly.BL performs better in terms of RMSE of simulated LSTs, which is 1.14 °C compared to 1.34 °C by Rou_1 and 2.46 °C by Rou_2.When the roughness lengths are fixed to 1 mm (Rou_1), an increase in sensible heat fluxes up to 30 W m -2 is caused in cold seasons but slight decreases are resulted in warm seasons (Figure 11a).A larger effect is observed in latent heat flux, which is increased throughout the year by up to 30 W m -2 (Figure 11b).Annual average LST is reduced by ~1 °C due to excessive outgoing heat fluxes, mainly the latent heat flux.These effects are amplified when fixing the roughness lengths at 10 10 mm (Rou_2), resulting in up to 100 W m -2 increases in winter and 50 W m -2 decreases in summer for sensible heat fluxes and up to 200 W m -2 increases throughout the year for latent heat fluxes.With Rou_2, annual average LST is accordingly reduced by ~2.5 °C compared to BL.

Effects of light extinction coefficient 5
Light extinction coefficient is the key parameter controlling the absorption and distribution of solar radiation in the lake body.BL, Ext_1, Ext_2, and Ext_3 form a group of experiments for extinction coefficient by prescribing it as 1.00 m -1 , 0.13 m -1 , 0.30 m -1 , and 3.00 m -1 respectively.The control of extinction coefficient over energy distribution is reflected in the monthly averaged vertical temperature profiles (Figure 12).In general, as the extinction coefficient increases, a larger portion of energy from solar radiation is maintained in the shallow layers, producing shallower stratification.As the extinction coefficient decreases, more solar radiation penetrates to depth, resulting in a better-developed epilimnion (the top-most and well mixed layer in a thermally stratified lake, occurring above the deeper hypolimnion).BL and Ext_3 produced very similar temperature profiles, which indicates that when the extinction coefficient is sufficiently large, increasing it merely brings in marginal influences in the vertical temperature profile.Specifically, for the Nuozhadu Reservoir, the threshold is ~1.0 m -1 .Among the four experiments, Ext_1 applied the default light extinction coefficient by WRF-Lake (the smallest).This allowed solar radiation to penetrated deeper and energy to be distributed more evenly in the lake body, resulting in less temperature stratification in the topmost 10 m almost all year around.The cumulative influence of lower lake opacity in Ext_1 is also well manifested in seasonal LST.In spring and summer, Ext_1 simulated lower LSTs than other scenarios as more solar radiation penetrated into and warmed up deeper layers.Further, with lower LSTs, less sensible and latent heat were dissipated and more energy was maintained within the lake body.By the end of the year, Ext_1 resulted in an integrally warmer water body than other scenarios with higher water temperatures of 1.5 -2 °C.

Uncertainties and limitations
Ice and snow processes could play a significant role in lake-atmosphere interactions (Brown and Duguay, 2010), especially for high-latitude lakes (e.g., North Eurasian lakes (Subin et al., 2012), Great Lakes (Xiao et al., 2016)).However, given the warm climatology of the Nuozhadu Reservoir, we only examined here the performance of WRF-rLake under icefree conditions.Future work should be carried out to assess WRF-rLake performance at more reservoirs or lakes with icecovered periods as well as different bathymetry and climate to evaluate the broader model applicability.
Due to the lack of Secchi disk measurements at the Nuozhadu Reservoir, we were unable to further improve the radiation scheme in WRF-rLake.However, we note that lake biology is a dominant factor influencing lake optical properties (Cristofor et al., 1994), which itself is affected by many processes, including lake hydrology and biogeochemistry.As modelling such processes is beyond the scope of the current WRF-rLake module, the parameterization of light extinction in WRF-Lake is simply based on lake depth (equation 3).Considering the remarkable impacts of lake opacity on lake energy distribution and mixing regime, future developments should include more sophisticated parameterizations of extinction coefficient with considerations of lake hydrology, chemistry, and biology.
Operation-induced inflows and outflows are key features of artificial reservoirs and can strongly affect seasonal and interannual evolution of reservoir surface water levels, intensity of thermal stratification, and thermal structure (Anohin et al., 2006;Çalışkan and Şebnem, 2009).Given that reservoirs are essential infrastructures for utilisation and management of water resources (Jain and Singh, 2003;Ahmad et al, 2014), the WRF-rLake framework should be extended to include reservoir operation features (e.g., inflow and outflow controls) to better characterize reservoir-atmosphere interactions.

Conclusions
In this study, we revised the WRF-Lake model by adding a new spatial discretization scheme, modifying surface property and vertical diffusivity parameterizations, and adopting a revised convection scheme.The revised lake model, WRF-rLake, was evaluated at the deep Nuozhadu Reservoir in southwestern China and demonstrated overall improved performance in simulating water temperatures in comparison with WRF-Lake, the current lake model in WRF.Compared to WRF-Lake, WRF-rLake reduced the RMSE against observed lake surface temperatures (LSTs) from 1.35 °C to 1.14 °C and that against monthly observed lake temperatures profiles from 1.51 °C to 1.13 °C.
Based on the comparison between WRF-Lake and WRF-rLake, we found that the coarse discretization of water layers in the current WRF-lake made it less able to accurately predict temperatures in the thermocline.An adaptive 25-layer scheme was thus introduced to WRF-rLake and demonstrated better performance than the original 10-layer scheme by reducing the RMSE against observed lake temperatures profiles from 1.64°C to 1.13 °C.
In WRF-Lake, the lake surface roughness length is prescribed to be 1 mm, which we found led to biased simulated surface fluxes and surface temperatures.Lakes have much smoother surfaces, and thus smaller roughness lengths, compared to land and oceans, and lake roughness lengths vary with wind and surface waves.Replacing the original parameterization of roughness lengths with our proposed parameterization (equations 5, 6, 7, 8) reduced latent heat flux by up to 30 W m -2 , and considerably improved LST simulations by reducing the RMSE against observations from 1.34 °C to 1.14 °C.
Simulations of temperature stratification and surface fluxes are sensitive to lake opacity, which modulates the absorption of solar radiation in the lake body.Considering that lake opacity may vary over more than two orders of magnitude (Subin et al., 2012), more detailed global datasets on lake opacity based on remote sensing should be developed.Field measurements of extinction coefficients will be critical for achieving high quality weather and climate simulations in lake-rich areas.
Previous studies have recognized that the wind-driven eddy diffusivity parameterization by Hostetler and Bartlein (1990) is insufficient to simulate large and deep lake surface fluxes and water temperatures.Gu et al. (2015) therefore increased eddy diffusivity for deep lakes with large multiplicative factors that depend on lake depth and surface temperature.However, we found such treatment resulted in unrealistically large eddy diffusivity in about the top 10 m of the Nuozhadu Reservoir and failed to compensate for the unresolved 3D mixing processes in deeper layers.WRF-rLake thus adopts a maximum of 10 2 cm 2 s -1 for eddy diffusivity to avoid overestimation in the surface layers and an enhanced diffusion term for deep lakes, resulting in more realistic eddy diffusivities in deeper layers.In the case of the Nuozhadu Reservoir, adopting all the modifications to eddy diffusivity produced overall similar diffusivity to those measured in Lake Zürich (a lake with similar topography and depth).Although considerable improvements have been brought in by WRF-rLake, the parameterization for vertical diffusivity should be further evaluated and improved to resolve remaining uncertainties in its performance at other lakes, especially for deep lakes.We note several other limitations that could guide future work, e.g., the need to evaluate ice and snow processes and to improve the radiation scheme and the in/outflow parameterization.Our future work will couple the WRF-rLake module with the WRF framework to examine the performance of the coupled system.Overall, we find that simulations of lake temperatures and surface energy balances were improved by WRF-rLake by modifications to the discretization scheme, lake opacity, parameterization for surface properties and vertical diffusivity, and the convection scheme.

Figure 1 .
Figure 1.Flow chart of the WRF-rLake model.The yellow star indicates the process is modified or newly added.

Figure 3 .
Figure 3.The year-long meteorological forcing data used include (a) shortwave and longwave radiation, (b) air temperature, (c) wind speed, (d) precipitation and (e) humidity.All data are averaged to produce mean daily values.

Figure 4 .
Figure 4. Measured monthly water temperature profile for Nuozhadu Reservoir in year 2015.The light gray line indicates the water level variation throughout the year.

Figure 5 .
Figure 5.Time series of lake surface temperatures by observation (red dots), baseline (BL: black line), control (CTL: blue line), and Diff_3 (green line), and air temperature (grey dashed line) of Nuozhadu Reservoir for the period 1 January 2015 through 31 December 2015.The configuration Diff_3 (with "half" modified diffusivity) overestimated LSTs by up to 5.55 °C, reaching an MBE of 0.61 °C.Seasonally, the LSTs are very well reproduced in the 1 st and 4 th quarter but are overestimated in the 2 nd and 3 rd quarter of the year.The vertical temperature profile (Figure6) shows that in the 2 nd and 3 rd quarters, Diff_3 simulated too much vertical mixing in the top 10 meters, resulting in warmer water temperatures in this zone.Meanwhile, a sharp temperature decline is observed between 10 and 20 m depth followed by an underestimation of temperature below 20 m, suggesting too

Figure 6 .
Figure 6.Monthly vertical temperature profile for the first 60 m water in year 2015 by observation (red dots), BL (black line), CTL (blue line) and Diff_3 (green line).

Figure 7 .
Figure 7. Monthly vertical temperature profile for the first 60 m water by observation (red dots), Lyr_1 (red line) and Lyr_2 (black line) in year 2015.4.3Effects of diffusivity BL, Diff_1, Diff_2, and Diff_3 form a group of sensitivity experiments for diffusivity.In the case of the Nuozhadu 5 Reservoir, Diff_2 applies the eddy diffusivity of Diff_1 with an increase of 100 times throughout all layers as is in Gu et al. (2015); Diff_3 additionally includes the enhanced diffusion term on top of Diff_2.

Figure 8 .
Figure 8. Monthly vertical temperature profiles for the first 60 m water by observation (red dots), BL (black line), Diff_1 (red line), Diff_2 (blue line), and Diff_3 (green line) in year 2015.Note Diff_2 and Diff_3 overlaps each other.

Figure 9 .
Figure 9. Monthly vertical diffusivity profile for the first 60 m water by BL (black line), Diff_1 (red line), Diff_2 (blue line), and Diff_3 (green line) in year 2015.The gray shading indicates the diffusivity range of Lake Zürich reported byLi (1973).

Table 2 .
Fixed factors (FF) for vertical discretization for different lakes based on depth.

Table 3 .
An overview of numerical experiments designed to demonstrate sensitivity in WRF-rLake.
*: configuration of Lyr_2 is the same as BL.

Table 4 .
Statistics of the discrepancy between simulated (CTL, Diff_3, and BL) and observed LSTs and monthly temperature profiles during 10 year 2015.Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) are calculated between each simulation and measurement.Mean Bias Error (MBE), Max Bias, and Min Bias are computed by simulation minus measurement.Coral and green shadings indicate the largest and smallest absolute values among three simulations, respectively.