Improvements of the coupled WRF-Lake model over Lake Nam Co, Central Tibetan Plateau

A series of model sensitivity simulations are carried out to calibrate and improve the Weather Research and Forecasting Model coupled with a one-dimensional lake model (WRF-Lake) based on observations over Lake Nam Co. Using the default lake model parameters, the solution of WRF-Lake exhibits significant biases in both the lake thermodynamics and regional climatology, i.e., higher lake surface temperature (LST), earlier onset of summer thermal stratification, and overestimated near-surface air temperature and precipitation induced by the lake’s excessive warming and moistening impacts. The performance of WRF-Lake is improved through adjusting the initial lake temperature profile, the temperature of maximum water density (Tdmax), the surface roughness length, and the light extinction coefficient. Results show that initializing the water temperature with spring observation mitigates the LST overestimation and reduces the timing error of the onset of thermal stratification. By further adjusting Tdmax from 4 °C to the observed value of 3.5 °C, the LST increase from June to mid-July is enhanced and the buildup of thermal stratification is more accurately predicted. Through incorporating the parameterized surface roughness length and decreasing the light extinction coefficient, the model better reproduces the observed daily evolution of LST and vertical lake temperature profile. The calibrated WRF-Lake effectively mitigates the overestimation of over-lake air temperature at 2 m height and precipitation over regions downwind the lake. This suggests that an improved lake scheme within the coupled WRF-Lake is essential for realistically simulating the lake–air interactions and the regional climate over the lake-rich Tibetan Plateau.


Introduction
The Tibetan Plateau (hereafter TP), known as the "Asian Water Tower", harbors more than 1400 lakes with an individual surface area above 1 km 2 (Ma et al. 2011;Zhang et al. 2019a, b). The total lake surface area over TP is approximately 47,000 km 2 and accounts for more than 52% of the total lake coverage in China (Jiang and Huang 2004;Song et al. 2013). As the water body is characterized by lower albedo, smaller surface roughness length, larger heat capacity and thermal conductance when compared with the other land types, the exchanges of water and energy at the lake-air interface distinctly differ from those at the land-air interface, causing the spatial heterogeneity of surface thermal conditions over TP (Gerken et al. 2013;Ma et al. 2014;Song et al. 2016). This would further influence the atmospheric circulation, the weather and climate at local and regional scales Zhu et al. 2017). In response to the ongoing climate change and anthropogenic forcing, the TP lakes experience variations in terms of water level/areas 1 3 (Jiang et al. 2017;Kleinherebrink et al. 2015;Zhang et al. 2017), aquatic biogeochemistry (Karl et al. 2018;Sahoo et al. 2016), and ice coverage Qi et al. 2019). The TP lakes also take significant feedbacks to the climate and environmental changes (Guo et al. 2019;Neckel et al. 2014;Liao et al. 2013;Yan et al. 2018). Hence, it is imperative to have a comprehensive understanding of the lake-air interactions over TP.
Owing to the high terrain and harsh climate conditions over TP, the in-situ lake observations are limited and inconsecutive, far from being adequate for studying the alpine lake-air interactions and the associated climatic effects Wang et al. 2009Wang et al. , 2015Wang et al. , 2019Biermann et al. 2014). As an alternative, the high-resolution weather/ climate models with the consideration of lake impacts are able to simulate the lake-induced spatiotemporal heterogeneity of meteorological elements (i.e., precipitation, temperature and turbulent fluxes), to elucidate the involved underlying mechanisms, and to project future lake variations under different climate change scenarios (Hwang et al. 2019;Wang et al. 2017;Wen et al. 2015;Xiao et al. 2018;Zhu et al. 2017). During the recent decades, many one-dimensional (1-D) and three-dimensional (3-D) lake models with different complexities have been developed and improved to study the lake-air interactions. In contrast to the 1-D models, the 3-D models can resolve both the lake hydrodynamics and thermodynamics explicitly, hence can accurately depicting the lake physical characteristics Notaro et al. 2013;Xue et al. 2016). However, the two-way coupling of a 3-D lake model with the other components of a regional climate model requires much more computational power due to the high horizontal resolution (around 2 km) and the complicated physical processes (Nyamweya et al. 2016;Pan et al. 2002;Song et al. 2004). Hence, at present, 1-D lake models are commonly used in the coupled regional weather prediction and climate models (Lofgren 1997;Martynov et al. 2012;Mironov et al. 2010;Samuelsson et al. 2010;Williams et al. 2015;Theeuwes et al. 2010;Thiery et al. 2015).
There are many 1-D lake models for lake simulations, including the Mixed-Layer Model (Goyette et al. 2000), the Freshwater Model based on self-similarity theory (Mironov 2008), the eddy-diffusive Hostetler model (Hostetler et al. 1993;Subin et al. 2012), the bulk formulation models (Imberger and Patterson 1981), and the finite-difference models with second-order (k-ε) turbulence closures (Perroud et al. 2009). According to the Lake Model Intercomparison Project (LakeMIP, Stepanenko et al. 2010Stepanenko et al. , 2013, these models involve different degrees of simplifications on lake thermo-and hydro-dynamics, and their performances vary with different model configurations. As the lake-air interactions and the lake mixing regimes are significantly influenced by different lake intrinsic features (i.e., depth, size and salinity), local orographic conditions and the background atmospheric circulation, it is essential to identify and calibrate the key model parameters or parameterization schemes for the lake processes within the specific research domains (Gu et al. 2015(Gu et al. , 2016Li et al. 2018;Wen et al. 2016;Thiery et al. 2014;Xiao et al. 2016;Zhang et al. 2019a, b). This is of great benefit to improve the representation of the lake-air interaction in climate models for better simulating the LST and regional climate over lake-dominant areas (Balsamo et al. 2012;Mallard et al. 2014;Martynov et al. 2013;Wang et al. 2018;).
In this study, we focused on the performances of the 1-D eddy diffusion lake scheme incorporated with the advanced Weather Research and Forecasting model (WRF-ARW V3.9.1, Skamarock et al. 2008) in simulating the thermal features of Lake Nam Co, which is the second largest brackish lake over the central TP (Fig. 1c), and the lake climatic effect. The default lake scheme originates from the basic concept of the Hostetler lake model and is further improved by Subin et al. (2012) and Gu et al. (2015). At present, the key parameters (i.e., lake temperature profile initialization, Tdmax, surface roughness length, light extinction coefficient and eddy diffusivity) were calibrated by Gu et al. (2015) based on the buoy observations over two plain freshwater lakes (Lake Erie and Lake Superior). However, as most TP lakes are situated in the alpine cold and semi-arid to arid climate zone and are less influenced by the anthropogenic activities, the lake physical properties and the related lake-air interactions show distinctive features from those of the lakes over lowlands and wet regions (Dai et al. 2018a, b;Lazhu et al. 2016;Li et al. 2015;Wang et al. 2017Wang et al. , 2019. This implies that the lake model with default parameter configurations may not appropriately represent the lake-air heat and water exchanges over TP and would cause large biases in simulating the LST or vertical thermal structures of these TP lakes . Additionally, as the LST and associated turbulent heat fluxes are the important boundary conditions for the atmospheric component, their simulation accuracy directly affects the performances of the lakeresolved climate model and is important for understanding the lake impacts on the regional climate over lake-rich areas (Sharma et al. 2018;Xiao et al. 2018;Xu et al. 2019;Zhang et al. 2016). Hence, based on the observed bathymetry and 1-year in-situ water temperature records of Lake Nam Co, we conduct a series of sensitivity simulations to calibrate the key lake model parameters in describing the lake thermodynamics. The comparisons between the experiments with the default and calibrated lake schemes are used to demonstrate that a better representation of lake-air interactions in the coupled WRF-Lake is indeed beneficial for improving its performance in reproducing the regional climate surrounding Lake Nam Co. Main findings of this work provide valuable information for further improvements and applications Improvements of the coupled WRF-Lake model over Lake Nam Co, Central Tibetan Plateau 1 3 of the coupled lake-air models over alpine and semi-arid to arid lake-rich regions such as TP.
The rest of this paper is organized as follows: The brief descriptions of the coupled WRF-Lake model and experimental designs are firstly introduced in Sect. 2.1. The following Sects. 2.2 and 2.3 present the lake scheme modifications, and the datasets and methods for evaluation. In Sect. 3, we evaluate the performances of the coupled WRF-Lake with the default and modified lake model parameters in simulating the temporal evolution of LST, the vertical water temperature profile, Fig. 1 a Geographic locations and terrain height of the three domains; b the spatial distributions of the lake depth in Lake Nam Co; c the sub-grid lake fraction and d depth of lakes in the nested Domain 3. In a, the blue curves outline the Yangtze and the Yellow River. The black solid and hollow aster symbol in b represent the locations of the lake temperature monitoring station and the automatic weather station at Lake Nam Co, respectively. The color scale of b is same as that of d 2 Model description, lake scheme modifications, evaluation data and methodology

Overview of the WRF-Lake model
This study evaluates the lake-air coupled WRF V3.9.1 model (Skamarock et al. 2008), whose lake component is a 1-D mass and energy balance lake model from the Community Land Model version 3.5 (Gu et al. 2015;Oleson et al. 2004). The lake model contains up to 5 snow layers on the lake ice, 10 or 25 lake water/ice layers, and 10 soil layers for the bottom sediment. We adopt a 25-layer vertical discretization for the lake column because this setup has been proved to possess advantages in simulating the variations of metalimnion temperature, ice coverage and other lake thermodynamic characteristics Xiao et al. 2016). The main physical processes can be divided into two parts: energy exchanges between the lake surface (snow, ice, or water) and the overlying atmosphere and the heat transfer within the lake body ). Following Subin et al. (2012), the surface heat fluxes are solved using the conditions at the atmospheric reference height as a top boundary condition and the top lake model layer temperature (in the lake body for snow depth < 40 mm and in the snow otherwise) as a bottom boundary condition. The lake surface temperature T s is solved simultaneously with the friction velocity, surface roughness, and aerodynamic resistances in order to balance the surface energy budget: where is the surface albedo for the downward solar radiation SW ↓ (W m −2 ), β = 0.4 is the surface absorption rate of the net solar radiation. In the net longwave radiation term LW ↓ − T 4 s , = 0.97 is the lake surface emissivity, LW ↓ is the incident longwave radiation (W m −2 ), σ = 5.67 × 10 −8 Wm −2 K −4 is the Stefan-Boltzmann constant. SH ↑ and LH ↑ are the turbulent fluxes of sensible and latent heat (W m −2 ), and G ↓ is the downward heat flux into the lake surface layer (W m −2 ). For the above variables, positive values mean downward for SW/LW/G, and upward for SH/LH. The subsurface temperature for the lake column including snow, ice, water, and soil is calculated by a diffusion equation: where T is the water temperature (K) at depth z (m); t is time (s); k m =1.43 × 10 -7 m 2 s −1 is the molecular diffusion coefficient; k e is the wind-driven eddy diffusivity for unfrozen lake (m 2 s −1 ); m d is the mixing factor with a range from 10 2 to 10 5 that is adopted to compensate for the unresolved vertical mixing processes in deep lakes; c w is the volumetric heat capacity (Jm −3 K −1 ); and is the net downward solar radiation penetrating to depth z following the Beer-Lambert law, where is the light extinction coefficient and z s is the thickness of surface layer.

Experimental designs
The WRF-Lake model encompasses three one-way nested domains with the horizontal resolutions of 50 km, 10 km, and 2 km from Domain 1 to Domain 3 (Fig. 1a). The model domains are centered at (30.7° N, 90.7° E) with 30 sigma levels in the vertical direction and consist of 100 × 100, 146 × 146, and 261 × 246 grid points (number of grids in the west-east direction multiplies that in south-north direction), respectively. The main physical parameterization schemes are listed in Table 1 and the cumulus convection scheme is turned off in Domain 3. The initial and lateral boundary conditions of atmosphere are generated from the 6-h ERA-Interim reanalysis data with a grid spacing of 0.75° × 0.75° in longitude/latitude (https ://apps.ecmwf .int/ datas ets/data/inter im-full-daily , Simmons et al. 2007). The daily real-time global sea surface temperature dataset with a 0.5° × 0.5° resolution (ftp://polar .ncep.noaa.gov/pub/histo ry/ sst/rtg_low_res/, Thiebaux et al. 2003) provides the oceanic boundary conditions. The Moderate Resolution Imaging Dudhia (Dudhia 1989) Longwave radiation Rapid radiative transfer model (Mlawer et al. 1997) Spectroradiometer (MODIS) land-use categories with a resolution of 10 arc seconds (Friedl et al. 2010) provide the land surface types at each model grid. The global lake database version 2 (GLDBv2) including the lake location and depth with the horizontal resolution of 1 km developed by Kourzeneva (2010), Kourzeneva et al. (2012) and Choulga et al. (2014) (available at https ://www.flake .igb-berli n.de/ ep-data.shtml ) gives the ratio of sub-grid lake area to the grid area (lake area/grid area, hereafter represented by lake fraction) and the lake depth at each domain grid ( Fig. 1c and d). As the depth of Lake Nam Co is missing in the GLDBv2 dataset, it is updated by the observed bathymetry data (Fig. 1b) provided by Wang et al. (2009). The air-lake coupled WRF-Lake with the default lake schemes is applied in the outer two domains, and the temperature, humidity and horizontal wind fields above the planetary boundary layer of Domain 1 are nudged to the ERA-Interim reanalysis. The model outputs from an outer domain are employed to generate the initial and boundary conditions for the next-level inner domain. The one-way nesting method inhibits the feedback from Domain 3 to the outer domains. The detailed parameter adjustments are presented in the following Sect. 2.2 and the related configurations of the sensitivity experiments are listed in Table 2. All the experiments run from May 20th to August 31st in 2013. The first 12 days are for spin-up and the results for summer months (June-July-August, JJA) are analyzed.

Lake scheme modifications
To improve the WRF-Lake model's capability in reproducing the lake thermodynamics and therefore the lake climatic effect over central TP (D03), we evaluated and modified the following lake model parameters based on the different properties between plain fresh lakes and Lake Nam Co. The alpine lakes over central Tibetan share similar characteristics that are different from those of plain freshwater lakes, i.e., the colder and thermally mixed lake temperature in early May, a smaller Tdmax influenced by water salinity and lower surface air pressure, and higher water clarity under the less anthropogenic contaminations. Hence, the following parameter settings determined by the observed lake properties of Lake Nam Co are transferred to the other lake basins in D03.

Lake temperature profile initialization
In the default lake model, the lake temperature of the top layer is initialized with the input surface skin temperature generated by the avg_tsfc module utility program in WRF (Maussion et al. 2013). The initial lake temperature is assigned with the Tdmax (4 °C) for layers below 50 m, and is calculated according to the surface skin temperature, Tdmax (4 °C) and different layer depth for layers above 50 m. This method is based on the shape of the observed lake temperature profile for Lake Superior in September 1st, 2008 (Gu et al. 2015). Based on the in-situ observation at Lake Nam Co provided by Wang et al. (2019), the entire lake body is well mixed and water temperature from top to bottom exhibits small differences at the start time of the simulations (May 20th, 2013). Hence, for the experiments except the control (CTRL) run, we initialize all vertical layers of Lake Nam Co with the value 2.5 °C, which is calculated by averaging the observed temperatures of different lake layers on May 20th, 2013.

Temperature of maximum water density
The temperature of maximum water density (Tdmax) directly corresponds to the heat capacity of hypolimnion and plays an important role in affecting the development of lake thermal stratification. Tdmax is set to 4 °C for plain freshwater lake. This default value may not be applicable to Lake Nam Co because it is an alpine brackish carbonate lake and Tdmax depends on salinity and atmospheric pressure. According to Wang et al. (2019), Tdmax is close to 3.4 °C as indicated by the observed stable deep-water temperature after the formation of thermal stratification, and is close to 3.6 °C based on the empirical calculation using air pressure of 0.57 bars and salinity of 1.7 g/L (Chen and Millero 1986). In our simulations other than CTRL, we set the Tdmax as 3.5 °C.

Lake surface roughness length
The roughness length is a key parameter for calculating the surface turbulent fluxes between the lake and the overlying atmosphere. In the default lake model, the surface roughness lengths for momentum ( z 0m ), sensible heat ( z 0h ) and latent heat ( z 0q ) are set as the same constant value according to different lake states, and for unfrozen lake, z 0m = z 0h = z 0q =0.001 m. The offline WRF-Lake evaluation on Lake Nam Co by Huang et al. (2019) demonstrated that replacing the constant lake surface roughness lengths with the parameterization scheme suggested by Subin et al. (2012) can improve the model performance in reproducing the lake thermodynamics during warm seasons. We test this roughness length parameterization scheme in the sensitivity simulations, represented as where is the kinematic viscosity of air (m 2 s −1 ), u * is the surface friction velocity (m s −1 ), C is the effective Charnock coefficient, g=9.80616 m s −2 is the gravitational acceleration, R 0 = max z 0m u * , 0.1 is the near-surface atmospheric roughness Reynolds number, P r =0.71 is the molecular Prandt number for air, S c =0.66 is the molecular Schmidt number for water in air. More details can be found in Subin et al. (2012).

Light extinction coefficient
The distribution of net downward solar radiation at different lake depth depends on the light extinction coefficient, which is related to the degree of lake clarity. In the default lake model, = 1.1925d −0.424 is expressed as an empirical function of lake depth (d) according to the data from 88 Swedish lakes (Hakanson 1995), implying that the light extinction coefficient experiences horizontal variability with the varying lake depth across the model grids. As the TP lakes are less contaminated by human activities and contain smaller amount of aquatic plants compared with plain lakes, we consider TP lakes more limpid and scale down the default by a factor of 0.8 to allow more solar radiation to penetrate into the deep lake layers. In this case, at the observation site over Lake Nam Co is decreased from 0.175 to 0.14 m −1 , which is closer to the mean estimated value 0.12 m −1 in Wang et al. (2009) and has also been used in Lazhu et al. (2016).

Thermal diffusivity
For mitigating the WRF-Lake's underestimation of deeplake vertical mixing strength and therefore improving the LST simulations over Lake Superior, Gu et al. (2015) increased the wind-driven eddy diffusivity k e by a factor of 10 2 -10 5 based on the different lake depth and surface temperature within the default lake model. For example, k e is multiplied by m d = 100 when lake depth exceeds 15 m and T s is greater than 4 °C. In our simulations, we enlarge k e using m d = 150.

Evaluation data and methodology
The model performances in lake thermal structure are evaluated with the daily lake temperature profile in 2013 measured at a station located at the southeastern part of Lake Nam Co (30°45.74′ N, 90°46.83′ E, denoted by a solid aster symbol in Fig. 1b). The water temperature were measured at depths of 3, 6, 16, 21, 26, 31, 36, 56, 66 and 83 m, with details documented by Lazhu et al. (2016). The simulated wind speed, downward shortwave and longwave radiation over the Lake Nam Co are compared with an hourly forcing data for offline lake modelling ). This dataset was generated by Lazhu et al. (2016) based on the 3-hourly, 0.1° × 0.1° dataset from the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS, available at https ://en.tpeda tabas e.cn/porta l/; Chen et al. 2011;He et al. 2020) and the observation data collected at the meteorological station on the northeastern shore of Lake Nam Co (Fig. 1b). The modelled LST is assessed against the MODIS product (MOD11A1) with a spatial resolution of 1 km (available at https ://modis .gsfc.nasa.gov/data/datap rod/ mod11 .php, Wan et al. 2004), which provides the LST over Lake Nam Co at approximately 11:00 and 21:00 local time. Following Dai et al. (2018a), we only employ the nighttime data for LST evaluation because the nighttime retrieval from satellite remote sensing is more accurate due to the absence of strong solar heating bias (Hook et al. 2003). The model simulated 2-m air temperature (hereafter denoted as T 2m ) is compared against the ITPCAS data, and the precipitation is verified by both the ITPCAS and the Tropical Rainfall Measuring Mission 3B42 datasets (hereafter abbreviated as TRMM3B42, available at ftp://disc2 .nasco m.nasa.gov/ data/ TRMM/Gridded, Huffman et al. 2007). In addition, the hourly, 0.25° × 0.25° ERA5 reanalysis data (available at https ://cds.clima te.coper nicus .eu) is utilized to provide the summer atmospheric circulation information over TP. Specifically, the model simulated LST and water temperature profile for evaluation are taken at the grid point nearest to the Nam Co Station, with the lake depth to be 91.3 m according to the model grid and to be 93 m at the station ). The simulated 25-layer water temperature is linearly interpolated onto the depths of the observational data. Four statistics including the mean bias (BIAS), the root-mean-square error (RMSE), the Pearson temporal/pattern correlation (TC/PC), and the Taylor score (TS) are utilized for the model assessment. The above verification statistics are estimated according to where S i (O i ) denotes the simulated (observed) values with a total of N samples in time or in space; − S ( − O) represents the average of S i (O i ) over N points; σ is the temporal/spatial standard deviation of the simulated values normalized by that of the observed values; and TC 0 = 1 is the achievable maximum temporal/spatial correlation. Higher TC indicates greater similarity in temporal/spatial variation between the simulated and observed values, while lower RMSE indicates that the simulated values are more consistent with the observed values in quantity. As documented in Taylor (2001), TS can quantify how the simulated values match the observed values in terms of both the temporal/spatial pattern and amplitude, and better model performance can be indicated by a large TS closer to the achievable maximum value 1.
Additionally, we utilize the equivalent potential temperature e to describe the atmospheric stability: where T is the air temperature, L v is the latent heat of evaporation (2.5 × 10 6 Jkg −1 ), q v is the specific humidity, C p is the specific heat capacity of air at constant pressure (1005 Jkg −1 K −1 ), P ref is the standard atmospheric pressure (1 × 10 5 Pa), and P is the air pressure.

Control simulation with the default lake scheme
During the simulation period (May 20th to August 31st, 2013), the observed epilimnion temperature almost exhibits an identical value from the lake surface to 6 m depth ). Thus, we employ both the LST and the lake temperature at 3 m depth (hereafter denoted as TLake 3m ) to evaluate WRF-Lake's performance in simulating the lake surface layer temperature. Figure 2a and b compares the CTRL simulated nighttime LST (daily TLake 3m ) versus that in MODIS data (in-situ observations). According to observations, the lake surface layer temperature of Lake Nam Co ranges from 3° to 12 °C and exhibits a continuous warming trend during summer. The WRF-Lake with default lake schemes reasonably reproduces the temporal evolution of near-surface lake temperature with a TC of 0.84 (0.96) in comparison with the MODIS (in-situ) data. However, CTRL overestimates the LST by a maximal BIAS of 6.07 °C (3.27 °C) relative to the MODIS (in-situ) data before mid-July, and underestimates the LST from then on. The RMSE against the MODIS (in-situ) data is 2.15 °C (1.13 °C). The radiative forcing is the primary input energy for the entire lake system and plays an important role in determining the equilibrium temperature at the lake-atmosphere boundary (Jimenez et al. 2014). From Fig. 2c and d, the biases between the modeled and observed downward shortwave radiation (longwave radiation) are 39.3 W m −2 (− 4.4 W m −2 ), suggesting that the atmospheric fluxes generated by WRF are within a plausible range (Xu et al. 2016) and the errors of the simulated LST are mainly resulted from the inaccurate lake physical processes depicted by the default lake scheme in the WRF-Lake model.

Lake surface temperature from simulations with improved lake scheme
To reveal the influences of the calibrated lake model parameters on LST simulations, the summer daily time series of TLake 3m from the in-situ observations and each model simulation are shown in Fig. 3, and the quantitative performances of the simulations are illustrated in Fig. 4. Obviously, the simulated TLake 3m is sensitive to changes in lake model parameters.
From Fig. 3, it is obvious that the overestimated TLake 3m in CTRL during June to mid-July is related to the unrealistically higher initial water temperature profile, which also results in the incorrectly earlier onset of the lake thermal stratification (Fig. 5). By initializing the water column temperature with the observed water temperature (2.5 °C) averaged over different lake water layers, the Sen_LTI experiment with less initial thermal storage obtains TLake 3m in better agreement with observations, as indicated by the slight increase in TC and the significant increase in TS by 0.19 (Fig. 4). Sen_LTI underestimates the summer TLake 3m with a BIAS (RMSE) of − 0.98 °C (1.13 °C).
The intensified solar heating in May-June leads to the springtime warming from the original thermally mixed state in Lake Nam Co. Before the LST reaching Tdmax, the eddy diffusivity is enlarged by 10 4 times (Gu et al. 2015), and thus a large portion of the incoming radiative heat flux is transferred into the deeper lake layers, causing the water column temperature increasing with a rather uniform rate until the hypolimnetic temperature reaches Tdmax. Thus, a decreased Tdmax can reduce the hypolimnetic thermal storage and also fasten LST's ascending speed in response to the radiative forcing at the beginning of springtime warming. From  Fig. 3, by adjusting the Tdmax from 4 °C for the freshwater lake to the observed 3.5 °C in Lake Nam Co, the Sen_Tdmax experiment effectively increases the warming amplitude of TLake 3m from early June to mid-July. Compared with Sen_ LTI, the Sen_Tdmax experiment obtains summer TLake 3m that agrees better with observations, with RMSE (BIAS) of 0.93 °C (− 0.68 °C), while 1.13 °C (− 0.98 °C) for Sen_LTI. and root-mean-square-error (RMSE) between the simulations and observations for LST and TLake 3m are presented in a and b, and the BIAS between modeled and observed SW ↓ and LW ↓ are presented in c and d Fig. 3 The daily TLake 3m (unit: °C) from the station measurements and the simulations at the observation site during June 1st to August 31st 2013 As discussed in many previous studies (Dai et al. 2018c;Xu et al. 2016;Wang et al. 2018), the default surface roughness length for unfrozen lake is often too large, thus inducing unrealistically strong wind-driven mixing right below the lake-atmosphere interface, excessive loss of latent heat flux and the underestimation of LST. In the Sen_SRL experiment, the parameterized roughness length with a dynamical range of 0.1-0.5 mm is used to better represent the lake-air interactions. Compared with Sen_Tdmax, the Sen_SRL experiment obtains a significantly increased TLake 3m (Fig. 3), which is largely due to the decreased latent and sensible heat fluxes. The RMSE (BIAS) of summer TLake 3m in Sen_SRL experiment is 0.74 °C (0.14 °C), reduced from 0.93 °C (− 0.68 °C) in Sen_Tdmax run (Fig. 4). Sen_SRL overcorrects TLake 3m , with a mean positive BIAS of 0.49 °C during July and August.
The light extinction coefficient and eddy diffusivity are two key parameters determining the distribution and transmission of the sub-surface heat fluxes, and therefore the vertical distribution of lake temperature. In general, more solar radiation penetrates into and heats the deep-water layers with the decreasing , resulting in lower LST, a warmer mixolimnion, and a deepened thermocline. With a decreased in the Sen_LEC experiment, the BIAS of TLake 3m in July significantly decreases to 0.09 °C, from 0.54 °C in the Sen_SRL experiment (Figs. 3 and 4). The RMSE (BIAS) of summer TLake 3m in Sen_LEC is 0.71 °C (− 0.07 °C), reduced from 0.74 °C (0.14 °C) in Sen_SRL. For depicting the strong turbulent water mixing in deep lakes, k e was significantly enlarged by a mixing factor ( m d ) in previous studies (Martynov et al. 2012;Subin et al. 2012;Xiao et al. 2016), and larger m d values lead to more heat transfer from the nearsurface water to the deep lake body. By enhancing the eddy diffusivity with m d = 150 in the Sen_MF run, the RMSE (BIAS) of summer TLake 3m is 0.94 °C (− 0.44 °C), larger than 0.71 °C (− 0.07 °C) in the Sen_LEC experiment using m d = 100 . This indicates that m d =100 better represents the vertical mixing strength at the deep layers of Lake Nam Co.

Vertical temperature profile simulations with improved lake scheme
To reveal the impacts of the tuned model parameters on simulating the onset of thermal stratification and the vertical distribution of water temperature, Fig. 5 shows the time-depth distributions of the daily observed and the According to observations (Fig. 5a), the lake water temperature slightly decreases from the lake surface to deep layers until the establishment of the lake thermal stratification on June 12th. As a result of the enhanced solar heating from then on, the upper mixed layer gradually warmed and deepened to ~ 26 m along with the development of thermocline. In general, the evolution features of the vertical temperature profile are reasonably reproduced by all the experiments except for the CTRL run (Fig. 5), which largely overestimates the whole lake body temperature and obtains an incorrectly earlier thermal stratification development due to the higher initial water column temperature. Compared to CTRL, the Sen_LTI experiment with a realistic initial water temperature (2.5 °C) postpones the onset of thermal stratification on June 17th (Fig. 5b and c). Comparison between Sen_LTI and Sen_Tdmax demonstrates that the tuned Tdmax, which made the model more realistic in depicting the stability threshold for lakes to adjust their subsurface thermal structure, results in a more accurate onset of summer thermal stratification on June 11th (Fig. 5c and d). Further model modifications lead to little improvement on the simulated onset of thermal stratification, but influence the water temperature simulations at different lake depth during the stable stratified period.
From Figs. 6 and 7, the lake model's accuracy in reproducing the amplitude and pattern of water temperature variability decreases with the increasing lake depth. At 6 m depth, all experiments perform well in the temporal evolution of water temperature and the pattern and amplitude of the temperature variability with the TC/TS values greater than 0.8, while their capability in quantitatively reproducing the lake temperature varies greatly. By adjusting the initial water column temperature, Tdmax, surface roughness length, and light extinction coefficient one by one from Sen_LTI to Sen_LEC, the WRF-Lake model's ability in simulating the epilimnetic temperature is significantly improved with the RMSE (BIAS) value reduced from 0.97 °C (− 0.84 °C) to 0.59 °C (0.11 °C). From Fig. 6b and c, the observed lake temperature evolutions at 21 m and 56 m depths show a fluctuant ascending trend, which is due to the simultaneous influences of the enhanced convective mixing and the active internal wave motions across the thermocline. However, due to the missing representations of the complex lake interior hydrodynamic processes in the 1-D WRF-Lake, all the simulations fail to depict the variations of metalimnion and hypolimnion temperature indicated by the lower TS values. In terms of the verification statistics including RMSE and BIAS, the Sen_SRL experiment shows relatively better performances in simulating the temperature at 21 m and 56 m depths.

Impacts of the improved lake scheme on atmospheric component
The above evaluation results reveal the improvements in simulating the lake thermodynamics, especially the nearsurface water temperature, through using the calibrated lake scheme in the coupled WRF-Lake model. This may lead to improvement in the simulated water and heat exchanges between the lake surface and the overlying atmosphere.
In this section, the T 2m and precipitation simulated by the CTRL and Sen_LEC runs are compared for demonstrating the effects of the improved LST on regional climate modelling. A sub-region bounded by 30° N-31.5° N and 90° E-91.5° E (i.e. the red rectangle in Fig. 8a) is selected for focusing on the comparisons surrounding Lake Nam Co. From Fig. 8a-c, both the CTRL and Sen_LEC experiments can reasonably capture the general spatial patterns of the summer mean T 2m over central TP, except for some cold (warm) biases in the southern mountainous areas (northern lower lands). As the CTRL run predicts an unrealistically higher near-surface water temperature, the related upward water and heat fluxes are overestimated, leading to positive T 2m simulation biases over the lake and surrounding regions ( Fig. 8b and d). With the improvement in the modelled LST, the summer mean over-lake T 2m in Sen_LEC is overestimated by 1.4 °C, reduced from 1.71 °C in the CTRL run ( Fig. 8d-f). The T 2m simulation over the sub-region is also improved in terms of PC, TS, and RMSE (Table 3).
During the boreal summer, the dominant southwesterly monsoon circulation plays an important role in modulating the TP climate and contributes to the southeast-northwest precipitation gradient over the TP (Fig. 9a) (Maussion et al. 2013). From Fig. 9b and c, both the CTRL and Sen_LEC experiments can generally reproduce the 500 hPa wind fileds over Domain 3, and predict the summer mean precipitation with the TS values of 0.79 and 0.80, respectively. From Fig. 9d and e, both experiments overestimate (underestimate) the rainfall surrounding (over) the Lake Nam Co. The overestimated rainfall around the Nyainqentanglha Mountain and the regions downwind the lake (marked by the blue rectangle in Figs. 9f and 10f) is reduced from CTRL to Sen_LEC, which is attributed to the disparate lake-effect strength. Previous studies have demonstrated that due to the large water vapor supply, the warm surface lake water and the small water surface roughness length in summer, the prevailing southwest wind passing the Lake Nam Co can be accelerated, moistened, and warmed by the underlying water surface (Dai et al. 2018a;Wu et al. 2019). Additionally, under the superimposed impact of mountain-valley wind induced by the sharp orographic contrast, the warm and moist wind would be further amplified and then converge and ascend over the lake and downwind regions, promoting more convective developments therein. By mitigating the LST overestimation in the CTRL run, the Sen_LEC experiment predicts reduced upward moisture and heat fluxes that support the convective stimulation, and thus the decreased lake-effect precipitation strength. The 500 hPa wind filed differences between Sen_LEC and CTRL experiments show divergent flow patterns over the regions downwind the lake, and a reduced positive rainfall bias in the Sen_LEC run (Fig. 8f). Further comparisons between the simulated results and the ITPCAS rainfall dataset also demonstrate that the Sen_LEC experiment can effectively weaken the lake's excessive warming and moistening impacts modeled by CTRL run, and thus mitigate the rainfall overestimation around the Nyainqentanglha Mountain and the regions downwind the lake (Fig. 10). Overall, the Sen_LEC experiment improves the simulation of the summer mean rainfall over the subregion in terms of PC, TS, and RMSE (Table 3).
As the precipitation over the regions downwind the lake (30.7-31.2° N, 90.81-91.33° E) is significantly influenced by the lake-induced downwind enhancement impacts (Dai et al. 2018a, b;Wu et al. 2019), we select this area to explore the influences of the improved lake-air interactions on modeling the rainfall diurnal variations. The intercomparisons between the diurnal cycles of the regional averaged rainfall from TRMM, ITPCAS, CTRL and Sen_LEC during summer of 2013 are presented in Fig. 11. From the 3-hourly TRMM and ITPCAS data, the rainfall diurnal variation downwind the Lake Nam Co is characterized by a pronounced lateafternoon peak (18 local solar hour, hereafter denoted as LSH) and the relatively larger values from night to early morning. Both the CTRL and Sen_LEC experiments predict a delayed dominant late-afternoon rainfall peak (19 LSH) and overestimate the rainfall during daytime and nighttime. Compared with the CTRL results, Sen_LEC tends to mitigate the general rainfall overestimation and better predict the nighttime rainfall (19-23 LSH) due to the more realistic representations of LST and turbulent water/heat exchanges. However, the WRF-Lake with a calibrated lake scheme still exhibits evident positive biases in modeling the daytime rainfall. Given that the rainfall diurnal variation is a synthesized climatic phenomenon modulated by  Table 3 The pattern correlation (PC), Taylor score (TS), and rootmean-square-error (RMSE) between the simulations and observations for the T 2m and precipitation over the sub-region (30° N-31 the diurnal evolution of boundary layer properties and largescale atmospheric forcing (Dai and Trenberth 2004;Yang et al. 2017), the current daytime rainfall overestimation may be attributed to the uncertainties from the other physical processes in the coupled WRF-Lake. This points out that more assessments on WRF-Lake's physical processes are needed for expanding the model's applicability on TP climatic simulations.
To further understand the underlying mechanisms of the improved T 2m and precipitation simulations, we investigate the disparate atmospheric responses to different LST patterns modeled by the CTRL and Sen_LEC experiments along two vertical transections (Fig. 12), which are indicated by the red lines in Fig. 9f. The southwestnortheast transection crossing the Lake Nam Co and the downwind regions is denoted as Lake-Downwind, and the northwest-southeast transection crossing the Lake Nam Co and the Nyainqentanglha Mountain is denoted as Lake-Mountain. Compared with CTRL, the Sen_LEC experiment mitigates the LST overestimation, weakens the upward water and heat fluxes, and obtains a older and drier overlying atmosphere around Lake Nam Co (Fig. 12a-d).
The lake-induced anomalous cooling and desiccating effects can reinforce each other in generating the stable boundary layer up to ~ 5.8 km and inhibit the convective developments around the Lake Nam Co (Fig. 12e and f). At the diurnal timescale, as the over-lake colder and drier air in Sen_LEC could be advected by the lake breeze and valley wind during daytime, the atmosphere column over the surrounding land and mountainous regions is featured by a stable boundary layer and anomalous katabatic motions. This explains the mitigated rainfall overestimation over these regions in the Sen_LEC run. During nighttime, the colder LST simulated by Sen_LEC weakens the lake-land thermal contrast and the over-lake ascending motions, resulting in less rainfall in the central lake regions (Fig. 9d-f). tions, and f Sen_LEC-CTRL differences, respectively. The black curves in (a-f) outline the locations and ranges of TP lakes. In (f), the blue dashed rectangular denotes the regions downwind the lake (30.7-31.2° N, 90.81-91.33° E), and the red lines represent the transections used for lake-downwind and lake-mountain vertical cross sections in Fig. 12

Summary and discussion
Based on the in-situ observation data at Lake Nam Co, we evaluate and calibrate the coupled WRF-Lake from the perspective of model's performance in simulating the LST and the evolutions of vertical temperature profiles during summer. The effects of the improved LST and lake-air interactions on regional climate modelling are also discussed. Main findings are summarized as follows: During summer, the upper-most water temperature (i.e., the lake temperature at 3 m and 6 m) gradually increases from 3.5° to 12 °C. The whole water column warms slowly and depicts a weak temperature gradient before the onset of stable lake thermal stratification (June 12th). In response to the enhanced solar heating from then on, the epilimnetic temperature rises rapidly with the mixed layer deepening to ~ 26 m, and the simultaneous metalimnion and hypolimnion temperature exhibit fluctuant rising features. However, the coupled WRF-Lake model with default lake scheme parameters predicts excessively warmer LST, earlier establishment of lake thermal stratification and a colder epilimnion during the late summer. By adjusting the initial lake temperature profile, Tdmax, surface roughness length and light extinction coefficient one by one from Sen_LTI to Sen_LEC experiments, WRF-Lake model's capability in reproducing the evolutions of LST and vertical temperature profiles in Lake Nam Co is improved. The RMSE (BIAS) Fig. 10 Same as Fig. 9, but for the comparisons of summer mean rainfall spatial distributions from ITPCAS datasets, CTRL, and Sen_LEC experiments Fig. 11 The diurnal cycles of rainfall averaged over the lake downwind regions from TRMM, ITPCAS, CTRL, and Sen_LEC during summer in 2013. The X axis is local solar hour (LSH) of summer TLake 3m in Sen_LEC experiment is 0.71 °C (− 0.07 °C), reduced from 1.13 °C (0.47 °C) in CTRL run. The onset of lake thermal stratification predicted by Sen_ LEC (June 11th) is also in better agreement with observations (June 12th). Benefited from the more realistic simulations in the near-surface water temperature and lake-air interactions, the calibrated WRF-Lake model mitigates the lake's excessive warming and moistening effects, weakens the convection over the downwind and mountainous regions around Lake Nam Co, and effectively reduces the over-lake T 2m overestimation and the positive precipitation biases therein.
By adopting the model configurations in CTRL and Sen_LEC, we conduct two additional experiments running from May 20th to November 31st 2013 to investigate the performances of the calibrated WRF-Lake during autumn, when the lake-air interactions feature strong upward vapor/ heat transport due to the large positive lake-air temperature/humidity gradients and the enhanced over-lake wind speed (Wang et al. 2020). The onset (end) of lake thermal stratification is defined as the first (last) date from which the temperature differences between 6 and 56 m depths are greater than 1 °C in the subsequent (previous) 4 days without exception. The overestimated TLake 3m during June to mid-July and the earlier onset of the lake thermal stratification in CTRL run are attributed to the unrealistically higher initial water temperature profile ( Fig. 13a and b), while the underestimated TLake 3m since mid-July and the earlier autumnal Fig. 12 The height-latitude lake-downwind (left) and lakemountain (right) cross sections of the simulated summer mean differences (Sen_LEC-CTRL) in the air temperature (top), specific humidity (middle), and equivalent potential temperature (shaded) with the longitudinal circulation (vectors, units: m s −1 ) (bottom) during 2013. The vertical velocity was enlarged by 70 times. The gray shaded areas in (a-f) denote the absolute terrain height above sea level. The sunken areas marked by the blue dashed lines and stars in (a-f) roughly represent the locations of Lake Nam Co destratification are mainly due to the excessive turbulent vapor/heat loss modulated by the large constant surface roughness length. Compared with CTRL, the Sen_LEC experiment well reproduces the surface layer lake temperature, the onset/end of thermal stratification, and the vertical water temperature variations throughout the simulation periods (Fig. 13). From Fig. 14, both CTRL and Sen_LEC generally reproduce the west-east precipitation gradient over the central TP and the westerly wind fields at 500 hPa. For the regions surrounding Lake Nam Co, the local precipitation from TRMM/ITPCAS (figure omitted) is characterized by a maximal center located to the west of the lake, namely the regions downwind the lake (marked by blue rectangle in Fig. 14f). Both experiments predict the maximal rainfall center downwind the lake despite that the precipitation intensity is overestimated. Compared with CTRL, the Sen_LEC experiment with a smaller parameterized surface roughness length effectively restricts the upward water/heat transport and thus mitigates the rainfall overestimation downwind the lake. The 500-hPa wind field differences therein is featured by divergent patterns and the regional-averaged rainfall overestimation is reduced by 32%, from 1.18 mm·day −1 in CTRL to 0.8 mm·day −1 in Sen_LEC (Fig. 14f). The above comparisons demonstrate that the calibrated WRF-Lake is superior in modeling the autumnal lake-air interactions and the associated lake-effected regional climate surrounding the large lake basins over the central TP.
Although substantial efforts have been devoted to obtain better LST and lake-air interaction simulations, the calibrated WRF-Lake model (Sen_LEC) still exhibits disadvantages in predicting the rapid LST warming from 4.42° to 9.29 °C during June 11th-June 23th and the high-frequency oscillating variations of the metalimnion temperature (Fig. 3). These discrepancies may be attributed to both the lake model's forcing uncertainties inherited from WRF and the 1-D physical simplification by neglecting the flowdependent heat redistribution processes within lakes. On one hand, during the rapid LST rising period in early-summer, the Sen_LEC experiment overestimates the 10-m wind speed at the Nam Co water temperature station with a maximal (mean) BIAS of 2.98 m s −1 (0.07 m s −1 ), which can result in the excessive upward heat loss and therefore the underestimation of LST increase. As documented by Xu et al. (2016), the surface wind speed overestimation in WRF-Lake can be partly attributed to the smoothing orography representation over complex-terrain regions (i.e., Nyainqentanglha Mountainous area). Hence, future work needs to parameterize the unresolved turbulent orographic effects to better reproduce the surface wind speed (Jimenez et al. 2014;Zhou et al. 2019) and thus mitigate the LST underestimation. Another possible reason for the inadequate modeled LST increase at the Nam Co water temperature station is the missing heat redistribution processes related to the 3-D thermohydrodynamics. When the Lake Nam Co begins to be heated from its original thermally mixed state in early spring, LST over the shallower shore areas rises faster in response to the radiative forcing because the heat capacitiy therein is smaller than that in the central deep lake. This would lead to the horizontal nearshore-offshore temperature gradient across the isobath (Fig. 15), which has also been reported in Fig. 14 Same as Fig. 9, but for the comparisons of mean rainfall and 500-hPa wind spatial distributions during autumn (September to November) Fig. 15 The mean water temperature at 0.05 m depth (TLake 0.05 m ) from Sen_LEC over Lake Nam Co during June in 2013. The vectors are the synchronous mean wind fields at 500 hPa (unit: m s −1 ) from ERA5 reanalysis. The black aster represents the location of the lake temperature monitoring site Lake Victoria, Laurentian Great Lakes, and many other large deep-water bodies with sloping bottom bathymetry (Bai et al. 2013;Nyamweya et al. 2016;Rao et al. 2004). The cross-isobath temperature gradient tends to transport vast heat from the warm shore to mid-lake areas through both the horizontal temperature diffusion and the gyre-related advection induced by the baroclinic and bathymetry effects (Beletsky and Schwab 2001). In addition, accompanied by the outbreak and enhancement of the southwest monsoonal circulations (Figs. 9 and 15), the warmer epilimnetic water in the southwest shallower lake shores may experience northward advection processes and thus helps to generate the rapid LST warming at the Nam Co water temperature station. Such horizontal heat transportation processes driven by the temperature-current interactions and the prevailing surface wind have been proved to play critical roles in determining the spatiotemporal variabilities of LST over many large lakes worldwide (Song et al. 2004;Thiery et al. 2014;Xue et al. 2016). However, constrained by the 1-D structure and its simplified representation in lake mixing mechanism, WRF-Lake is unable to depict the above lake thermohydrodynamics and the rapid LST warming at deep-water Nam Co buoy since mid-June. Regarding the pronounced highfrequency oscillation of metalimnion temperature, i.e., the abrupt thermocline jump episodes during August 10th-20th, it is likely attributed to the intensive internal wave activities in the upper portion of thermocline Wang et al. 2019), but the current 1-D lake scheme does not include any representation of internal waves. It remains interesting to demonstrate the benefits of employing a fully coupled 3-D lake component in the regional coupled climate models over the lake-rich TP.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.