Characteristics of land-atmosphere interaction parameters in hinterland of the Taklimakan Desert

The importance of the energy exchange between the land surface and the atmosphere can be characterized by bulk transfer coefficients for momentum, Cd, and heat, Ch. The diurnal and monthly variations of both bulk transfer coefficients and lengths of surface roughness are analyzed. Based on observed data from January to December 2009 in hinterland of the Taklimakan Desert, the characteristics of aerodynamic roughness length, z0m, and thermal roughness length, z0h, are discussed. It should be noted that the diurnal and monthly variations of the parameters are fundamentally different from those reported in vegetated areas. Specifically, four unique features can be identified in the surface layer. First, in Taklimakan Desert, z0m does not vary with seasons; however, it significantly depends on wind speed. Second, z0h is higher in the daytime and lower at night, showing obvious diurnal characteristics. The high values appear at sunrise and sunset. Third, both Cd and Ch have two peaks, one peak at sunrise, and another one at noon. Fourth, both Cd and Ch have larger values in winter season and smaller values in summer season.

turbulent momentum and energy transported from the surface of land to the atmosphere. They directly reflect the land-atmosphere coupling strength. Thus, investigating their characteristics and variations is essential for calculating the energy exchange of land-atmosphere interactions.
Taklimakan Desert is the driest area in China. It is also the second largest flow desert in the world, following the Sahara desert in Africa. Therefore, investigating C d and C h in the hinterland of Taklimakan Desert can reveal the physical processes governing the regional climates. The characteristic variations of C d and C h in the northern edge of TD have been investigated. However, in hinterland of the TD, C d and C h exhibit significant differences from other regions. This study defines C d and C h , using data directly measured by the atmospheric environment observation station, located in Tazhong in the hinterland of the TD (Fig. 1, hereinafter Tazhong Station).
The defined coefficients can be applied to similar land surfaces where direct measurements are not possible or not practical. Prior to the process, two key parameters, i.e., aerodynamic roughness length (z 0m ) and thermal roughness length (z 0h ) should be determined. They are crucial parameters for calculating the turbulent flux via bulk transfer equations. The aerodynamic roughness length and thermal roughness length are defined as the surface nonslip condition and surface temperature, respectively, which can be applied in the framework of the Monin-Obukhov Similarity (MOS) theory. Based on the MOS theory, z 0m is the height at which the speed of the extrapolated wind vanishes, and z 0h is the height at which the temperature of the extrapolated air is equal to the surface temperature. Generally speaking, z 0h is significantly different from z 0m because the momentum transport is partly dependent on the turbulent resistance on roughness obstacles, whereas the heat transport is not related to these obstacles. z 0h refers to the parameterization of heat transport mechanisms near the surface, where the molecular viscosity and the molecular thermal diffusivity of air may have significant impacts on the heat transport 26 . The results from both experiments and theoretical analyses have demonstrated that z 0m ≠ z 0h in many cases 25,27 . Many published works 23,[28][29][30][31][32][33] have shown that z 0h is not constant. Instead, it has been demonstrated to have great diurnal and monthly variations.

Results and Discussions
Determination of aerodynamic roughness length. The continuously obtained data, which was recorded at 30 min intervals, were used to estimate z 0m in 2009. Due to the observation errors in the data or the fact that the meteorological conditions did not meet the similarity theory, the resulting ln(z 0m ) may be multiple values. Therefore, the distribution of ln(z 0m ) was calculated by SAS (Statistics Analysis System). The probability distribution function (PDF) of ln(z 0m ) was obtained. The value of z 0m at the peak frequency of the histogram of ln(z 0m ) was considered as the optimal value. Figure 2 shows at the peak frequency, ln(z 0m ) is −5.77. Thus, ln(z 0m ) = −5.77 and z 0m = 3.11 × 10 −3 m.
In general, z 0m is a function of the vegetation that changes with the season. In addition, z 0m can also be affected by the speed of wind, as discussed in the following section 'Surface roughness lengths' . However, the "seasonally varying vegetation" plays no role when analyzing a bare desert surface. Therefore, according to Eq. (2), the variation of z 0m is thus primarily related to u/u * and z/L. Further analyses shown (Fig. 3) that the u * and z/L are closely related to the wind speed. Liner regression equation, u* = 0.05 u + 0.026, indicates u * linear dependence with u, and z/L tends to zero with wind speed increased. So, in the hinterland of TD, the Eq. (2) can be identified as: . Hence, z 0m is primarily relate to wind speed. As indicated by the results in Fig. 4a,b, ln(z 0m ) is the reciprocal of u, but the integrated stability correction term for wind results in a time lag between the peaks of ln(z 0m ) and wind speed, on both diurnal and seasonal time scale. Thus, during the warm period from April 1 to September 30, u is stronger and z 0m is lower. During the cold period from October 1 to March 3, u is weaker and z 0m is higher. In addition, it varies monthly, from 1.95 × 10 −3 m to 9.68 × 10 −3 m; the variation is much smaller than those over vegetated regions. The monthly changes in percentile are shown relative to the annual mean 3.74 × 10 −3 (Fig. 4d). Hence, the scatter plot of daily mean ln(z 0m ) and u demonstrates that ln(z 0m ) is inversely correlated with u (Fig. 4c). Their relationship can be fitted as: Determination of thermal roughness length. Compare with z 0m , z 0h displays more significant diurnal variations (Fig. 5a). The mean value of ln(z 0h ) in the daytime is considerably higher than at night, and the daytime variations are larger than those at night. z 0h exhibits diurnal variation, which is 9.78 × 10 −7 m during the night time and 3.28 × 10 −4 m in the daytime. Two peaks are observed for ln(z 0h ), i.e., one at sunrise and the other one at sunset. There is a slight trough between the two peaks. The two peaks occur at times when the atmospheric stability is near-neutral. In other words, the atmospheric stability will change from stable to unstable at sunrise, and change inversely at sunset. Under near-neutral conditions, the second term, -(θ a -θ 0 )/θ * , i.e., in Eq. (3), is close to peak (Fig. 6), while the third term Ψ h (z/L), which is the integral stability correction term of temperature profile, is close to zero. Therefore, ln(z 0h ) reaches maximum values at sunrise and sunset. In addition, ln(z 0h ) exhibits a greater monthly variation (Fig. 5b), with a peak occurring in July and a valley in January. It varies from 2.52 × 10 −6 m to 2.31 × 10 −5 m, with a mean of 7.93 × 10 −6 m. Therefore, z 0h has a larger value in summer and a smaller value in winter. Since the variation of z 0h in the daytime is considerably higher than at night, and the daytime in summer is longer than in winter, accordingly, which result in the peak values occurred in July and valley in January. (7)   (2) using all data from 2009. The value of ln(z 0m ) was calculated to be −5.77 at the peak frequency in the curve, which was considered as the optimal value. The solid curve is fit for the histogram. www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ station than in the north edge of the TD 20 . The mean value of C d was larger than in Gobi and other desert areas (e.g., Dunhung and HEIFE) 9,34,35 , but C h similar to those areas.

Determination of bulk transfer coefficients. Equation
As is evident in Fig. 7, C d and C h are considerably higher during the daytime than the night. When z/L at the station changed from −0.94 to 0.36, C d varied from 2.95 × 10 −3 to 5.81 × 10 −3 and C h varied from 1.24 × 10 −3 to 2.56 × 10 −3 . During the daytime, C d had the mean value of 5.05 × 10 −3 and C h had the mean value of 2.40 × 10 −3 . During the nighttime, the mean values of C d and C h were 3.06 × 10 −3 and 1.31 × 10 −3 , respectively. Both C d and C h are inversely correlated with the atmospheric stability, z/L (Fig. 7c). This result is consistent with the results reported in previous studies for mid-latitudinal areas 25,36,37 . But their peaks do not match closely to the atmospheric stability, because the values are also influenced by surface aerodynamic and thermal roughness lengths. In general, as the unstable surface layer (z/L < 0) increase, the bulk transfer coefficients will also increase, while as the stable surface layer (z/L > 0) increase, the bulk transfer coefficients will decrease. Eq. (10) indicates clearly that the bulk transfer coefficients are relate closely to atmospheric stability and surface roughness lengths. According to the significant diurnal variations of roughness lengths (Fig. 5a), two peak values appear at sunrise and sunset, and z/L changes dramatically from stable (z/L > 0 at 6:00) to unstable (z/L < 0 at 7:00) (Fig. 7c, Fig. 8c), and Ψ m (z/L) and Ψ h (z/L) exhibit the largest rates of variation, as their values increase rapidly from negative to positive (Fig. 10c). As a result, C d and C h also showed diurnal changing patterns, and they responded to changes in the surface roughness lengths and atmospheric stability. Thus, C d has one peak at sunrise and another highest peak at noon, because the largest value of z 0m appears at sunrise, but the variation of z 0m is small, the highest peak dominated by z/L. Similar to C d , C h also has two peaks at sunrise and at noon, respectively, but the highest peak appears at sunrise due to the z 0h has dramatically highest values at sunrise, thus, z 0h dominates the highest peak, and z/L determines another peak at noon. www.nature.com/scientificreports www.nature.com/scientificreports/ We also examined the seasonal variations of C d and C h , and their relationships with the atmospheric stability. Figure 8 shows that the aforementioned "sunrise effect" is more significant in the cold season than the warm season, because the surface layer at sunrise in the cold season is more unstable than the warm season (Fig. 8c).
As we known, increased wind speed enhances atmospheric turbulence intensity, consequently decreasing the unstable conditions. The wind speed in the cold season is lower than in the warm season (Fig. 11), which results in increased instability in the cold season.
We note that the "sunrise effect" may be unique for deserts because the works from vegetated areas reveal different temporal patterns 28,37,38 . Both C d and C h values in vegetated areas start to rise at sunrise and gradually reach peaks at noon. Another unique feature is that both C d and C h values have large values in winter and smaller values in summer because of inverse atmospheric stability (Fig. 9). This phenomenon is nearly completely opposite to that of vegetated areas in southern China 37 , where there are little dust storms happen. C d varied from 4.08 × 10 −3 in April to 6.16 × 10 −3 in January, and C h varied from 1.99 × 10 −3 in April to 2.70 × 10 −3 in January when z/L ranged from −1.05 in January to −0.16 in April. The variation patterns exhibit extremes in April because the amount of dusts and the number of dust storm days are higher than any other month. This is due to the dramatic interaction of cold and warm air, which enhances the atmospheric turbulence intensity, and decreases the unstable conditions. In addition, if the atmospheric stability is near-neutral, from Eq. (10), z/L = 0 and Ψ m (z/L) = Ψ h (z/L) = 0. Therefore, the surface roughness lengths only determine the bulk transfer coefficients in the near-neutral conditions. Thus, the means of the bulk transfer coefficients were obtained under different conditions (Table 1).

conclusions
Momentum fluxes and sensible heat are two key parameters that characterize the energy exchange process between the land surface and the atmosphere. These parameters can be obtained by calculating the bulk transfer coefficients for momentum (C d ) and heat (C h ) using the bulk transfer equations. Two important parameters, i.e., the aerodynamic roughness length (z 0m ) and the thermal roughness length (z 0h ), need to be determined in the calculation process. Based on the analysis results from the observed data, the following conclusions can be obtained:  www.nature.com/scientificreports www.nature.com/scientificreports/ (1) The aerodynamic roughness length, z 0m , in the Taklimakan Desert does not have seasonal variation, instead it is strongly correlated with wind speed. On the contrast, in the vegetated areas, z 0m is greatly dependent on the seasonal variations of vegetation. Therefore, a parameterization scheme of z 0m , as it relates to wind speed, is obtained. z 0m decreases at the warm temperature and increases at the cold temperature, with an highest value of 3.11 × 10 −3 m. (2) The thermal roughness length, z 0h , in the Taklimakan Desert is considerably higher during the daytime than at night. On the contrast, in the vegetated areas, z 0h is higher during the night than the day. Two peaks of z 0h , i.e., at sunrise and sunset, were observed. These peaks occur when the atmospheric stability reaches near-neutral conditions. The mean value of z 0h is 7.93 × 10 −6 m. (3) The bulk transfer coefficient values (C d and C h ) have two peaks, one peak appears at sunrise (i.e., so-called "sunrise effect"), due to the instability occurring at sunrise, and the largest roughness lengths for aerodynamic and thermal occurring at sunrise, another peak appears at noon. The "sunrise effect" of the Taklimakan Desert contrasts with the "noon effect" found in vegetated areas, where the values gradually reach their peaks at noon. C d and C h had the mean values of 4.68 × 10 −3 and 2.26 × 10 −3 , respectively. (4) Another unique feature is that both C d and C h have larger values in winter and smaller values in summer, being nearly completely opposite to those of vegetated areas. C d varied from 4.08 × 10 −3 in April to 6.16 × 10 −3 in January, and C h varied from 1.99 × 10 −3 in April to 2.70 × 10 −3 in January. In addition, the aforementioned "sunrise effect" is much stronger in the cold season than the warm season, due to the more unstable surface layer at sunrise in the cold season.
Some limitations of this study include that the results are only derived from the Taklimakan Desert. We are not sure whether the characteristics of these parameters agree well with other deserts. However, the topic will be the focus of future studies.

Data and Methods tazhong station. The Taklimakan Desert is located in the Tarim Depression surrounded by mountains,
which is a large internal drainage. The Tianshan Mountains and the inland Tarim River from east to west are located in the north. The Pamir Plateau is located in the west and the Kunlun Mountains are located in the south. The east side of Tarim Depression opens into the Lop Nor Depression. The Taklimakan Desert covers an area of 337,000 km 2 and is the largest desert in China. It also has the greatest aridity in China. Shifting dunes cover 85% of the total desert, with dunes reaching 20-200 m in height. From the observed data between 1996 and 2013 in the Tazhong Station (Fig. 1a), the annual average temperature is 12.4 °C, the annual average precipitation is only 23.0 mm, and the annual average potential evaporation is 3800.0 mm. The highest temperature in the record is 45.6 °C while the lowest temperature in the record is −32.7 °C. In the summer, the surface temperature can be as high as 80.0 °C. The wind speed has the average value of 2.5 m s −1 and with the highest instantaneous value of as high as 24.0 m s −1 . There are 260 days with sand and dust weather annually.
The data in this paper were collected in the Tazhong Station, which is the only site in the desert hinterland away from the surrounding cities of the Tarim Depression. It is at a linear distance of approximately 300 km from Luntai in the north, 550 km from Sache in the west, 220 km from Niya County in the south and 380 km from Ruoqiang in the east. This unique environment provides good conditions to investigate the atmospheric boundary layer in a desert. Tazhong Station is comprised of a main station (38°58′05″N, 83°39′35″E, altitude 1104 m) and a subsidiary station (38°58′52″N, 83°38′29″E, altitude 1103 m). The main station, located at the oilfield operation area, includes an 80 m gradient detection tower, a three-layer eddy-covariance (EC) system and a radiation observation system. The subsidiary station is located 2.2 km northwest of the main station. The systems in the subsidiary station include an EC system, radiation sensors, soil heat fluxes sensors, and an automatic weather station (AWS). The subsidiary station is located in an open environment with a relatively flat underlying surface in a natural quick sand area. The prevailing wind direction is northeast. The complex sand dunes are present around the subsidiary station, including 850 m in the east, 1600 m in the south and 1700 m in the west. All sand dunes are below 40 m, and run in a northeast-southwest direction, in agreement with the prevailing wind direction. Our study site in this paper is based on the subsidiary station.
instrumentation. The main instrumentations at the study site include an eddy-covariance (EC) measurement system, a radiation observation system (Fig. 1b) and an automatic weather station (AWS). In the EC system, a 3D sonic anemometer (CSAT3, Campbell Scientific Inc., USA) is used to measure the velocity and sonic virtual temperature in three dimensions. The EC system is installed 3 m above the ground. The raw data were collected by a CR5000 data logger (Campbell Scientific Inc., USA) at a sampling rate of 10 Hz. The radiation observation system is comprised of four components (CNR-1, Kippp & Zonen, The Netherlands) and used to measure solar radiation and far infrared radiation, which are the downward and upward shortwave and the downward and www.nature.com/scientificreports www.nature.com/scientificreports/ upward longwave radiation fluxes, respectively. These components were installed at a height of 1.5 m on the same mast as the EC system. The raw radiation data were stored by a CR1000 data logger (Campbell Scientific Inc., USA) at 1 Hz. The AWS, located 30 m northeast of the EC system, was used to measure wind speed/direction profile, air temperature/humidity profile, air pressure, and surface infrared temperature. The AWS data were detected at the sampling rate of 0.1 Hz, and stored by the CR1000 data loggers at 1-minute interval. The instruments are described in details in Table 2. All instrumentations were powered by solar panels and batteries. Raw data were stored on CF cards and output to the laboratory for post-processing every month. The data from January 1 to December 25, 2009 were strictly processed and the average processing time was 30 minutes.
Data processing. The post-processing software EdiRe (University of Edinburgh, http://www.geos.ed.ac. uk/abs/research/micromet/EdiRe) was used to acquire the raw data at 10 Hz. Post-processing included spike removal, correction of sonic virtual temperature, the performance estimation of the planar fit coordinate rotation, corrections of density fluctuation (WPL-correction) 39 and correction of frequency response 40 . In addition, quality control of the half-hour flux data 41 was conducted based on the following criteria: (1) Data during sensor malfunction (e.g., when there was a faulty diagnostic signal) were rejected; (2) Data within 1 hour window of precipitation were rejected; (3) When the missing data constituted more than 3% of the 30 minutes raw recording, the incomplete 30 minutes data were rejected; (4) Data recorded at the friction velocity below 0.01 m s −1 during the nighttime were rejected 42 ; and (5) When the wind speed was below 1.0 m s −1 , and the sensible heat flux had a value below 10 W m −2 or the opposite sign to the temperature difference between surface and air (surface minus air), data were rejected. The average value of 30-minute interval radiations was also calculated by removing the out-of-range data. AWS data beyond physical possibilities were rejected. In fact, the AWS data were not used other than as a reference.

Theory and Methodology
Surface roughness lengths. Neither z 0m nor z 0h has clear physical meaning. Thus, their values cannot be directly measured. Instead, two methods, i.e., the flux method and profile method, are typically used to estimate z 0m and z 0h from the observed data. In the flux method, the lengths are obtained by applying the bulk transfer equations at a level of observation using eddy-correlation fluxes with specific stability functions. In the profile method, the lengths are obtained from the observed wind and temperature profiles. Both the flux method and the profile method are based on the MOS theory. The integral gradient of the wind and temperature profiles in a horizontally homogeneous surface layer with a stability correction function can be expressed as: where ⁎ u (m s −1 ), u (m s −1 ), θ ⁎ (K) and θ a (K) refer to the observed frictional velocity, wind speed, temperature scale and air potential temperature at height z (m), respectively. They are known and obtained using the EC system. K (with the value of 0.4) is the von Kármán constant, θ 0 (K) is surface temperature, L (m) is the Obukhov length, z/L is the dimensionless stability parameter (A negative value means unstable while a positive value means stable) and Pr is the turbulent Prandtl number (Pr is 1 when z/L > 0 and 0.95 when z/L < 0). Pr indicates the ratio of the eddy diffusivities between momentum K m and heat K h , i.e., Pr = K m /K h . Because the surface of the Taklimakan Desert is bare, sandy, fairly homogenous and smooth, the displacement height, d, is thus negligible, i.e., d = 0.
Ψ m (z/L) and Ψ h (z/L) are the integrated stability functions for the wind and temperature profiles, respectively. Their formulae were first established based on the data collected from Kansas wheat-farming land by Businger et al. 43 . Afterwards, many studies [44][45][46][47][48][49][50][51] have conducted extensive research on stability functions. It is unknown if those formulae remain valid for the arid surface of the Taklimakan Desert. Thus, the verification and evaluation of empirical stability functions was investigated using the observed data of both flux and gradient. Appropriate Ψ m (z/L) and Ψ h (z/L) were determined for stable (z/L > 0) and unstable (z/L < 0) surface layers, as given by: h where x = (1-12z/L) 1/4 and y = (1-13z/L) 1/2 . In many observational studies, the surface temperature is considered to be equal to the surface radiation temperature 27,30,32,37,52,53 . Based on the Stefan-Boltzmann law, the surface radiation temperature, θ 0 , can be obtained from the upward longwave radiation, ↑ R lw (W m −2 ), and the downward longwave radiation, ↓ R lw (W m −2 ): where σ (=5.67 × 10 −8 W m −2 K −4 ) is the Stefan-Boltzmann constant, and ε 0 is the surface emissivity derived from observations. From Eq. (2), ln(z-d) is a variable value on vegetable-covered surfaces due to seasonal changes. For a bare, sandy surface, the value will remain constant. Hence, for bare, sandy surfaces, ln(z 0m ) is a function of u/u * and Ψ m (z/L). Similarly, from Eq.(3), ln(z 0h ) is a function of (θ a -θ 0 )/ θ * and Ψ h (z/L).