In this paper, the impact of urban surface roughness length parameterization scheme on the atmospheric environment simulation over Beijing has been investigated through two sets of numerical experiments using the Weather Research and Forecasting model coupled with the Urban Canopy Model. For the control experiment (CTL), the urban surface parameterization scheme used in UCM is the model default one. For another experiment (EXP), a newly developed urban surface parameterization scheme is adopted, which takes into account the comprehensive effects of urban morphology. The comparison of the two sets of simulation results shows that all the roughness parameters computed from the EXP run are larger than those in the CTL run. The increased roughness parameters in the EXP run result in strengthened drag and blocking effects exerted by buildings, which lead to enhanced friction velocity, weakened wind speed in daytime, and boosted turbulent kinetic energy after sunset. Thermal variables (sensible heat flux and temperature) are much less sensitive to variations. In contrast with the CTL run, the EXP run reasonably simulates the observed nocturnal low-level jet. Besides, the EXP run-simulated land surface-atmosphere momentum and heat exchanges are also in better agreement with the observation.

1. Introduction

With the intensified human activities and urbanization processes, the impacts of urbanization on atmospheric environment (e.g., urban heat islands and urban breezes) at the local to regional scale have drawn considerable attention during recent years [1, 2]. Many a mesoscale modeling study has been performed to investigate the urban issues and much progress has been made in understanding the physical processes such as the turbulent exchanges over urban surfaces [3, 4]. The variation of turbulent exchanges, which largely determines the structure of urban boundary layer and the diffusion of contaminant plumes, is especially important for the evolution of urban atmospheric environment [5, 6]. However, partly owing to the complexity of urban landscapes, the turbulent exchange between urban surface and atmosphere has not yet been fully understood and properly represented in most mesoscale atmospheric models.

Generally, turbulent exchanges over urban areas are much more complicated and geometry-dependent [7], as cities have distinctive roughness element geometries like high-rise buildings and crisscrossed pavements [8]. To better describe the turbulent exchange process over urban surfaces, the effects of surface morphology (i.e., urban form) must be considered. The aerodynamic roughness length (), which gives a measure of the capacity of the surface elements in absorbing momentum, is one of the fundamental parameters in atmospheric models to link up the turbulent exchange process with surface morphology. Moreover, is also a crucial parameter for the calculation of aerodynamic resistance to momentum, heat, and moisture transfers, which in turn affects the land surface-atmosphere interactions [9].

A number of numerical simulations have demonstrated the importance of implementations in the atmospheric environment. Early in 1980s, atmospheric modeling studies have shown that changes in surface roughness length could lead to changes in surface sensible and latent heat fluxes, which later become a trigger for changes in near-surface climate and atmospheric circulation patterns [1013]. For cases with mesoscale atmospheric models, the simulation accuracy of urban atmospheric environment is also improved by assigning a more reasonable value to over urban areas [1417].

Theoretically, has been found to be closely correlated with the geometry of surface elements including mean height, frontal area index (frontal area of roughness elements per unit ground area), plan area index (fraction of unit ground area covered by roughness elements), and the arrangement of roughness elements and is much heightened over urban surfaces because of the drag exerted by buildings and other infrastructures [1822]. Among these parameterization schemes, the scheme proposed by Macdonald et al. [21] (hereafter as MGH98 scheme) has been incorporated into the Urban Canopy Model (UCM) in the Weather Research and Forecasting model (WRF). Noting the importance of the variability of building height in describing the geometrical complexities for real urban areas [23], Cao et al. [24] proposed a new parameterization scheme which is basically from Shao and Yang [22] but additionally considers the dependency of upon building height variability (hereafter as SY08N scheme). Building height variability makes tall buildings experience less sheltering by other buildings, which in turn induces large drag forces above the urban surface. It has been verified by the observational data collected in a densely built-up area of Beijing that, despite being substantially underestimated by other schemes, estimates from the SY08N scheme are comparable to the observationally obtained counterparts when at high building densities [24].

This paper aims to evaluate the role of the improved parameterization scheme in enhancing the simulation accuracy of urban atmospheric environment via WRF. After introducing the SY08N scheme into the UCM within WRF, two sets of numerical simulations are conducted, with parameterization scheme using the MGH98 scheme as the control experiment and that adopting the SY08N scheme as the sensitivity experiment. By comparing the model results with the observation, we evaluate the role of the SY08N scheme in improving WRF’s potentiality of representing the land surface-atmosphere exchange processes. The rest of this paper is organized as follows: Section 2 describes the coupled WRF-UCM model, the parameterization scheme, and the experiment design. Section 3 presents the model results and the mechanism behind results. Finally, we conclude this paper in Section 4.

2. Model and Experiment Design

2.1. The WRF-UCM Model

The mesoscale atmospheric model used in this study is the WRF model with Advanced Research WRF (ARW) dynamic core version 3.2 [25]. In this version, WRF has been coupled with the Noah land surface model and UCM. The UCM is a single-layer model developed for parameterizing the effects of urban canopy geometry such as the increased mechanical drag, turbulence production, and the trapping of radiation [26]. Such coupled WRF-UCM model has been applied to a dozen of metropolitan regions (e.g., Beijing, Guangzhou, Taiwan, Houston, New York City, and Salt Lake City), and its performance was well validated by observational data from multiple sources [2732].

In the framework of UCM, the considerable complexity of the urban land surface is reduced to a street canyon, where a road is bordered by two facing building walls. As shown in Figure 1, the urban surface is composed of horizontal (floor and roof) and vertical surfaces (wall). Buildings are located along identical roads, the length of which is hypothesized as far greater than their width. For the geometrical characteristics, all buildings are assumed to have the same width , so it is the same to all roads with a width of . Building heights among different street canyons are varied, and their average and variability are marked as and , respectively. Based on the four parameters, the plan area index (), the canyon (floor and wall) frontal area index (), and the roof frontal area index () which are required in parameterization schemes can be easily derived:

2.2. Parameterization of

As described in Section 1, the default parameterizations scheme used in the UCM within WRF model is the one proposed by Macdonald et al. [21] (MGH98 scheme). Another parameterization which will be used in this study is the one proposed by Cao et al. [24] (SY08N scheme), which is developed from Shao and Yang [22] scheme. Each parameterization scheme will be addressed in the following two subsections, respectively.

It is worth mentioning that the calculation of urban energy budget is split into two parts in UCM: one for the canyon (floor and wall) and another for the roof. In this context, there are in total three parameters in UCM reflecting the effects of : the canyon aerodynamic roughness length (), the canyon zero-plane displacement height (), and the roof aerodynamic roughness length .

2.2.1. The MGH98 Scheme

Macdonald et al. [21] proposed a semiempirical parameterization scheme intended for use in urban areas, which accounts for both the decline of at higher roughness densities and the drag differences caused by different obstacle layouts. It requires inputs of mean height, frontal area index, and plan area index of roughness elements: where is an empirical coefficient, is the drag coefficient and equals 1.2, is the Von Karman constant (usually taken to be 0.4), and is a correction factor for the drag coefficient. Macdonald et al. [21] have “calibrated” this parameterization using wind tunnel data and recommended that for staggered array of obstacles, and = 1.0 (also used here).

2.2.2. The SY08N Scheme

Based upon Shao and Yang [22], Cao et al. [24] proposed a new parameterization scheme for application in heterogeneous urban areas, which additionally takes into consideration the effects of building height variability upon variation. Concretely, the four geometrical parameters discussed in Section 2.1 (i.e., mean height, frontal area index, plan area index, and height variability of roughness elements) are used as input into the SY08N scheme: where is reference height (taken to be ), is the mean speed at the reference height, and is the friction velocity. According to Shao and Yang [22], is a weighted average of and , where is the value of roughness element drag coefficient at and is the value of surface drag coefficient at . The formulation of is written as where , , and represent the correction of roughness geometry to drag coefficient: where , , and are empirical constants which equal 3, 5, and 0.1, respectively. Following the definition of , is the value of surface drag coefficient at . Mathematically, and , where is the roughness length for bare surface (i.e., ) and is the roughness length for fully covered surface (i.e., ). Because building height variability in real cities remains to affect the flow field even when the underlying surface is totally occupied by buildings, the effects of height variability ought to be embodied in the expression for : where is an empirical constant (taken to be 193 here). Essentially, (3)–(6) are derived from the classical drag partition theory; that is, the total drag is partitioned into a pressure drag (), a ground-surface drag (), and a roughness-element surface drag (): where is a comprehensive demonstration of the influence of roughness element properties (e.g., aspect ratio) upon drag partition and thus . The mathematical formulation for is where is an empirical constant (taken to be 3.67 here). The calculation of is analogous to that of , except that in all equations is replaced with .

2.3. Experiment Design

The WRF-UCM model was run at 9,  3, and 1 km horizontal grid spacings with , dimensions of , , and grid cells, respectively. Figure 2 illustrates the nested domain configuration: the largest domain (D1) extends from 113°E to 120°E and from 37°N to 43°N, the second largest domain (D2) covers the greater Beijing area, and the innermost domain (D3) corresponds to the city center of Beijing. All of the domains use 41 layers in the vertical direction with 15 layers in the lowest 1 km. The setting of higher vertical resolution in the lower layer is primarily to improve the simulated accuracy of boundary-layer structures. The National Centers for Environmental Prediction (NCEP) operational Global Final (FNL) Analyses were used for the model initial and outermost lateral boundary conditions. The land use and land cover were characterized by the MODerate resolution Imaging Spectroradiometer (MODIS) for 2004. We use the Mellor-Yamada Janjic boundary layer scheme, which predicts turbulent kinetic energy and allows vertical mixing between individual layers within the boundary layer [33]. Other physical parameterization schemes include the Purdue Lin Microphysics scheme [34], Rapid Radiative Transfer Model longwave radiation scheme [35], the Goddard shortwave radiation scheme [36], Eta similarity theory [33], Noah land surface model [37], and an ensemble cumulus parameterization [38].

Two numerical experiments using different parameterization schemes are designed. One experiment was the control run using the MGH98 scheme (denoted as CTL run) and the other was the sensitivity run using the SY08N scheme (denoted as EXP run). Except for the selection of parameterization scheme over urban areas, the rest of the model settings were kept the same in both runs. Both simulations were initiated at 0000 UTC (i.e., 0800 LST) 27 February and ended 0000 UTC 1 March 2001. Hereafter, only model results from the last 24 hours (0000 LST 28 February-0000 LST 1 March) were analyzed.

Besides, as our main objective is to demonstrate the role of the improved parameterization scheme in improving WRF’s performance in urban atmospheric environment modeling, model results from the CTL run and the EXP run are validated against the observational data from BECAPEX [39, 40]. Here, data from two observing sites, Beijing 325 m Meteorological (B325) Tower and Fangzhuang which are located in the city center of Beijing (see Figure 2) are selected out for validation.

3. Model Results

Generally, there are two mechanisms for parameterization scheme affecting the simulation of urban atmospheric environment. One is the dynamic effect; that is, variation in causes variation of drag force which leads to variation in turbulence characteristics and momentum exchange. The other is the thermal effect; that is, variation in can also contribute to the variation of sensible and latent heat exchanges.

Figure 3 shows the three roughness parameters (, , and ) computed from the two numerical experiments. As can be seen, all the roughness parameters (, , and ) in the EXP run (0.626 m, 5.819 m, and 0.643 m) are larger than those (0.333 m, 5.718 m, and 0.126 m) in the CTL run, indicating that SY08N scheme with the inclusion of building height variability gives a higher estimation of than the MGH98 scheme. We will further investigate in detail how the enhanced roughness parameters would affect the diurnal variation of land surface-atmosphere turbulent exchanges and boundary-layer structures.

3.1. Friction Velocity

Friction velocity is a scaling variable for characterizing the near-surface friction stress (i.e., momentum flux). Figure 4 gives the comparison of simulated and observed diurnal variation of at the B325 Tower. As seen, the diurnal fluctuation of the observed is quite remarkable, varying from 0.2 m  to about 1.0 m . Two peak values of are found, with one peak value (0.91 m ) appearing near 1400 LST and another peak value (0.73 m ) at around 2300 LST. The enhanced suggests strong turbulent mixing occurred during those two periods. For the simulation, before 1400 LST, the variation of the observed with time is basically captured by the two numerical experiments, but the magnitude of the simulated is systematically lower than the observation. After 1400 LST, the performance of the WRF-UCM model becomes less satisfactory; the simulation results tend to lag behind the observation and the simulated magnitudes remain lower than the observation for most of the time.

The main differences between the results from the EXP run and the CTL run lie in the magnitude. As depicted in Figure 4, the differences between both runs are unanimously positive, inferring simulated by the CTL run is always weaker than that from the EXP run. The most significant positive difference is about 0.25 m , occurring at 2000 LST when simulated by the CTL run hits the valley of the day. Comparatively, the EXP run slightly improves the model performance for it gives a closer estimation of to the observation.

3.2. Sensible Heat Flux

For urban grids, sensible heat flux is approximately equal to the weighted average of sensible heat flux from roof (), wall (), and floor (), of which , , and are estimated separately. The influence of parameterization on is reflected primarily through its impacts upon and , which are small and less pronounced than that upon . Figure 5 presents the simulated and observed variation of with time at the B325 Tower. For the observed , significant diurnal variation can be found in the daytime, with values of ranging from near 0 W  at night to about 200 W  in the noon. Starting from 0800 LST, the observed begins to increase with time and reaches the peak at near 1300 LST, and then it drops quickly to almost 0 W  at around 1600 LST. Generally, both the magnitude and phase of the observed variation are well captured by the numerical experiments, although the simulation overestimates the observed sensible heat flux in the day and slightly underestimates the observation at the night.

Similar to the change of , the main differences between the simulated from the EXP run and the CTL run also exhibit in the magnitude, with remarkable difference appearing in the daytime and smaller difference in the nighttime. Overall, simulated by the EXP run is larger than that by the CTL run, with the most substantial difference (approximately −15 W ) identified at 1400 LST. The reason is primarily because the SY08N scheme in the EXP run gives higher roughness parameters (Figure 3) which leads to enhanced blocking effect of buildings. As a result, the resistance to heat transfer is strengthened and the sensible heat exchange is weakened.

3.3. Wind Speed Profile

At Fangzhuang, a tethered balloon for obtaining detailed meteorological soundings in the lower 1 km of the atmosphere was made every 3 hour to measure the wind, temperature, and humidity profiles. Figure 6 gives the comparison of the simulated and observed wind speed profiles at 0800 LST (i.e., daytime) and 2000 LST (i.e., nighttime). At 0800 LST (Figure 6(a)), the observed wind speed increases rapidly with height at the lower boundary layer and reaches the maximum value (about 10 m ) at around 450 m. The two simulation experiments reasonably reproduce the observed increase of wind speed with height and also the maximum wind speed near 450 m. However, the simulated magnitude of wind speed is much smaller, of which the simulated maximum wind speed at the height of 500 m is only about 4.5 m . The reason for model’s underestimation in velocity may be from the imperfection of PBL scheme [41], observation uncertainty [42], and mismatch of spatial scales between the model grid data and observations. As shown from Figure 6(b), the differences of the simulated wind speed between the EXP run and the CTL run are small, which are in the range of −0.4 m  and 0.2 m .

At 2000 LST (Figure 6(c)), the observed wind speed is found to increase rapidly with height in the lower level and reach the maximum value at the level of 200 m (i.e., the so-called low-level jet); then it gradually decreases with height and gets to the minimum value at about 500 m and again increases with height from 500 m above. The observed low-level jet is well captured by the EXP run at the height of around 100 m, though the magnitude is lower than the observation. However, the CTL run fails to reproduce the observed vertical wind profile, and it even simulates a minimum value at the height of 200 m. Besides, we find that the simulated wind speed in both experiments increases sharply with height of 400 m and above, which is nearly opposite to the observation. These biases may suggest that the WRF-UCM model needs to improve its performance in simulating the stable nocturnal boundary-layer structures in the future.

3.4. Temperature Profile

Figure 7 gives the comparison of the simulated and observed temperature profiles at 0200 LST (i.e., nighttime) and 1400 LST (i.e., daytime). At 0200 LST (Figure 7(a)), a suspended inversion layer is detected at the level of 100–200 m. Above 300 m, the observed temperature decreases monotonically with height at a rate of 0.86°C per 100 m. Such observed characteristics of temperature profile are all well reproduced by the two numerical experiments, except that the simulated temperatures below 500 m are slightly lower (approximately 2°C) than the observation. In agreement with the small differences between the EXP run and the CTL run at 0200 LST, temperature differences of both runs are also insignificant (within ±0.15°C) (Figure 7(b)).

At 1400 LST (Figure 7(c)), because of the strong daytime turbulent mixing which makes the atmosphere unstable, the observed temperature decreases monotonically with height at a rate of 1.07°C per 100 m. Likewise, the simulated temperature lapse rate is 1.11°C per 100 m, which agrees closely with the observation. Similar to what has been identified at 0200 LST, the differences between the results from the EXP run and the CTL run are also neglectable (within ±0.1°C) (Figure 7(d)).

3.5. Bulk Transfer Coefficient

Bulk transfer coefficients are key parameters for determining the transfer efficiency of momentum, heat, and moisture from the underlying surface to atmospheric reference height. In the WRF-UCM model, the expression for the bulk transfer coefficient of heat () is formulated as where is the 2 m friction velocity, is the roughness length for heat [43], is the universal function of heat [44], is the background surface roughness length (about 0.15 m), is an empirical constant (taken to be 0.1), Re* is the roughness Reynolds number, and is the kinematic molecular viscosity. Clearly, is not directly correlated with the roughness parameters for urban surfaces (, , and ).

For the expression of the bulk transfer coefficient of momentum (), where is the atmospheric reference height and is the universal function of momentum [44]. As can be seen, is directly tied to parameters such as and .

Figure 8 compares the differences of the simulated and between the EXP run and the CTL run at 1400 LST. Over rural areas (e.g., croplands; see Figure 2), differences between both experiments are nearly ignorable (within ±0.005), and it is the same with differences. However, differences of are fairly conspicuous over urban areas, although differences are still neglectable (within ±0.005). Averaged over the urbanized region, is 0.059 in the CTL run and increases to 0.081 (by a percentage of 34.2%) in the EXP run. Because the only difference between both runs rests with the roughness parameters over the urban areas rather than those over the rural areas, the simulated differences over urban areas are reasonable. To pull it a bit further, the large sensitivity of urban dynamic variables (friction velocity and wind speed) to parameterization can be partly connected attributed to the large sensitivity of to parameterization schemes for urban areas. Likewise, the small sensitivity of urban thermal variables (sensible heat flux and temperature) to parameterization is also partially attributable to the small sensitivity of to parameterization scheme over the urban land surface.

3.6. Turbulent Kinetic Energy

Turbulent kinetic energy (TKE) is a measure of turbulence intensity, which is closely related to the transport of momentum, heat, and moisture through the boundary layer. Figure 9 vividly shows the variation of the simulated TKE with height and its evolution with time at the B325 Tower. Basically, both the CTL run and the EXP run are consistent with each other in modeling the diurnal variation of TKE. The simulated evolution of TKE is as follows: before 0900 LST, TKE is very weak; at 1000 LST, TKE is gradually enhanced with the formation of turbulent mixing layer; at 1500 LST, the development of TKE reaches the pinnacle and stretch over the whole boundary layer which makes for the full mixing of momentum, heat, and moisture; at 1600 LST, TKE below 500 m continues to increase, with the maximum TKE value appearing near the surface; after 2100 LST, TKE is weakened again over the whole boundary layer.

On the whole, the EXP run simulates a stronger development of TKE than the CTL run. The differences between both runs become apparent after 1200 LST, with the most substantial positive difference (0.33 m2) emerging at around 1800 LST. As known, TKE is produced by shear, buoyancy, or through external forcing at low-frequency eddies scale. Usually, buoyancy predominates the development of TKE after sunrise (i.e., daytime), and shear plays the dominant role after sunset (i.e., nighttime). As learned from the previous analysis, parameterization affects mainly the dynamical variables rather than the thermal variables. Hence TKE differences should be apparent (or ignorable) when shear (or buoyancy) is dominating the development of TKE. Accordingly, TKE differences of both runs are prominent at nighttime while quite small in daytime.

4. Conclusions

By introducing a new parameterization (SY08N scheme) into the WRF-UCM model, which accounts for the combined effects of urban building morphology (height, frontal area index, plan area index, and height variability) upon variation, two sets of numerical experiments in the framework of the WRF-UCM model were conducted to investigate the impact of parameterization scheme upon atmospheric environment simulations over the urban area and its surrounding areas. The only difference between the two numerical experiments lies in the usage of parameterization for urban areas, of which the CTL run uses the default parameterization (MGH98 scheme) and the EXP run adopts the SY08N scheme. The evaluation and discussion focus on the simulation of friction velocity, sensible heat flux, wind speed profile, temperature profile, bulk transfer coefficient, and turbulent kinetic energy.

The comparison between simulated results shows that all the roughness parameters in the EXP run are larger than those in the CTL run, which indicates that, with the inclusion of building height variability, the SY08N scheme gives a higher estimation of than the MGH98 scheme. The enhanced roughness parameters in the EXP run result in strengthened drag and blocking effects exerted by buildings, which in turn lead to increased friction velocity, decreased wind speed in daytime, and boosted turbulent kinetic energy after sunset. However, as the influence of parameterization has not been incorporated in the computation of bulk heat transfer coefficient, the simulation of temperature and sensible heat flux are hardly sensitive to variations.

The comparison with observations demonstrates relatively better simulation of urban boundary-layer structures and land surface-atmosphere exchanges by the EXP run. For the boundary-layer temperature structures, both the EXP run and the CTL run reproduce the observed vertical profile and also the inversion layer at nighttime. For the boundary-layer wind structures, the EXP run reasonably reproduces the nocturnal low-level jet at around 200 m while the CTL run fails. Besides, the magnitude of land surface-atmosphere exchanges (e.g., friction velocity and sensible heat flux) as simulated by the EXP run is much closer to the observation. All these indicate that the SY08N scheme can improve the model performance in simulating the urban atmospheric environment.

Nevertheless, it is worth mentioning that we used the default settings for urban land use information in the CTL run and the EXP run. The urban land use information given by the default settings is “imaginary” and might not be representative of the Greater Beijing area. For example, roughness parameters were homogeneously distributed over the entire urban region in both runs (Figure 3). Because the verification of the SY08N scheme is also dependent on the input of “true” surface information, we plan to replace the default urban land use data with what was representative of the Greater Beijing area in future work and to further investigate the combined usage of the SY08N scheme and real urban land use data in improving the simulation accuracy of urban atmospheric environment.


:Aerodynamic roughness length (m)
:Building width (m), as illustrated in Figure 1
:Road width (m), as illustrated in Figure 1
:Average building height (m), as illustrated in Figure 1
:Building height variability (m), as illustrated in Figure 1
:Plan area index
:Frontal area index
:Frontal area index
:Canyon (floor and wall) aerodynamic roughness length (m)
:Canyon (floor and wall) zero-plane displacement height (m)
:Roof aerodynamic roughness length (m)
:Von Karman constant (—)
:Friction velocity
:Sensible heat flux
:Sensible heat flux from roof to the atmosphere
:Sensible heat flux from wall to the canyon space
:Sensible heat flux from floor to the canyon space
:Bulk transfer coefficient of heat (—)
:Mean wind speed
:Friction velocity at the height of 2 m 
:Roughness length for heat (m)
:Universal function of heat (—)
:Background surface roughness length (m)
:Empirical constant in the expression for (—)
Re*:Roughness Reynolds number (—)
Kinematic molecular viscosity
:Bulk transfer coefficient of momentum (—)
:Atmospheric reference height (m)
:Universal function of momentum (—)
:Drag coefficient (—) in the MGH98 scheme
:Empirical coefficient (—) in the MGH98 scheme
:Empirical coefficient (—) in the MGH98 scheme
:Reference height (m) in the SY08N scheme
:Mean speed at the reference height in the SY08N scheme
:Friction velocity in the SY08N scheme
:Roughness element drag coefficient at  (—) in the SY08N scheme
:Surface drag coefficient at   (—) in the SY08N scheme
:Surface drag coefficient at  (—) in the SY08N scheme
:Roughness length for bare surface in the SY08N scheme
:Roughness length for fully covered surface in the SY08N scheme
:Correction factor for drag coefficient (—) in the SY08N scheme
:Correction factor for drag coefficient (—) in the SY08N scheme
:Correction factor for drag coefficient (—) in the SY08N scheme
:Empirical constant (—) in the SY08N scheme
:Empirical constant (—) in the SY08N scheme
Empirical constant (—) in the SY08N scheme
Total drag in the SY08N scheme
:Pressure drag in the SY08N scheme
:Ground-surface drag in the SY08N scheme
:Roughness-element-surface skin drag (kg m ) in the SY08N scheme
Ratio of to   (—) in the SY08N scheme
:Empirical constant (—) in the SY08N scheme for
:Empirical constant (—) in the SY08N scheme for .

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This research is jointly supported by the “Strategic Priority Research Program” of the Chinese Academy of Sciences (Grant no. XDA05110000) and the National Natural Science Foundation of China (Grant no. 41305093). The authors thank Professor Zhiqiu Gao for providing the observational datasets for validation.