Evaluation of the Effect of Stability Schemes on the Simulation of Land Surface Processes at a Western Tibetan Site

The surface fluxes calculated in land surface models (LSMs) are sensitive to the determination of the stability parameter. Further, calculation of the surface fluxes over the Tibetan Plateau (TP) is crucial in the simulation of regional and global weather and climate. In this study, we use 2-year micrometeorological data measured from Shiquanhe, located in the western TP, to evaluate the performance of the widely used Noah LSM with five stability parameterization schemes. Results show that all five stability parameterization schemes can generally reproduce the observations, but the scheme proposed by Li has the smallest bias. The reason is that Li’s scheme is more accurate under the unstable condition, and the surface layer at Shiquanhe is mostly unstable. Further, the four non-iterative schemes show an advantage in terms of their computational efficiency compared to the iterative scheme adopted by the Noah LSM.


Introduction
The Tibetan Plateau (TP) plays an important role in regional and global weather and climate through its thermodynamic and mechanical forcing [1][2][3]. To understand the thermodynamic forcing of the TP and to accurately estimate the exchange of mass and energy over the TP, field experiments have been conducted in recent decades, and many studies revealed that the land-atmosphere interactions over the TP greatly affected the climate [4][5][6][7][8][9][10]. In recent years, studies showed that the TP was a "hot spot" for landatmosphere interaction. Therefore, accurately estimating the exchange of mass and energy has great importance for understanding the TP and global climate. A land surface model (LSM) alone, or coupled to an atmospheric model, has been widely used to simulate the surface mass and energy flux and the climate over the TP. However, apparent biases in LSMs can be detected [11][12][13]. These surface flux biases in LSMs constitute one of the uncertainty sources in numerical weather or climate models, which can then lead to low skill in the simulation of rainfall and temperature over the TP [14]. Hence, improving the performance of the parameterization of surface layer fluxes in LSMs is essential to the accuracy of weather and climate model simulations over the TP.
By comparing different land surface process models and related parameterization schemes, it was found that the near-surface parameterization scheme was the most essential for the simulated flux over the TP. Therefore, how to reasonably parameterize the nearsurface layer in the model was particularly important. For parameterizing the near-surface layer, the Monin-Obukhov similarity theory (MOST) [15] has been widely applied in LSM [16]. According to MOST, land surface fluxes of momentum and heat can be calculated via the surface wind speed, temperature gradient, and bulk transfer coefficient [17]. While the first two variables can be provided by either observation or numerical modeling, estimation of the bulk transfer coefficient is crucial for the calculation of surface fluxes. In an LSM, the bulk transfer coefficient is usually determined by the aerodynamic roughness length, thermal roughness length, profile function, and stability parameter. The bulk transfer coefficient is not sensitive to the profile function, but is sensitive to roughness lengths, as shown in a real case simulation using the Noah LSM [18]. The value of the aerodynamic roughness length is dependent on the land cover properties, while various methods have been proposed to calculate the thermal roughness length [19][20][21][22].
The last parameter, i.e., the stability parameter, is important in determining the bulk transfer coefficient, and various schemes have been proposed to calculate its value [23]. A numerical iteration is usually used to solve stability equations [24,25]. In recent years, various non-iterative schemes have been proposed to replace the iterative method [26][27][28][29][30]. For example, Yang [26] derived an exact solution of the stability parameter equation for a stable surface layer and proposed an approximate analytical solution for an unstable surface layer by using a bulk Richardson number and profile functions. Li [27], Sharan [30], and Wouters [29], using theoretical calculations or observational data from certain stations, performed non-iterative calculations or linear fitting of stability parameters and found that the stability can be calculated more accurately, thereby improving the flux simulation. However, the research of Sharan et al. is only an offline calculation, and it is not really applied to the land surface model, so it is difficult to evaluate the impact on the overall land surface process. In addition, the fitting coefficients given in this type of study may be related to the site, making it difficult to assess its applicability over the TP. Non-iterative schemes are computationally efficient theoretically, but how well the surface layer over the western TP can be simulated by these schemes remains unknown.
The objective of this work is twofold: (1) to assess the sensitivity of the land surface flux and other variables to different parameterizations of stability in the Noah LSM and (2) to evaluate the performance and applicability of four currently developed non-iterative schemes on land surface processes in the Noah LSM against field observations at Shiquanhe in the western Tibetan Plateau. The Noah LSM is selected because it is widely used and has been adopted for operations and research in National Centers for Environmental Prediction (NCEP) weather and climate predictions models and relevant data assimilation systems. Following the introduction, descriptions of the data, the model, and the experimental design are all presented in Section 2. The main results are given in Section 3, with some further discussion provided in Section 4. Finally, a summary is presented in Section 5.

Data
The micrometeorological data obtained from Shiquanhe (32.50 • N, 80.08 • E, 4279.3 m above sea level), located in the western TP, are used to evaluate the performance of the Noah LSM. The west and south sides of Shiquanhe are surrounded by the Himalayas, while north of Shiquanhe is the Kailash Range. The site is located in the transition zone between the monsoon region and the non-monsoon region, affected by the India monsoon to some extent, categorized as a mid-latitude arid continental temperate climate. The annual mean surface air temperature in Shiquanhe is 0.35 • C, and the total precipitation is 73 mm. Figure 1 shows the site position, surface condition, sensors, and the meteorology element variations during the experiment. The land surface consists of sand and some small areas of gravel, which is typical in the western TP. The site was established in September 2013 and continues to operate today. The data from 1 September 2013 to 1 September 2015 are used in this study, including the 30-min mean downward shortwave and longwave fluxes, sensible heat flux, latent heat flux, soil moisture and temperature, wind speed, relative humidity, air pressure, air temperature, and precipitation. Information on the sensors used -to measure these parameters is provided in Table 1.   Time series of wind speed, wind direction, air temperature, specific humidity, air pressure, and precipitation at Shiquanhe from 1 September 2013 to 1 September 2015, respectively. EddyPro 6. 2.1 (LI-COR Inc., Lincoln, USA, 2017) was used in quality control and turbulent flux calculation to obtain good performance fluxes (such as the momentum flux, the sensible heat flux, and the latent heat flux) from the turbulent data collected by EC150. The steps and flux corrections are as follows: (1) set 30 min as the flux averaging interval to compute turbulent fluxes, (2) remove the spikes of the raw time series with a criterion of (( − 4 ) < ( ) < ( + 4 )), where x(t) denotes the measured variables (i.e., three-dimensional wind speed u, v, and w), and and σ are the mean and standard deviation of the data block, (3) skip the turbulent data block that missed more than 10%, (4) use the double rotation method for the sonic anemometer tilt correction, (5) use the WPL (Webb, Pearman, and Leuning) method for density corrections, (6) use high-/low-frequency spectral corrections to compensate for flux losses, (7) remove the flux data in rainy (snowy) or foggy (hazy) weather, as under such conditions, eddy covariance systems work improperly [31].
In addition, the stability parameter ξ (= z/L) based on observation was computed with Equations (1)-(3), and it is used as the criterion to identify the stable/unstable condition when the sensitivity of the simulations to the atmospheric stability is discussed. ξ > 0 was considered as stable stratification, while ξ < 0 was considered unstable.    interval to compute turbulent fluxes, (2) remove the spikes of the raw time series with a criterion of ( X − 4σ < x(t) < X + 4σ ), where x(t) denotes the measured variables (i.e., three-dimensional wind speed u, v, and w), and X and σ are the mean and standard deviation of the data block, (3) skip the turbulent data block that missed more than 10%, (4) use the double rotation method for the sonic anemometer tilt correction, (5) use the WPL (Webb, Pearman, and Leuning) method for density corrections, (6) use high-/lowfrequency spectral corrections to compensate for flux losses, (7) remove the flux data in rainy (snowy) or foggy (hazy) weather, as under such conditions, eddy covariance systems work improperly [31].
In addition, the stability parameter ξ (= z/L) based on observation was computed with Equations (1)-(3), and it is used as the criterion to identify the stable/unstable condition when the sensitivity of the simulations to the atmospheric stability is discussed. ξ > 0 was considered as stable stratification, while ξ < 0 was considered unstable.
where z is the height of EC150 installed, u * is the friction velocity, θ * is the potential temperature scale, L is the Monin-Obukhov length, u', v', w', and θ' are the deviations of the observed horizontal and vertical wind speeds and potential temperatures, respectively, k is the von Kármán constant, and g is the acceleration of gravity.

Model and Experimental Design
Our aim in this paper is to evaluate the performance of four non-iterative stability parameterization schemes in simulating the surface layer at Shiquanhe, located in the western TP, by using the Noah LSM. Specifically, the Community Noah Land Surface Model, version 2.7.1, was used in this study ( [32]), which we refer to simply as the Noah LSM. The Noah LSM has been tested for different land surfaces and has also been coupled to several numerical weather models, including the Weather Research and Forecasting model [33][34][35][36]. This version of the Noah LSM is a standalone, uncoupled, 1-D column version, used to execute single-site land surface simulations at Shiquanhe. It requires nearsurface atmospheric forcing data, initial soil moisture, and surface parameters. In addition to stability parameters, vegetation and soil parameters need to be input. In this study, the default parameters of the model were directly used to obtain the surface parameters, which are obtained by looking up the table. Among them, SOILTYP (variable name for soil type) and VEGTYP (variable name for vegetation) are the soil and vegetation parameters, respectively. There are many studies on the importance of soil and vegetation parameters. To avoid the influence of soil and vegetation parameters on the research results, this study directly used the method of a look-up table. The land surface at Shiquanhe consists of sand and some small areas of gravel. Thus, the values of soil type (SOILTYP) and vegetation type (VEGTYP) were set to 1 and 11, respectively. The soil has five layers, at 0.05, 0.1, 0.2, 0.4, and 0.8 m, respectively. To directly compare with the observation, the soil stratification is consistent with the observation. The value of the aerodynamic roughness length is 0.108, which was calculated from measured data at Shiquanhe. The values of soil parameters were obtained from the model's look-up tables. The initial soil temperature and moisture were provided by site observation at 00Z 1 September 2013. The near-surface atmospheric forcing data are composed of the observed 30-min mean surface wind speed, temperature, relative humidity, pressure, downward shortwave and longwave flux, and precipitation. The simulation period was two years, from 1 September 2013 to 1 September 2015. The simulation of the first year was regarded as a "spin-up" period, and thus only the simulation of the last year is analyzed in this paper.
The surface layer parameterization in the Noah LSM is based on MOST, and the following equations are used to calculate surface fluxes: where τ and H are surface momentum and sensible heat fluxes, respectively; u * and θ * are friction velocity and temperature, respectively; ρ and c p are the air density and the air heat capacity, respectively; u and θ are wind speed and air temperature at reference height, respectively; C D and C H are surface exchange coefficients for momentum and heat, respectively; and θ s is the surface air potential temperature. Based on MOST, u * and θ * are written as where k is the von Kármán constant; ψ M and ψ H are profile functions for momentum and heat, respectively; L is the Monin-Obukhov length; ξ = z/L is the stability parameter; z is the reference height; ξ 0 (=z 0m /L) and ξ T (=z T /L) are both temporary variables and have no clear physical meaning, while z 0m and z T are aerodynamic and thermal roughness lengths, respectively; g is the acceleration of gravity; and T is absolute mean air temperature.
In the Noah LSM, the stability parameter is estimated by an iterative computational method. We adopted the following four non-iterative schemes to replace the iterative method in the Noah LSM in this study: (1) Yang [26] et al. (Y01) proposed the following equations for the stability parameter by using the bulk Richardson number: For a stable surface layer, the exact solution of the stability parameter equations can be derived directly from the linear function, while only an approximate analytical solution is obtained for an unstable surface layer.
(2) Li [27] et al. (L10) suggested the following form of ξ for unstable, weak-stable, and strong-stable conditions: where A i and B i (i = 0, 1, 2) are regression coefficients whose values are affected by the observation height as well as the aerodynamic and thermal roughness lengths.
(3) Wouters [29] et al. (W12) proposed the following scheme for calculating the stability parameter after considering the effect of the roughness sublayer: where f i is a multiple regression function of the bulk Richardson number, reference height, and aerodynamic and thermal roughness lengths, Ri b , and t is the value to separate strongstable and weak-stable conditions.
(4) Sharan and Srivastava [30] (S14) proposed a semi-analytical scheme to parametrize the stability parameter: where R i is the gradient Richardson number; a i = a i (Ri), i = 0, 1, 2; and g is a multiple regression function.
To evaluate the effect of the stability parameter on the simulation of the surface layer at Shiquanhe, we ran the Noah LSM with its default iterative scheme, and the four noniterative schemes of Y01, L10, W12, and S14. Thus, there are five 1-year simulations, which are referred to as Noah, Noah_Y01, Noah_L10, Noah_W12, and Noah_S14, according to the scheme used.

Results
To evaluate the simulation effect of the Noah MP model more accurately, the statistical analysis method is used in this study. The statistics used are as follows: (1) Coefficients of correlation (R). (2) Model efficiency index (NSE, Nash-Sutcliffe forecasting efficiency) [37,38], which is used to evaluate the predictability of the model. Its value ranges from negative infinity to 1.0. When it is equal to 1.0, the simulated value is exactly the same as the observed value, indicating that the model is perfect. (3) The mean bias error (MBE). (4) Root mean square error (RMSE), which can directly give the deviation between the simulated value and the observed value. The smaller the value, the closer the simulated value to the observed value. The calculation formulas of the statistics are as follows: where N, P i , and O i are the sample number, the simulated value, and the observed value, respectively. Figure 2 shows the observed and simulated monthly mean diurnal cycles of the surface fluxes for momentum and both sensible and latent heat. To show the monthly variation of the whole year in turn, the data from January to August in the figure are from 2015, and the data from September to December in the figure are from 2014. Consistent with people's general thinking habits, the chart is arranged in the order of January to December. Table 2 lists R, NSE, MBE, and RMSE of the simulated fluxes of momentum and both sensible and latent heat with observations. All the simulations capture the basic features of the diurnal cycle and seasonal variation in the observed fluxes. All the correlations of the observations with the simulated fluxes are higher than 0.6. The correlation coefficients indicate that the Noah LSM shows the best skill for simulating the diurnal cycle and seasonal variation in sensible heat, but the lowest skill for latent heat ( Table 2). All the simulations overestimate the momentum flux, especially for spring. The maximum bias of the simulated momentum flux is produced by Noah_Y01, with an MBE of 0.02 kg m −1 s −2 , while the minimum bias is found for Noah_L10, with the lowest values for both MBE and RMSE. All simulations overestimate sensible heat for all months. Noah_L10 also produces the minimum bias of sensible heat. There are no large differences among the other four simulations of sensible heat, whose MBE and RMSE are higher than 20 and 60 W m −2 , respectively. The differences in latent heat flux are small among the five simulations, whose MBE and RMSE for latent heat flux are higher than 0.1 and 10 W m −2 , respectively. While the simulations capture the general feature of sensible heat flux, they show a large bias for some months-June, for example. The reason may be that the sensible heat flux simulated by the model is not only related to the near surface parameterization scheme, but also closely related to other physical processes, such as soil water and heat transfer. Moreover, the effects of surface parameters such as albedo, emissivity, and thermal roughness cannot be ignored. These physical processes and parameterization schemes are not the focus of this study, so they are all default parameters of the Noah model, which may cause large deviation. Especially in June, the western Tibetan Plateau also entered the rainy season, the precipitation increased, and the soil moisture increased, which led to great changes in the thermodynamic parameters, which may have a large deviation from the default parameters of the model, and eventually led to the deviation of the sensible heat flux simulation.    The above analyses indicate that estimating the stability parameters accurately can improve the simulated sensible and latent heat fluxes, and Noah_L10 performs the best. On the other hand, there are no large differences in the simulated latent heat. One possible reason is that latent heat is small at Shiquanhe. Noah outperforms Noah_Y01 for all the turbulent fluxes, but it is beaten by Noah_L10. Noah's skill is compared to other non-iterative schemes. Figure 3 shows the observed and simulated monthly mean diurnal cycles of upward shortwave and longwave radiative fluxes, and Table 3 lists R, NSE, MBE, and RMSE of the simulations of radiative flux with observations of radiative flux. All the simulations perform very well for the variation in upward shortwave and longwave fluxes, with correlation coefficients nearing 1.0. The differences in upward shortwave flux among the simulations are small, which is ascribed to the similar simulation results of the surface albedo. As the surface albedo is significantly affected by snowfall in winter and rainfall in summer, the Noah LSM produces a larger bias in winter and summer (Figure 3a). There are large diversities in the simulated upward longwave flux, including large positive and negative biases. Noah_Y01 has the largest bias, with an MBE of −19.48 W m −2 and an RMSE of 34.58 W m −2 ; Noah_L10 produces the minimum bias, with an MBE of 3.19 W m −2 and an RMSE of 12.62 W m −2 . The non-iterative schemes outperform Noah_L10 in simulating upward longwave flux, except for Noah_Y01. The Noah LSM calculates the upward longwave radiation (ULR) from the surface temperature using the Stefan-Boltzman law and the modification of surface emissivity and surface temperature directly influences the ULR [39]. It can be inferred that the emissivity and surface temperature are sensible to stability parameter schemes. The non-iterative schemes outperform Noah_L10 in simulating upward longwave flux, except for Noah_Y01. The Noah LSM calculates the upward longwave radiation (ULR) from the surface temperature using the Stefan-Boltzman law and the modification of surface emissivity and surface temperature directly influences the ULR [39]. It can be inferred that the emissivity and surface temperature are sensible to stability parameter schemes.   Figure 4 presents the observed and simulated monthly mean diurnal cycles of the soil temperature at 10 and 40 cm, separately. The soil temperature at shallow depths shows an obvious diurnal cycle, which is not observed at the deeper depth (40 cm). Soil temperature is higher in summer but lower in winter. All the simulations capture these features, and the correlations between the observation and simulations are higher than 0.88. All the simulations produce a cold bias of soil temperature, except for Noah_L10 for  Figure 4 presents the observed and simulated monthly mean diurnal cycles of the soil temperature at 10 and 40 cm, separately. The soil temperature at shallow depths shows an obvious diurnal cycle, which is not observed at the deeper depth (40 cm). Soil temperature is higher in summer but lower in winter. All the simulations capture these features, and the correlations between the observation and simulations are higher than 0.88. All the simulations produce a cold bias of soil temperature, except for Noah_L10 for the soil temperature at 5, 10, and 20 cm ( Table 4). The simulations with the non-iterative schemes generally outperform Noah, except Noah_Y01. Noah_L10 performs best. These results also indicate that the stability parameter affects the simulation of soil temperature. Theoretically, soil temperature can affect surface sensible heat flux, latent heat flux, and net longwave radiation, and vice versa. When the change in stability parameters leads to the change in sensible heat flux, the surface energy balance also changes, which will cause the change in soil temperature. Li et al. [40] evaluated the Noah LSM with multi-parameterization (Noah-MP) for soil temperature simulation and pointed out that Noah-MP generally underestimates soil temperature, especially during the cold season. Further, the simulation uncertainty is greater in the cold season and for deep soil layers. At Amdo in the central Tibetan Plateau, Noah-MP underestimates soil temperature from the top to the bottom layers in the monsoon season, and the underestimation might relate to the initiation of soil state variables [41].  the change in soil temperature. Li et al. [40] evaluated the Noah LSM with multi-parameterization (Noah-MP) for soil temperature simulation and pointed out that Noah-MP generally underestimates soil temperature, especially during the cold season. Further, the simulation uncertainty is greater in the cold season and for deep soil layers. At Amdo in the central Tibetan Plateau, Noah-MP underestimates soil temperature from the top to the bottom layers in the monsoon season, and the underestimation might relate to the initiation of soil state variables [41].   The observed and simulated monthly mean diurnal cycles of soil moisture indicate that all the simulations have a positive bias in simulated soil moisture (figures not shown). They also cannot capture the seasonal variation in soil moisture. There are no large differences in the simulated soil moisture among the simulations. Thus, the stability parameter is not a key parameter for the simulation of soil moisture at Shiquanhe by the Noah LSM. The observed and simulated monthly mean diurnal cycles of soil moisture indicate that all the simulations have a positive bias in simulated soil moisture (figures not shown). They also cannot capture the seasonal variation in soil moisture. There are no large differences in the simulated soil moisture among the simulations. Thus, the stability parameter is not a key parameter for the simulation of soil moisture at Shiquanhe by the Noah LSM. In theory, soil moisture can affect surface evaporation, or latent heat flux. According to the energy balance in the land surface model, soil moisture affects latent heat flux, changes the surface energy balance, and causes the change in sensible heat flux, and vice versa, that is, the change in energy balance will affect the latent heat flux, and then may change the soil moisture. At Shiquanhe, the difference in sensible heat flux caused by the difference in stability parameters can change the surface energy balance, but the area is extremely dry and the soil moisture is very small, so the change in energy balance has little effect on the soil moisture. Ye et al. also found that the simulation error of the Noah LSM is larger when the soil moisture is small and considered that it is due to the observation error and the inaccurate root system distribution parameters of vegetation [42].

Sensitivity of the Simulations to the Stability Condition
As stability is an important factor for estimating stability parameters, how well do the schemes perform for the stable and unstable conditions? Figure 5 shows the observed and simulated fluxes for momentum (upper panels) and sensible heat (lower panels) for the stable condition. The simulations generally have better skill for momentum than that for sensible heat, except for Noah_Y01, which has the largest bias for both momentum and sensible heat fluxes. There is no large difference in both momentum and sensible heat fluxes among the other simulations, while Noah_S14 and Noah have the highest correlation with the observation for momentum and sensible heat fluxes, respectively. Based on both regression coefficients and correlation coefficients, Noah_S14 and Noah_L10 perform best for momentum flux and sensible heat flux, respectively. Overall, these three schemes (Noah_L10, Noah_S14, Noah) are potentially applicable to the simulation of a stable surface layer at Shiquanhe. sensible heat fluxes. There is no large difference in both momentum and sensible heat fluxes among the other simulations, while Noah_S14 and Noah have the highest correlation with the observation for momentum and sensible heat fluxes, respectively. Based on both regression coefficients and correlation coefficients, Noah_S14 and Noah_L10 perform best for momentum flux and sensible heat flux, respectively. Overall, these three schemes (Noah_L10, Noah_S14, Noah) are potentially applicable to the simulation of a stable surface layer at Shiquanhe.  Figure 6 shows the observed and simulated fluxes for momentum (upper panels) and sensible heat (lower panels) for the unstable condition. Compared to sensible heat flux, the simulations have better skill for momentum flux, except for Noah_Y01. All the simulations overestimate sensible heat flux. Noah and Noah_L10 perform very well for momentum flux, as the regression line is generally identical to the 1:1 line. Noah_L10 also  Figure 6 shows the observed and simulated fluxes for momentum (upper panels) and sensible heat (lower panels) for the unstable condition. Compared to sensible heat flux, the simulations have better skill for momentum flux, except for Noah_Y01. All the simulations overestimate sensible heat flux. Noah and Noah_L10 perform very well for momentum flux, as the regression line is generally identical to the 1:1 line. Noah_L10 also performs best for sensible heat flux, while some apparent bias is noticed. Thus, Noah_L10 is suitable for an unstable surface layer.
It can be seen from Figures 2b and 6 that although the sensible heat fluxes simulated by the model consistently exceed the observation, there are still weak differences between the simulated sensible heat fluxes. It seems that the stability parameter schemes can affect the simulated surface heat flux. However, the sensible heat flux simulated by the model is not only related to the near surface parameterization scheme, but also closely related to other physical processes, such as soil water and heat transfer. Moreover, the effects of surface parameters such as albedo, emissivity, and thermal roughness cannot be ignored. These physical processes and parameterization schemes are not the focus of this study, so they are all default parameters of the Noah model, which may cause large deviation. Especially in June, the western Tibetan Plateau also enters the rainy season, the precipitation will increase, and the soil moisture will increase, which lead to great changes in the thermodynamic parameters, which may have a large deviation from the default parameters of the model, and eventually lead to the deviation of the sensible heat flux simulation. Hu et al. [43] pointed out that vegetation land surface processes are the primary source of uncertainty in heat flux simulation, and that the more arid the area is, the larger the deviation of sensible heat flux simulated by Noah. Gao et al. found that with Noah-MP, the surface sensible and latent heat fluxes are better simulated in the monsoon season as well in the central Tibetan Plateau [41]. performs best for sensible heat flux, while some apparent bias is noticed. Thus, Noah_L10 is suitable for an unstable surface layer. It can be seen from Figures 2b and 6 that although the sensible heat fluxes simulated by the model consistently exceed the observation, there are still weak differences between the simulated sensible heat fluxes. It seems that the stability parameter schemes can affect the simulated surface heat flux. However, the sensible heat flux simulated by the model is not only related to the near surface parameterization scheme, but also closely related to other physical processes, such as soil water and heat transfer. Moreover, the effects of surface parameters such as albedo, emissivity, and thermal roughness cannot be ignored. These physical processes and parameterization schemes are not the focus of this study, so they are all default parameters of the Noah model, which may cause large deviation. Especially in June, the western Tibetan Plateau also enters the rainy season, the precipitation will increase, and the soil moisture will increase, which lead to great changes in the thermodynamic parameters, which may have a large deviation from the default parameters of the model, and eventually lead to the deviation of the sensible heat flux simulation. Hu et al. [43] pointed out that vegetation land surface processes are the primary source of uncertainty in heat flux simulation, and that the more arid the area is, the larger the deviation of sensible heat flux simulated by Noah. Gao et al. found that with Noah-MP, the surface sensible and latent heat fluxes are better simulated in the monsoon season as well in the central Tibetan Plateau [41]. Figure 7 shows the frequency of the stability parameter ξbased on turbulence data at Shiquanhe. The values of ξ with a frequency higher than 0.001 range from −2.0 to 4.0, and frequent values range from −0.1 to 0.1; the total frequency of the unstable condition (ξ<0, the red bar) is 69.3%, and the total frequency of the stable condition (ξ> 0, the blue bar) is 30.7%, which indicates the surface layer at Shiquanhe is mostly unstable. As Noah_L10 performs best for the unstable condition, it is reasonable that it has the best skill in simulating the surface layer at Shiquanhe, as discussed in the above section.  Figure 7 shows the frequency of the stability parameter ξbased on turbulence data at Shiquanhe. The values of ξ with a frequency higher than 0.001 range from −2.0 to 4.0, and frequent values range from −0.1 to 0.1; the total frequency of the unstable condition (ξ < 0, the red bar) is 69.3%, and the total frequency of the stable condition (ξ > 0, the blue bar) is 30.7%, which indicates the surface layer at Shiquanhe is mostly unstable. As Noah_L10 performs best for the unstable condition, it is reasonable that it has the best skill in simulating the surface layer at Shiquanhe, as discussed in the above section.

How Stability Schemes Affect Simulation Results
The simulation results show that varied stability schemes have important effects on the simulated turbulent flux, net radiation, and soil temperature and humidity. How does

How Stability Schemes Affect Simulation Results
The simulation results show that varied stability schemes have important effects on the simulated turbulent flux, net radiation, and soil temperature and humidity. How does stability affect these variables? According to the stability parameterization scheme, the stability directly affects the turbulent exchange coefficient, which in turn affects the turbulent flux. The difference between the turbulence exchange coefficients calculated with different stability schemes will lead to differences in the simulated turbulent fluxes (Figure 2). The difference will be conducted for radiation and soil temperature and moisture by the energy balance between land and air in the LSM. However, due to the extremely dry climate at Shiquanhe selected in this study, the impact on soil temperature is obvious, and the impact on soil moisture is weak.

Execution Time
As a climate model usually runs for many decades, such a model needs to adopt a scheme with a short execution time. To test the execution time for the five schemes, we ran a 2-year simulation forced by the observed surface atmospheric elements, repeated 10 times, with the five schemes. The setups for the simulations were the same. The execution times for the 20-year simulations were 20.46, 17.37, 17.81, 17.76, and 20.06 s for Noah, Noah_Y01, Noah_L10, Noah_W12, and Noah_S14, respectively. It can be seen the noniterative schemes are generally more computationally efficient compared to the iterative scheme. It seems that under the existing calculation conditions, there is no significant benefit to replace the iterative scheme with the non-iterative scheme in the land surface model. However, from a large number of model evaluation results, we can see that the existing land surface model still has defects, and its parameterization scheme needs to be further improved. Different researchers hold different opinions on parameterization schemes of which physical processes should be improved and how to improve them. For the Tibetan Plateau, the current simulation studies confirmed that the surface parameterization scheme has the most important influence on the flux simulation. Therefore, improving the surface parameterization scheme is a hot topic in the study of plateau land surface processes and climate.

Conclusions
The TP plays an important role in regional and global weather and climate via its thermodynamic and mechanical forcing. Thus, the exchange of energy and mass between the surface and the atmosphere over the TP is vital for understanding its weather and climate effects. Compared to the central and eastern TP, only a few stations exist in the western TP and there is a lack of long-term flux observations. Thus, a high-quality land surface simulation is imperative for both research and operational communities. In this study, we used 2-year flux data obtained from a site at Shiquanhe, located in the western TP, to evaluate the performance of the widely used Noah LSM with different schemes to estimate the stability parameters for the surface layer.
The Noah LSM does a good job in simulating the surface fluxes of momentum, sensible heat, and upward longwave radiation, as well as soil temperature, at Shiquanhe, which are sensitive to the scheme used to estimate the surface layer stability parameter; thus, the simulation ability of the land surface model can be improved by improving the stability parameterization scheme. The four non-iterative schemes employed in this study show an advantage in terms of their computational efficiency compared to the iterative scheme adopted by the Noah LSM. According to the Nash efficiency coefficient, the scheme proposed by Li [27] performs best and is better than the iterative scheme, and its simulation efficiency of momentum flux, sensible heat flux, upward longwave radiation, and 10 and 40 cm soil temperatures is 1%, 88%, 8%, 6%, and 13% higher than Noah's default scheme, respectively, but its simulation efficiency of the 5 cm soil temperature is lower than Noah's default scheme. Further, it is suitable for an unstable surface layer, which occurs frequently at Shiquanhe, and has potential application value to improve the land surface model of the Tibetan Plateau.
This study indicates that a non-iterative scheme for estimation of the stability parameter can improve the performance of an LSM. The values of the coefficients in the equations used to estimate the stability parameter are mostly based on empirical or regression methods, and thus profile functions may play a role in estimation of the stability parameter. However, the underlying surface at Shiquanhe is the typical underlying surface in the western Tibetan Plateau, which is well representative of the western Tibetan Plateau. The conclusions obtained in this study may be limited to the western TP, which is a typical alpine desert region. Whether it can be applied to the central and eastern TP with its alpine meadow environment is deserving of further study.