#### Abstract

The Coupled Ocean-Atmosphere-Wave-Sediment Transport model has been used to simulate Super Typhoon Yutu (2018). The impacts of four momentum transfer parameterization schemes (COARE, TY, OT, and DN) and three heat transfer parameterization schemes (COARE, GR, and ZK) on typhoon modelling have been studied by means of the track, intensity, and radial structure of typhoon. The results show that the track of Yutu is not sensitive to the choice of parameterization scheme, while the combinations of different parameterization schemes affect the intensity of Yutu. Among the four momentum flux parameterization schemes, three wave-state-based schemes (TY, OT, and DN) provide better intensity results than the wind-speed-based COARE scheme, but the differences between the three wave-state-based schemes are not obvious. Among the three heat flux parameterization schemes, the results of the GR scheme are slightly better than those of the COARE scheme, and both the GR and COARE schemes are significantly better than the ZK scheme, from which the intensity of Yutu is underpredicted obviously. The influence of the combination of different parameterization schemes on the intensity of the typhoon is also reflected in the radial structure of the typhoon, and the radial structure of the typhoon simulated by experiments with stronger typhoon intensity also develops faster. Differences of intensity between experiments are due mainly to the differences in sea surface heat flux, the enthalpy transferred from sea surface to the atmosphere has a significant impact on the bottom atmosphere wind field, and there is a strong correspondence between the distribution of enthalpy flux and the bottom wind field.

#### 1. Introduction

Tropical cyclones are one of the most severe weather phenomena, which cause huge loss of lives and property to coastal areas [1, 2]. Accurate predictions of the track and intensity of tropical cyclones, especially typhoons and hurricanes, are of great significance for reducing hazards of tropical cyclones and emergency management.

Due to the development of numerical forecasting models, typhoon track forecasts have been greatly improved during the past few decades. However, because of the complexity of the internal physical mechanism of typhoons and the difficulty of accurately describing some small-scale processes related to the intensification of typhoon in numerical forecasting models, the improvement of typhoon intensity forecasting is relatively slow [3].

As the energy source of typhoons, the ocean plays an important role in the generation and intensification of typhoons. An accurate description of the energy exchange between the atmosphere and the ocean is essential to the prediction of typhoon intensity [4–7]. In current models, the energy exchange between the atmosphere and the ocean is described by momentum flux and moist enthalpy flux, where the moisture enthalpy flux is the sum of sensible heat flux and latent heat flux. Momentum flux , sensible heat flux , and latent heat flux between air and sea can be expressed aswhere is the air density, is the wind speed at a height of 10 m, is the specific heat capacity of air, and is the latent heat of vaporization, , , and are momentum exchange coefficient (also called drag coefficient), sensible heat flux exchange coefficient, and latent heat flux exchange coefficient, respectively. In neural stratifications, , , and can be calculated fromwhere is the Von Kármán constant, is the reference height, often taken as 10 m, , , and are the roughness length for momentum, sensible heat, and water vapor, respectively.

For the calculation of momentum flux, early applications treated or as a function of wind speed. In low and moderate wind conditions (), many observations show that increases linearly with wind speed [8–10]. Therefore, the function of in low and moderate wind conditions can be expressed asFrom different observations, the coefficients and are fitted to different values (Table 1). It can be seen that the functions of fitted by different observations all reveal the monotonically increasing relationship of with , but the values of coefficients and differ significantly. The difference indicates that may not only depend on wind speed. Many studies have shown that wave state has an important influence on the surface roughness [11, 12].

Due to the lack of observations under high wind speed conditions, in some early models [21, 22], the linear relation between and is extrapolated to high wind speed conditions. However, some recent field experiments and laboratory experiments show that when the wind speed reaches a certain level, no longer increases with the increase of , but reaches a saturation value [23, 24], or decreases with the increase of [25]. This saturation effect is related to the wave breaking and the generation of sea foam under high wind speeds [26, 27]. Due to the important influence of wave state on momentum transfer, many wave-state-based sea surface roughness parameterization schemes have been proposed to calculate the momentum flux between air and sea, such as the wave-steepness-based parameterization scheme proposed by Taylor and Yelland [28] (equation (8)) and the wave-age-based scheme proposed by Drennan et al. [29] (equation (9)):where is the significant wave height, is the wave length at the spectral peak, denotes the wave steepness, is the wave speed at the spectral peak, is the friction velocity, and denotes the wave age. Compared with the scheme based on wind speed, schemes based on wave state describe the physical state of the sea surface more directly and are often used in the atmosphere-ocean-wave coupling model, in which the wave state is provided by the ocean model or the wave model.

The heat transported from ocean to atmosphere is the main energy source for typhoon [30, 31]. Due to the difficulties in measuring heat flux directly, the mechanism of air-sea heat transfer is relatively dubious. According to the method adopted, the parameterization schemes of heat flux can be roughly divided into the following categories: (1) based on the observations in low and moderate wind conditions, sensible heat flux exchange coefficient and latent heat flux exchange coefficient are treated as a constant [32]; (2) by fitting observations directly to dimensionless parameters, such as roughness Reynolds number , the sensible heat roughness and water vapor roughness can be denoted as a function of dimensionless parameters [33]; (3) describe the ratio between or and as a function of roughness Reynolds number and Prandtl number : [34]. As mentioned above, there are differences between different heat flux parameterization schemes, which will lead to differences in the typhoon modelling results.

This study aims to improve the understanding of the influence of air-sea energy exchange on typhoon modelling, especially the influence on intensity forecast. To examine the impact of momentum flux and heat flux on typhoon intensity, we evaluated the performance of four momentum flux parameterization schemes and three heat flux parameterization schemes in the simulation of Super Typhoon Yutu.

#### 2. Methods

##### 2.1. Momentum Flux Parameterization Schemes

In the numerical model, the momentum flux of the sea surface is determined by the drag coefficient (equation (1)), and can be calculated from (equation (4)). Therefore, different momentum flux parameterization schemes in the model are realized by different calculation schemes for .

###### 2.1.1. COARE

COARE (Coupled Ocean-Atmosphere Response Experiment) scheme was proposed based on the flux measurements collected from the TOGA (Tropical Ocean Global Atmosphere) COARE experiment [35], which is widely used in flux calculations in various numerical models [36, 37]; here we adopt the v3.5 version of COARE proposed by Edson et al. [10].

COARE divides into two parts:where represents the roughness corresponding to the part of the momentum transported entirely by the viscous effect of the sea surface when the sea surface is flat; in addition to the momentum transfer caused by the viscous effect, due to the actual unevenness of the sea surface, the wind produces horizontal wind pressure on the sea surface, which causes the horizontal momentum to be transferred to waves, and this part of the momentum transfer corresponds to roughness . is calculated fromwhere is the kinematic viscosity of air and is taken as 0.11 by COARE from observations. is often calculated from Charnock relation [38]:where is the Charnock coefficient and is the gravitational acceleration. Then, (10) is written as

Charnock coefficient is denoted as a function of wind speed by COARE based on observations:where m^{−1}s and .

###### 2.1.2. TY

Based on the observations from HEXOS (the Humidity Exchange over the Sea) [39], RASEX (the Risø Air-Sea Exchange) [40], and Lake Ontario [41] experiments, Taylor and Yelland [28] proposed a sea surface roughness parameterization scheme considering significant wave height and wave steepness :

Compared with the wind-speed-based COARE scheme, TY scheme describes the relationship between the physical properties of the sea surface and more directly.

###### 2.1.3. OT

From the momentum flux observations collected from ASGMAGE (an union of ASGASEX (Air Sea GAS EXchange) and MAGE (Marine Aerosol and Gas Exchange)) experiment, Oost et al. [42] proposed a scheme considering wave length and wave age :

Wave age represents the relative magnitude of wave speed and wind speed, and the ability of winds to transfer momentum to waves.

###### 2.1.4. DN

By filtering out momentum flux data from the pure wind wave field and deep water conditions from AGILE (measured from the 15-m research vessel *AGILE*) [43], FETCH (Flux, sea state, and remote sensing in conditions of variable fetch) [44], HEXOS, SWADE (Surface Waves Dynamics Experiment) [45], and WAVES (Water-Air Vertical Exchange Study) [46] experiments, Drennan et al. [29] proposed a scheme considering wave age :

Among the above four schemes, except for the COARE scheme, which uses wind speed to calculate , the other three schemes are all based on wave parameters. Since the atmospheric model generally does not directly provide wave parameters, TY, OT, and DN schemes are mainly used in the coupled model, in which the wave parameters can be provided by wave model.

##### 2.2. Heat Flux Parameterization Schemes

###### 2.2.1. COARE

In addition to the calculation scheme of , COARE also provides the calculation scheme of and . In COARE scheme, and are taken as the same value:where is the roughness Reynolds number.

###### 2.2.2. GR

GR scheme proposed by Garratt [47] calculate and from and ; values of and are slightly different:

###### 2.2.3. ZK

Based on the observations from TOGA COARE and SCOPE (San Clemente Ocean Probing Experiment) [48], Zilitinkevich et al. [34] proposed a scheme to calculate and in different sea surface roughness conditions:

= 0.1 is chosen as the threshold value because the fitting curves of observational data match the curves of smooth regime equations at that point.

##### 2.3. Evaluation Parameters

To quantitatively compare the results of the typhoon simulations, we introduced three parameters to evaluate the simulation results: the root mean square error (RMSE), Pearson correlation coefficient (*R*), and model skill (*S*), given aswhere and represent the observed and simulated values, respectively, with respect to time.

#### 3. Case Introduction and Experimental Design

##### 3.1. An Overview of Super Typhoon Yutu

Super Typhoon Yutu is the 26th named tropical cyclone in the Western North Pacific in 2018. At 1200 UTC 21 Oct 2018, Yutu formed as a tropical depression on the tropical ocean near 158°E and 8°N, it upgraded to a tropical storm at 0000 UTC 22 Oct, and then it reached the intensity of typhoon at 0000 UTC 23 Oct. At 2000 UTC 24 Oct, its intensity exceeded the Super Typhoon Mangkhut (22nd named tropical cyclone in the Northwest Pacific in 2018) and thus became the strongest tropical cyclone in 2018. Figure 1 shows the track and intensity changes of Yutu from the Best-Track data released by the JTWC (Joint Typhoon Warning Center).

##### 3.2. Experimental Design

Experiments in this study are conducted on the COAWST (Coupled Ocean–Atmosphere-Wave-Sediment Transport) [49] model. Super Typhoon Yutu is simulated in the atmospheric model, and wave model is activated to provide wave parameters to atmospheric model for the calculation of , because the activation of ocean model will cause a cold deviation of the sea surface temperature in the simulation of tropical cyclones over the Western North Pacific [50], and the ocean model was not activated in our experiments.

Weather Research and Forecasting (WRF‐ARW 4.1.5) model is the atmospheric model used in the coupled model. 6-hourly 0.25° × 0.25° GDAS (Global Data Assimilation System) Final Analysis data from NCEP (National Centers for Environmental Prediction) were used as the initial and boundary conditions. There are three two-way nested grids in the present work-D01, D02, and D03 (Figure 2). The outermost grid D01 has a horizontal spacing of 27 km and a time step of 90 s; the second grid D02 has a horizontal spacing of 9 km and a time step of 30 s; the innermost grid D03 is a vortex following grid moving along the typhoon center, the vortex is tracked at 500 hPa level, its position is calculated every 15 minutes, and the max vortex speed is taken as 40 m/s for the calculation of the new vortex center position; its initial position is drawn as the red box in Figure 2, the horizontal spacing of D03 is 3 km, and its time step is 10 s. A total of 44 vertical layers with a pressure top of 10 hPa are adopted. Purdue Lin’s [51] scheme is used as the microphysics scheme; both shortwave and longwave radiation are resolved by Rapid Radiative Transfer Model for Global Circulation Models (RRTMG) [52], Mellor-Yamada Nakanishi and Niino level 2.5 (MYNN) [53] is adopted as the planetary boundary layer scheme, surface layer model is MYNN, land surface model is set to Unified Noah [54], and Kain-Fritsch [55] cumulus scheme is only activated in D01 and D02.

Wave parameters needed (, , and ) in the calculation of are calculated by wave model SWAN (Simulating Waves Nearshore, 41.31) and are transferred to WRF through the coupler MCT (Model Coupling Toolkit); SWAN model also receives the 10 m wind field from WRF as the forcing field. Horizontal spacing of SWAN is about 9 km, and the time step is 180 s; the domain of SWAN is shown as the blue box in Figure 2.

Both atmospheric and wave model are initialized at 1200 UTC 22 Oct 2018 and are integrated for 5 days to 1200 UTC 27 Oct 2018. The simulation period covers the entire enhancement process of Yutu; in this period, track of Yutu is all over the sea, and the influence of different parameterization schemes on the momentum and heat transfer between ocean and typhoon can be fully analysed.

There are 12 experiments in this study: four momentum flux parameterization schemes and three heat flux parameterization schemes described in Section 2 are called pairwise. Except for the flux parameterization scheme, other settings are all the same. The parameterization schemes used in each group of experiments are summarized in Table 2.

#### 4. Results and Discussion

##### 4.1. Track and Intensity

Comparison between simulated tracks of Yutu and Best-Track data is plotted in Figure 3. Overall, the movement of typhoon from the southeast to the northwest is reproduced by every experiment, but the simulated moving speeds are slower than that of the Best-Track data. Differences of tracks between 12 experiments are not obvious, and the tracks of Yutu before 0000 UTC 25 Oct from 12 experiments are generally the same. In general, simulated tracks of Yutu are not sensitive to the flux parameterization schemes, which is in agreement with many typhoon simulation studies [3, 56], because typhoon track is mainly affected by the large-scale steering flow, while the small-scale flux transport has limited impact on it.

**(a)**

**(b)**

Comparisons between simulated intensities and Best-Track data are shown in Figure 4. Results of MSLP (minimum sea level pressure) are plotted in Figure 4(a). From the Best-Track data, we can see that Yutu intensified rapidly during 12–48 h of the simulation period, and its MSLP reduced by about 80 hPa, but this rapid intensification process has not been reproduced well by experiments; the reduction speeds of simulated MSLP are slower than that of Best-Track data. Although the mechanism of the rapid intensification of tropical cyclone is still controversial [57], it is widely accepted that the thermodynamic and kinematic properties of the TC core play an important role in rapid intensification [58, 59]. In the numerical simulations, how the model calculates convective process matters a lot for the rapid intensification of TC. Hence, the failure in reproducing the rapid intensification of Yutu may be due to the inability of the current convection calculation scheme to reproduce the convection processes in the TC core. Results simulated by different experiments differ a lot; MSLP from four experiments using ZK scheme as the heat flux parameterization scheme (CR_ZK, TY_ZK, OT_ZK, and DN_ZK) is much larger than others, which denotes the typhoon intensities are affected by the adoptions of ZK scheme. The performances of COARE, GR, and ZK schemes in numerical simulations have also been evaluated by Bao et al. [4] using the National Atmospheric and Oceanic Administration/Environmental Technology Laboratory (NOAA/ETL) regional air-sea coupled modelling system, from their results (Figure 1 in Bao et al. [4]). ZK scheme showed the worst performances, hurricane simulated by ZK scheme did not intensify at all, and GR scheme was slightly better than COARE scheme. Due to the different MSLP reduction speeds simulated by different experiments, the simulation results of MSLP showed a big difference at 120 h: the smallest is about 900 hPa (DN_CR) and the largest is about 940 hPa (CR_ZK); the difference exceeds 40 hPa.

**(a)**

**(b)**

Figure 4(b) presents the results of UMAX (maximum wind speed at 10 m level); the increasing trend of UMAX corresponds to the decreasing trend of MSLP, the rapid increasing of UMAX during 12–48 h is also not reproduced well by experiments, and UMAX simulated by the four experiments using ZK scheme as the heat flux parameterization scheme also showed obvious weak deviations.

Listed in Table 3 are the RMSE, Pearson correlation coefficient, and model skill of simulated MSLP from 12 experiments. Among the 12 experiments, RMSE of TY_CR is the smallest (17.30 hPa), and the largest is that of CR_ZK (34.49 hPa). The difference in Pearson correlation coefficient between different experiments is not obvious. The model skill also shows significant differences: the highest is TY_CR (0.8447) and the lowest is CR_ZK (0.4198), which is consistent with the results of RMSE.

To compare the performances of different parameterization schemes on MSLP simulation more intuitively, Table 4 lists the average RMSE, R, and S of experiments that use the same schemes. From the results of four momentum flux parameterization schemes, it can be seen that the performance of wave –state-based TY, OT, and DN schemes is generally better than the wind-speed-based COARE scheme. Compared with the COARE scheme, RMSE of MSLP simulated by three wave-state-based momentum flux parameterization schemes is reduced by about 3 hPa on average. COARE is one of the best wind-speed-based parameterization schemes, which has been widely used in momentum flux calculations [60, 61]; the better results from wave-state-based schemes reveal that the characteristics of momentum flux at typhoon condition are better captured by the wave-state-based schemes; this denotes that considering the impacts of wave state can provide more information for the parameterization of . Similar conclusions are presented by Drennan et al. [62] and Prakash et al. [63]. From the results of three heat flux parameterization schemes, it is demonstrated that the MSLP results simulated by the ZK scheme are worse than the other two schemes. The MSLP RMSE is significantly higher, and the model skill is significantly lower. Results of the GR scheme are slightly better than the COARE scheme. The relative performances of three heat flux parameterization schemes are consistent with Bao et al. [4].

RMSE, Pearson correlation coefficient, and model skill of simulated UMAX from 12 experiments are given in Table 5. Among 12 experiments, RMSE of TY_CR is the smallest (11.47 m/s), and the largest is that of CR_ZK (21.35 m/s). Results of model skill are similar to RMSE; the highest is TY_CR (0.7725), and the lowest is CR_ZK (0.3624). Table 6 lists the group averaged root mean square error (RMSE), Pearson correlation coefficient (*R*), and model skill (*S*). From the results of four momentum flux parameterization schemes, it can be seen that the performance of wave-state-based schemes (TY, OT, and DN) is generally better than the wind-speed-based scheme (COARE), which are consistent with the results in Table 4. RMSE of UMAX simulated by TY, OT, and DN is 1.2 m/s smaller than COARE on average. By comparing the results of COARE, GR, and ZK, it is demonstrated that the UMAX results simulated by the ZK scheme are worse than the other two schemes; results of the GR scheme are slightly better than the COARE scheme.

##### 4.2. Time Evolution of Radial Structure

To analyse the time evolution of the typhoon radial structure, the Hovmöller diagrams of azimuthally averaged SLP (sea level pressure) and are plotted in Figures 5 and 6, respectively. The intensification and the structure development of typhoon are presented in Figure 5; results in Figure 5 are consistent with the results in section 4.1, and experiments using ZK as the heat flux parameterization scheme tend to underestimate the intensity of Yutu. The faster the typhoon intensifies, the faster the radial pressure gradient increases. Here we take CR_ZK (Figure 5(c)) and TY_CR (Figure 5(d)) as examples. Typhoon simulated by TY_CR intensifies rapidly; it takes 42 h for TY_CR to intensify to 950 hPa, while typhoon simulated by CR_ZK intensifies slower, and it takes about 75 h for CR_ZK to intensify to 950 hPa; in addition to the difference in the speed of the central pressure drop, the development speed of their radial structure is also significantly different. At 72 h, the 980 hPa isobar of CR_ZK is located about 70 km from the center of typhoon, while the 980 hPa isobar of TY_CR is located about 90 km from the center of typhoon, indicating that the radial structure simulated by TY_CR is stronger. It is worth mentioning that typhoon structure simulated by the four experiments that use the ZK heat flux parameterization scheme (CR_ZK (Figure 5(c)), TY_ZK (Figure 5(f)), OT_ZK (Figure 5(j)), and DN_ZK (Figure 5(l)) is not sufficiently developed, which is consistent with the results in Section 4.1.

Figure 6 shows the time evolution of azimuthally averaged tangential winds and radial winds. It can be seen that the wind speeds are relatively small near the center and reaches the maximum value as radius increases and then decreases with radius. The values of tangential winds can reach to more than 50 m/s, while the maximum value of radial winds is only 25–30 m/s, indicating that the bottom wind field structure of a typhoon is dominated by tangential rotation, accompanied by a weaker convergence effect. The maximum value of radial wind generally appears at a distance of 60 km from the center, while the maximum value of tangential wind is probably distributed between 40 and 60 km, indicating that the region with strongest convergence effect is slightly outside the region with the strongest rotation effect. This feature has also been presented by Green and Zhang [6] (see Figure 7 therein). From the results of experiments with strong typhoon intensity (CR_GR (Figure 6(b)), TY_CR (Figure 6(d)), TY_GR (Figure 6(e)), OT_CR (Figure 6(g)), and DN_CR (Figure 6(j)), the maximum tangential wind speed radius at 120 h is roughly 40–50 km, while for experiments with weaker typhoon intensity (CR_ZK (Figure 6(c)) and DN_ZK (Figure 6(l)), the maximum tangential wind speed radius at 120 h is about 60 km, indicating that not only the horizontal scale of typhoon, but also the radial distributions of winds are affected by momentum and heat transfer between air and sea.

##### 4.3. Heat Flux at Sea Surface

To determine the causes of differences in typhoon intensity from different experiments, we analysed the heat flux at sea surface simulated by different experiments. Time evolution of azimuthally averaged latent heat flux and sensible heat flux at sea surface is plotted in Figures 7 and 8, respectively. By comparison, heat transported by latent heat flux is stronger than the sensible heat flux. Similar to the distribution of wind speed in Figure 6, the distributions of two heat fluxes both show the characteristics of increasing at first, and then decreasing along the radial direction. Heat transfer in Figures 7 and 8 corresponds to the development of wind field in Figure 6, and typhoon simulated by experiments with stronger heat flux can acquire more energy from the ocean, which lead to the development of wind field. The poor results from ZK scheme in the intensity simulation are caused by the poor results in the calculation of heat flux, which are caused by the unreasonable calculation of and . The enthalpy exchange coefficient calculated by ZK decreases with increasing wind speed (cf Figure 2(b) in Bao et al. [4]), which is not consistent with the consensus that enthalpy exchange coefficient increases or keeps constant with increasing wind speed [64–66].

To verify this result, we plotted the distribution of enthalpy flux (sum of latent heat flux and sensible heat flux) and 10 m level wind speed at 120 h (1200 UTC 27 Oct 2018) in Figure 9. It is demonstrated that there is a strong correlation between distributions of wind speed and distributions of enthalpy flux; the contours of wind speed are parallel to the contours of enthalpy flux in general, and the distribution of horizontal gradient of enthalpy flux is consistent with that of wind speed. Based on the strong correlation between enthalpy flux and wind speed, it can be concluded that different heat flux parameterization schemes lead to differences in sea surface heat flux, which in turn has an important impact on the intensity of typhoons.

#### 5. Conclusions

It is of great significance to accurately forecast the intensity of tropical cyclones, especially typhoons. Energy transfer between air and sea is crucial for the evolution of typhoon intensity. In this paper, we use the coupling model COAWST to study the impacts of four momentum flux parameterization schemes and three heat flux parameterization schemes on the simulation of Super Typhoon Yutu. From the results, we concluded the following.(1)Track simulations of Yutu are not sensitive to the flux parameterization schemes. Differences between schemes are mainly presented in the intensity simulations. Among the four momentum flux parameterization schemes, three wave-state-based schemes (TY, OT, and DN) provide better intensity results than the wind-speed-based COARE scheme. This is caused by the different performances of them in parameterization of . Three wave-state-based schemes (TY, OT, and DN) provide better intensity results than the wind-speed-based COARE scheme because wave state has a nonnegligible impact on the sea surface roughness.(2)Among the three heat flux parameterization schemes, the results of the GR scheme are slightly better than that of the COARE scheme; both the GR and COARE schemes are significantly better than the ZK scheme, from which the intensity of Yutu is under predicted obviously. The poor result of ZK is because it calculates an unreasonable and . The enthalpy exchange coefficient calculated by ZK decreases with increasing wind speed, which is not consistent with the consensus that enthalpy exchange coefficient increases or keeps constant with increasing wind speed.(3)Differences of intensity between different schemes are related to the simulated sea surface heat flux. Heat transferred at sea surface has a significant impact on the wind field. There is a strong correlation between the distributions of wind speed and distributions of enthalpy flux; area with large enthalpy flux can acquire more energy from the ocean, which leads to the development of wind field.

#### Data Availability

The JTWC TC Best-Track data used to support the findings of this study have been deposited in the Joint Typhoon Warning Center website (https://www.metoc.navy.mil/jtwc/jtwc.html?western-pacific). The NCEP analysis data are obtained from https://rda.ucar.edu/datasets/ds083.3/.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the National Key R&D Program of China (Grant no. 2018YFB0203801) and the National Nature Science Foundation of China (Grant no. 41605070).