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, 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 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 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 modeled lake processes. We found (1) large changes in surface energy balance fluxes (up to 60 W m⁻²) 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 and reservoirs.

1 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 influence on local and regional climate and are important in understanding lake–atmosphere 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., 2012b; Deng et al., 2013; Xiao et al., 2016).

Lakes interact with atmosphere through energy, mass (mostly water), and momentum exchanges (Lerman and Chou, 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., 2012a). 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, 2013; Subin et al., 2012a, 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 (1-D) models have been widely used for coupling with atmospheric models, because 1-D 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 1-D 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, 1991), CLM4-LISSS (Subin et al., 2012a), 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 k–ε 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 (k) 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 two-layer model based on the similarity theory; Mironov et al., 2010) as a typical example. Gula and Peltier (2012) and Mallard et al. (2014) have coupled FLake to WRF in one-way and two-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-modeling 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, 2016; Mallard et al., 2015; 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 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.

Table 1 Common 1-D lake models.

2 Model description and improvements

2.1 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 1-D 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., 2012a). 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.

2.1.1 Surface energy balance

In the WRF-Lake model, the energy balance at the lake surface is as follows:

$$\begin{array}{}\text{(1)}& S+L-G=H+\mathit{\lambda }E,\end{array}$$

where S is absorbed shortwave radiation, L is net longwave radiation, G is downward heat flux into the lake, H is sensible heat, and λE is latent heat. The units for all variables in Eq. (1) are watts per square meter (W m⁻²). In this equation, S is regulated by albedo (a); G, H, and λE are largely dependent on the thermal conductivity of the top lake layer (κ), aerodynamic resistance for heat (r_ah), and that for vapor (r_aw), respectively. A more thorough discussion of the relationship between lake surface fluxes and aerodynamic resistances is provided by Sect. 2.1.8 in Subin et al. (2012a).

2.1.2 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:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\varphi }=\left(\mathrm{1}-a\right)(\mathrm{1}-\mathit{\beta })S_i{e}^{-\left[\mathit{\eta }\left(z-z_{\mathrm{a}}\right)\right]},\end{array}$$

where ϕ (W m⁻²) is the distribution of solar radiation with depth z (m; positive downward), β is set to a constant 0.4 (Oleson et al., 2010), z_a is the base of the surface absorption layer (0.6 m), and η is the light extinction coefficient (m⁻¹) which is calculated as a function of the lake depth according to Håkanson (1995):

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\eta }=\mathrm{1.1925}{d}^{-\mathrm{0.424}},\end{array}$$

where d (m) is the lake depth.

Though there exist more sophisticated radiation schemes in other lake models (e.g., the nine-band scheme by Paulson and Simpson, 1981), we kept the current WRF-Lake radiation scheme since it includes the essential components in a water-body 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 1-D lake models (Fang and Stefan, 1996; Stepanenko and Lykosov, 2005). To improve the model performance, we tentatively set the cutoff depth z_a 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 z_a value (0 m) demonstrates acceptable performance in this work, a more lake-specific cutoff depth may be needed for better model performance.

2.1.3 The heat diffusion solution

After the surface energy balance and radiation absorption are calculated, the following 1-D thermal diffusion equation is solved:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(k{\displaystyle \frac{\partial T}{\partial z}}\right)-{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{w}}}}{\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}z}},\end{array}$$

where c_w is the volumetric heat capacity of water (J K⁻¹ m⁻³), T is water temperature (K), t is time (s), z is depth (m), k is the thermal diffusivity (m² s⁻¹), and ϕ is the solar radiation as in Eq. (3).

The thermal diffusivity k in WRF-Lake comprises molecular diffusivity (k_m) and wind-driven eddy diffusivity (k_e). k_m is a function of thermal conductivity and volumetric heat capacity of water ($\mathrm{1.433}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹). k_e follows the Henderson-Sellers parameterization (shown in Appendix A) and is determined by wind speed at 2 m above the water surface, a latitude-dependent Ekman profile, and a 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. (2012a) 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² to 10⁵ depending on lake depth and surface temperature (Gu et al., 2015).

2.1.4 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).

2.2 Improvements by WRF-rLake

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

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

2.2.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 nine 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 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 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, and 44.7 m. Similar 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.

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

2.2.2 Surface properties and parameters

As discussed in Sect. 2.1.1, the aerodynamic resistances for heat (r_ah) and vapor (r_aw) heat fluxes are critical for surface energy balance predictions. The aerodynamic resistances are functions of momentum (z_0 m) and scalar roughness lengths (z_0 h for sensible heat and z_0 q for latent heat). In WRF-Lake, z_0 m is set to 1, 5, and 2.5 mm for unfrozen lakes, frozen lakes without snow, and frozen lakes with snow, respectively; while z_0 h and z_0 q are always kept equal to z_0 m.

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 z_0 m=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 z_0 m value. For example, measurements on Lake Washington in the US show z_0 m is generally below 1 mm and ranged between 0.01 and 1 mm (Atakturk and Katsaros, 1999); measurements on Lake Ngoring, a high-altitude lake in the Tibetan Plateau, found z_0 m ranged between 0.001 and 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 z_0 m based on the forcing wind, friction velocity, fetch, and lake depth. Adopting these parameterizations has produced more accurate surface heat fluxes and lake surface temperatures (LSTs) over many natural lakes (Subin et al., 2012a).

We therefore adopted the CLM4-LISSS parameterization of roughness lengths with some further modifications:

• z_0 m is set to 2.4 mm for frozen lakes with snow and 1 mm for frozen lakes without explicit snow layers, and z_0 h and z_0 q are computed as functions of z_0 m and friction velocity u_* (m s⁻¹).

• For unfrozen lakes, z_0 m is parameterized as follows:

  $$\begin{array}{}\text{(5)}& {\displaystyle }{z}_{\mathrm{0}m}=max\left({\displaystyle \frac{\mathit{\gamma }\mathit{\nu }}{u_{*}}},\mathit{\alpha }{\displaystyle \frac{u_{*}}{g}}\right),\end{array}$$

  where γ is a dimensionless empirical constant (0.1), α is the dimensionless Charnock coefficient described below, and g is the acceleration of gravity (Fairall et al., 1996; Charnock, 1955; Smith, 1988):

  $$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }={\mathit{\alpha }}_{min}+\left({\mathit{\alpha }}_{max}-{\mathit{\alpha }}_{min}\right)\exp \left[-min\left(A,B\right)\right],\text{(7)}& {\displaystyle }& {\displaystyle }A={\left({\displaystyle \frac{Fg}{u^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{3}}/f_{\mathrm{c}},\text{(8)}& {\displaystyle }& {\displaystyle }B=\mathit{\epsilon }{\displaystyle \frac{\sqrt{dg}}{u}},\end{array}$$

  where F (m) is the lake fetch (assumed to be 25 times the lake depth when observations are unavailable), u (m s⁻¹) is 2 m wind speed, d (m) is lake depth, ${\mathit{\alpha }}_{min}=\mathrm{0.01}$, ${\mathit{\alpha }}_{max}=\mathrm{0.11}$, and ε=1. A and B account for influences from fetch and depth, respectively.

We also made further improvements in roughness length parameterizations for unfrozen lakes. In LISSS, f_c should be 100 but was tentatively set to 22, corresponding to the use of u instead of u_* in Eq. (7). Subin et al. (2012a) recommended that future lake models should relax this assumption by setting f_c=100 and directly applying u_* rather than u, which we have done in WRF-rLake. Because u_* depends on surface roughness lengths, which in turn depend on u_*, we introduced a fixed-point 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.

2.2.3 Improved diffusivity parameterization

The Henderson-Sellers parameterization (see Appendix A) 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 k_e declines exponentially with depth (Eq. A1), it is more likely to be underestimated in deeper layers than in the topmost one. Thus, enlarging k_e by the same factor for the whole lake may introduce new problems in two ways.

First, k_e 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² s⁻¹ and in oceans, where the forcing winds are usually stronger and moving surface waves are better developed, vertical heat diffusivity exceeds ∼10² cm² s⁻¹ (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² cm² s⁻¹ for wind-driven eddy diffusivity (corresponding to the maximum value for open seas) to avoid overestimation in surface layers.

Second, for deep layers, k_e may still be underestimated by WRF-Lake because k_e 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, D_ed (m² s⁻¹), from CLM4-LISSS (Ellis et al., 1991; Fang and Stefan, 1996; Subin et al., 2012a), which was originally suggested by Hondzo and Stefan (1993) to compensate for unresolved turbulence (even below ice or at large depth) and is given as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{D}_{\mathrm{ed}}=\mathit{\delta }{\left(N^{\mathrm{2}}\right)}^{-\mathrm{0.43}},\end{array}$$

where δ is the level of turbulence which can be related to lake area, and N (s⁻¹) is the Brunt–Väisälä frequency (a minimum N² is suggested to be $\mathrm{7.5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻²). CLM4-LISSS adopted the values measured at an ice-covered lake (Fang and Stefan, 1996) for Eq. (9) by setting $\mathit{\delta }=\mathrm{1.04}\times {\mathrm{10}}^{-\mathrm{8}}$. However, as discussed by Subin et al. (2012a), this suggested value for D_ed is only several times larger than k_m 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 but admit that this coefficient may vary from lake to lake (partially explained by lake area) and thus needs to be tuned under specific scenarios.

2.2.4 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⁻³ K m⁻¹, 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⁻⁴ kg m⁻³ m⁻¹ (which is equivalent to a temperature gradient threshold of about 10⁻³ K m⁻¹) to avoid this unphysical convective mixing. In our tests at Nuozhadu Reservoir, convective mixing can penetrate as deep as 5 m below the surface in WRF-Lake, while it is almost always restricted within the top 1 m in WRF-rLake.

3 Data and modeling experiments

3.1 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 (Fig. 2). The Nuozhadu Reservoir has a dam height of 262 m and a normal water level¹ of 812 m a.s.l. The water depth upstream of the dam in the 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.

Figure 2 Location of the Nuozhadu Reservoir.

3.2 Forcing data

Our study period covers from 1 January to 31 December 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 were repeated 7 times to form a 1-week spinup. The downward shortwave radiation and downward longwave radiation were obtained from the China Meteorological Forcing Dataset (Kun et al., 2010; Chen et al., 2011), which has a temporal resolution of 3 h and spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$. 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 considered to be an acceptable approximation given no data of higher temporal resolution are 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 h temporal resolution from a meteorological observatory run by China Huaneng Group Co., Ltd., the construction unit of the Nuozhadu Reservoir (Fig. 3). The station (22^∘40^′ N, 100^∘23^′ E) is located about 5 km upstream (or northwest) of the Nuozhadu Dam and has an observational height of 10 m.

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.

3.3 Observations for model initialization and evaluation

Water temperature was measured hourly from 712 m to 804 m a.s.l. at an interval of 2 m near the Nuozhadu dam during 2015 (Fig. 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 temperature on 1 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 (Shen, 2017), a three-dimensional hydrodynamic model proven to be sufficiently accurate at the Nuozhadu Reservoir.

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

3.4 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 terms, and the scheme as discussed in Sect. 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, 10 mm, calculated as in Subin et al. (2012a), or calculated with further modification as in Sect. 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 modeling deep lakes (Hocking and Straškraba, 2015). 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⁻¹ (default), 0.30, 1.00, and 3.00 m⁻¹. Although measurements have reported larger variability in the light extinction coefficient (e.g., 0.05 to 7.1 m⁻¹ in Subin et al., 2012a), we found simulated temperature profiles were insensitive to values outside of the 0.13 to 3.0 m⁻¹ range. We concluded that the best performance could be achieved by increasing the light extinction coefficient to ∼1.00 m⁻¹, 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 (i.e., our proposed new model structure) was configured with all modifications described in Sect. 2.2 and a tuned light extinction coefficient to demonstrate the effects of each improvement in WRF-rLake.

Table 3 An overview of numerical experiments designed to demonstrate sensitivity in WRF-rLake.

^a Configuration of Lyr_2 is the same as BL.
^b As the reservoir was ice-free in the year 2015, the constants given in this column refer to the roughness lengths for unfrozen lakes.

4 Results and discussion

4.1 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 (Fig. 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 first 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 overly large outgoing latent and sensible heat fluxes (Sect. 4.2).

Figure 5 Time series of lake surface temperatures by observation (red dots), baseline (BL: black line), control (CTL: blue line), Diff_3 (green line), and air temperature (gray 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 first and fourth quarter but are overestimated in the second and third quarter of the year. The vertical temperature profile (Fig. 6) shows that in the second and third quarters, Diff_3 simulated too much vertical mixing in the top 10 m, 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 overly weak vertical mixing simulated below 20 m. Further discussion in terms of diffusivity is shown in Sect. 4.3. Here the results of other diffusivity experiments (i.e., Diff_1 and Diff_2) are not shown.

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

The simulation by BL is better in terms of root-mean-square error (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 m 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.

Table 4 Statistics of the discrepancy between simulated (CTL, Diff_3, and BL) and observed LSTs and monthly temperature profiles during 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. Bold and italic font indicate the largest and smallest absolute values among three simulations, respectively.

Applying all the modifications to vertical diffusivity (discussed in Sect. 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.

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 the year 2015.

4.2 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 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 (Fig. 7). In the first 2 months, the temperature difference in the 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 and 20 W m⁻², 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 100-layer 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.

4.3 Effects of diffusivity

BL, Diff_1, Diff_2, and Diff_3 form a group of sensitivity experiments for diffusivity. In the case of the Nuozhadu Reservoir, Diff_2 applies the eddy diffusivity of Diff_1 with an increase of 100 times throughout all layers, as in Gu et al. (2015); Diff_3 additionally includes the enhanced diffusion term on top of Diff_2.

For the monthly-averaged vertical temperature profile (Fig. 8), the BL scheme best captures the decrease in 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, 1.47, and 1.47 ^∘C, respectively. Diff_1 yields strong stratification within 10 m 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 LSTs that are biased high in summer and LSTs that are biased low in winter.

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 the year 2015. Note Diff_2 and Diff_3 overlap each other.

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 uniform temperatures in the top 10 m (Fig. 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 (Fig. 9). The overall diffusivity by Diff_2 in the top 10 m could be as large as 10³ cm² s⁻¹, which is too large compared to estimates in the literature (Sect. 2.2.3), supporting the necessity for limiting eddy diffusivity (Sect. 2.2.4).

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 the year 2015. The gray shading indicates the diffusivity range of Lake Zürich reported by Li (1973).

The overall vertical diffusivity below 20 m by Diff_2 or Diff_3 is of the same magnitude as molecular diffusivity ($\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{3}}$ cm² s⁻¹). By increasing the enhanced diffusion term by 100 times for deep lakes, BL yielded more reasonable vertical diffusivity below 20 m (Fig. 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 and 1 cm² s⁻¹.

With the constrained eddy diffusivity and “enhanced diffusion term”, BL produced the best temperature profiles compared to observations (Fig. 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.

4.4 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, 10 mm, and the parameterization from the Subin et al. (2012a) 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 (Fig. 10) and latent and sensible heat fluxes (Fig. 11), suggesting that the modification added in Sect. 2.2.2 (i.e., setting f_c equal to 100 and applying u_* rather than u in Eq. 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.

Figure 10 Time series of LST by observation (red dots), BL (black line), Rou_1 (red line), Rou_2 (blue line), Rou_3 (green line), and air temperature (gray dashed line) of Nuozhadu Reservoir in the year 2015. Note that BL and Rou_3 overlap most of the time.

When the roughness lengths are fixed to 1 mm (Rou_1), an increase in sensible heat fluxes up to 30 W m⁻² is caused in cold seasons but slight decreases are seen in warm seasons (Fig. 11a). A larger effect is observed in latent heat flux, which is increased throughout the year by up to 30 W m⁻² (Fig. 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 mm (Rou_2), resulting in up to 100 W m⁻² increases in winter, 50 W m⁻² decreases in summer for sensible heat fluxes, and up to 200 W m⁻² increases throughout the year for latent heat fluxes. With Rou_2, annual average LST is accordingly reduced by ∼2.5 ^∘C compared to BL.

Figure 11 Time series of (a) upward sensible heat flux, (b) upward latent heat flux by BL (black line), Rou_1 (red line), Rou_2 (blue line), and Rou_3 (green line) of Nuozhadu Reservoir in 2015.

4.5 Effects of light extinction coefficient

The 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 the extinction coefficient by prescribing it as 1.00, 0.13, 0.30, and 3.00 m⁻¹ respectively. The control of the extinction coefficient over energy distribution is reflected in the monthly averaged vertical temperature profiles (Fig. 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⁻¹.

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

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 year-round. 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.

4.6 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., 2012a; 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 ice-free conditions. Future work should be carried out to assess WRF-rLake performance at more reservoirs or lakes with ice-covered 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 modeling 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 (Eq. 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 Elçi, 2009). Given that reservoirs are essential infrastructures for the utilization 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.

5 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 to 1.14 ^∘C and that against monthly observed lake temperatures profiles from 1.51 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. A 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 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 (Eqs. 5, 6, 7, 8) reduced latent heat flux by up to 30 W m⁻², and considerably improved LST simulations by reducing the RMSE against observations from 1.34 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 by more than 2 orders of magnitude (Subin et al., 2012a), 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 3-D mixing processes in deeper layers. WRF-rLake thus adopts a maximum of 10² cm² s⁻¹ 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 inflow–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 convection scheme.

Appendix A: Parameterization for wind-driven eddy diffusivity

In the Henderson-Sellers parameterization scheme, wind-driven eddy diffusivity (k_e) is determined by wind speed at 2 m above the water surface, a latitude-dependent Ekman profile, and a lake-stratification-dependent Brunt–Väisälä frequency (Henderson-Sellers et al., 1983; Henderson-Sellers, 1985):

$$\begin{array}{}\text{(A1)}& {\displaystyle }{k}_{\mathrm{e}}=\left(k{u}_{*}z/P_{\mathrm{0}}\right)e^{\left(-k^{*}z\right)}{\left(\mathrm{1}+\mathrm{37}{\mathit{\text{Ri}}}^{\mathrm{2}}\right)}^{-\mathrm{1}},\end{array}$$

where k(=0.4) is von Kármán's constant, u_* is the surface friction velocity (m s⁻²), P₀(=1.0) is the neutral value of the turbulent Prandtl number, k^* is a latitudinally dependent parameter of the Ekman profile, and Ri is the gradient Richardson number.

The parameter u_* is given by the following:

$$\begin{array}{}\text{(A2)}& {\displaystyle }{u}_{*}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{3}}u,\end{array}$$

where u (m s⁻¹) is wind speed at 2 m a.g.l.

The parameter k^* is expressed as follows:

$$\begin{array}{}\text{(A3)}& {\displaystyle }{k}^{*}=\mathrm{6.6}\times {\left(\sin \mathit{\varphi }\right)}^{\mathrm{1}/\mathrm{2}}{u}^{-\mathrm{1.84}},\end{array}$$

where ϕ is the latitude of the lake being modeled.

The parameter Ri is determined as follows:

$$\begin{array}{}\text{(A4)}& {\displaystyle }\mathit{\text{Ri}}={\displaystyle \frac{-\mathrm{1}+{\left\{\mathrm{1}+\mathrm{40}{N}^{\mathrm{2}}{k}^{\mathrm{2}}{z}^{\mathrm{2}}/\left[u_{*}^{\mathrm{2}}\exp \left(-\mathrm{2}{k}^{*}z\right)\right]\right\}}^{\mathrm{1}/\mathrm{2}}}{\mathrm{20}}},\end{array}$$

where N (s⁻¹) is the Brunt–Väisälä frequency specified as follows:

$$\begin{array}{}\text{(A5)}& {\displaystyle }N={\left[-g/\mathit{\rho }\left(\partial \mathit{\rho }/\partial z\right)\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$