Observed Key Surface Parameters for Characterizing Land–Atmospheric Interactions in the Northern Marginal Zone of the Taklimakan Desert, China

An observational data set of the year 2010 at a site in the northern marginal zone of the Taklimakan Desert (TD) was used to analyse the key surface parameters in land–atmospheric interactions in the desert climate of northwest China. We found that the surface albedo (α) and emissivity (ε) were 0.27 and 0.91, respectively, which were consistent with the values obtained based on observations in the hinterland of the TD as well as being similar to the dry parts of the Great Basin desert in North America, where they were comparable to the α and ε values retrieved from remote sensing products. Peak frequency value of z0m was 5.858 × 10−3 m, which was similar to the Mojave Desert, Peruvian desert, Sonoran Desert, HEIFE (Heihe region) Desert, and Badain Jaran Desert. The peak frequency value of z0h was 1.965 × 10−4 m, which was different from those obtained in the hinterland of the TD. The average annual value of excess resistance to heat transfer (kB−1) was 2.5, which was different from those obtained in the HEIFE Gobi and desert, but they were similar to those determined for the Qinghai–Tibetan Plateau and HAPEX-Sahel. Both z0m and z0h varied less diurnally but notably seasonally, and kB−1 exhibited weak diurnal and seasonal variations. We also found that z0m was strongly influenced by the local wind direction. There were many undulating sand dunes in the prevailing wind and opposite to the prevailing wind, which were consistent with the directions of the peak z0m value. The mean values calculated over 24 h for Cd and Ch were 6.34 × 10−3 and 5.96 × 10−3, respectively, which were larger than in the Gobi area, hinterland of the TD and semiarid areas, but similar to HEIFE desert. Under the normal prevailing (NNE–ESE) wind, the mean bulk transfer coefficient Cd and Ch were of the same order of magnitude as expected based on similarity theory. Using the data obtained under different wind directions, we determined the relationships between Cd, Ch, the wind speed U, and stability parameter z/L, and the results were different. Cd and Ch decreased rapidly as the wind speed dropped below 3.0 m s−1 and their minimum values reached around 1–2 m s−1. It should also be noted that the ε values estimated using the sensible heat flux (H) were better compared with those produced using other estimation methods.


Introduction
Interactions between the atmosphere and the underlying surface play key roles in the land-air exchange of mass and energy for weather and climate changes [1,2].In particular, many studies have observed and modeled land-atmospheric interactions in semiarid and arid regions because their natural environments and human populations are especially vulnerable to anomalous weather and climate conditions.
The physical processes of surface mass, energy and water exchange in the atmosphere are affected by several key parameters, such as the surface albedo (α), surface emissivity (ε), aerodynamic and thermal roughness lengths (z 0m and z 0h ), and momentum and sensible heat drag coefficients (C d and C h ) [3][4][5][6].Therefore, it is very important to study the variations in these parameters under different climate conditions.Thus, Boussetta et al. [7] showed that weather forecasts are sensitive to the surface albedo because it controls the partition of surface energy fluxes.In addition, ε is a critical variable for calculating surface longwave radiation and it is used in both land and atmospheric remote sensing.Studies have shown that the surface temperature, outgoing longwave radiation, cooling rates, and frozen surface extent are sensitive to far-infrared surface emissivity [8,9].z 0m and z 0h are defined as the height at which the wind speed becomes zero and at which the extrapolated air temperature is identical to the surface temperature over a homogeneous surface under neutral and thermally stratified conditions [10][11][12], respectively, and they are very important parameters for estimating the momentum, heat, and mass exchange between the surface and atmosphere [13][14][15][16][17][18][19].They are generally estimated directly using the eddy covariance method and satellite data inversion technique [20,21].Many previous studies have shown that z 0h and kB −1 exhibit diurnal and monthly variations [22][23][24][25][26].In surface energy balance calculations, C d and C h are usually obtained from the bulk transfer equations [4].Therefore, accurate estimates of these key parameters are important for land-air interaction observation experiments.However, we still lack a thorough understanding of the key parameters for representing land-atmospheric interactions and their underlying processes, thereby limiting our capacity for modeling and forecasting in desert.
In deserts, low vegetation density, soil surface organic carbon contents, and soil moisture lead to high α values (With more carbon and more soil moisture, the surface be darker, reducing albedo), which lower the net radiation compared with other ecosystems and with the same input of incoming solar radiation [27,28].In addition, the surface energy residual exhibits strong diurnal variations due to the dominance of sensible heat.The α values also exhibit high seasonal variations due to frequent sand events during the sandstorm seasons [29][30][31][32][33][34].The thermal effects of land surface processes play important roles in the formation of sandstorm weather, and dust aerosols in the atmospheric boundary layer can then have significant thermal dynamic and dynamic effects on the atmosphere [35].The TD is an arid region of China and the second largest shifting sand desert on Earth; it has been studied widely based on variations in the characteristic parameters of land-atmospheric interactions and parameterization schemes [8,[36][37][38], although we have limited knowledge of the spatial variations in the key parameters that dominate land-air coupling across this giant desert.
The TD (total area = ca 338,000 km 2 ) is located in the Tarim Basin, west China (Figure 1).It is the largest desert in China and it comprises 85% shifting sand with most types of sand dunes [39][40][41].The TD strongly influences the climate and environment in the northwest arid region of China.Dust storms in this region also alter the energy balance and atmospheric chemistry [42,43].In a recent study, Liu et al. [8,44] estimated the values of α, ε, z 0m , and z 0h in the hinterland of the TD.However, the information on key surface parameters in the northern marginal zone of the TD remains limited including α, ε, z 0m , z 0h , C d , and C h .Thus, in this study, we aimed to analyze the key surface parameters (α, ε, z 0m , z 0h , kB −1 , C d , and C h ) at Xiaotang experimental field observation station in the northern marginal zone of the TD.In contrast to Tazhong (with an underlying surface of fine sand) where extensive observations have been made [8] and the site has slight distortions due to surrounding vegetation (Gongliu), the station (the underlying surface is an ancient river bed with a crust of clay and sand) is affected more evidently by nearby vegetation (Populus diversifolia).The important parameters ( α , ε, z0m, z0h, Cd, and Ch) were analyzed in a spatially inhomogeneous environment and compared with data obtained for other deserts or similar The important parameters (α, ε, z 0m , z 0h , C d , and C h ) were analyzed in a spatially inhomogeneous environment and compared with data obtained for other deserts or similar ecosystems.The remainder of this paper is organized as follows.In Section 2, we describe the observational site, instruments and data collection processes, quality control procedures, and methods used for deriving the key surface parameters.In Section 3, we present the observational analysis and comparisons with results obtained in Tazhong.Finally, we give our conclusions based on this analysis in Section 4.

Site Description
Observations obtained from the Xiaotang Land-Atmosphere Interaction Observation Station (933 m a.s.l., 40 • 48 N, 84 • 18 E) were used in this study.This station is located in the transitional zone (from oasis to desert ecosystem) over the northern margin of TD (Figure 1a-c).The station is in a sandy area that comprises sand and ancient river beds (Figure 1b,c).
According to observations obtained from Xiaotang meteorological station in 2010, the annual average temperature was 12.0 • C, and the average temperatures were 28.6 • C in July and −16.3 • C in January.This area is dominated by an extremely dry climate where the annual average precipitation is only 86.5 mm but the potential evaporation is 3025.2mm.The prevailing wind direction is ENE (Figure 1d).

Data
The measurement levels, instruments, near-surface meteorological variables, turbulent fluxes, and surface radiation spectrum at the station are shown in Table 1.The observations of meteorological variables and turbulent fluxes used in this study covered the period from 5 January to 31 December 2010.The observations of wind speed and direction at 30-min intervals were acquired by a Metone 010C/020C sensor at a height of 10 m and using a CR1000 data logger.Observations of radiation at 30-min intervals were obtained used a CRN-1 mounted at a height of 1.5 m and with a CR1000 data logger.Observations of the air temperature at 30-min intervals were made using a Vaisala HMP45D sensor at heights of 0.5, 1, 2, 4, and 10 m, and with a CR1000 data logger.Observations of the soil temperature were acquired at 30-min intervals using Campbell 109L sensors at the surface soil (half of the sensor was above the surface) at heights of 10, 20, and 40 cm above the surface, and with a CR1000 data logger.Raw turbulent flux data were obtained using a Campbell CSAT3/Licor 7500 sensor at a height of 3.0 m and with a CR5000 data logger.
The raw turbulent data were processed to obtain 30-min interval turbulent flux averages using EddyPro (Eddy Covariance Processing Software) software (the parameters employed in EddyPro are shown in Table 2).In addition, a series of quality control procedures were applied to the turbulent flux data, i.e., noise removal, coordinate rotation, stability testing, variance similarity analysis, and removal of out-of-correct-range data.The statistical analyses showed that the proportions of missing data and unqualified data (abnormally high or low data) were 8.2% and 28.5%, respectively.The unqualified data were excluded.The Fourier transform infrared spectrometer (FTIR) was produced by Designs and Prototypes Company (USA).The instrument comprised of one Michelson interferometer, one detector, one black body, one electric current converter, one block reflector gold board, and one built-in individual computer.The computer program was used for sampling, processing, and the storage of samples on an individual computer.The FTIR had a fast recording speed, high signal-to-noise ratio, good sensitivity, and very little stray radiation, where the light flux = 0.016 cm 2 sr, spectral range = 2-16 µm, and spectral resolution = 2-24 cm −1 .The measurement results obtained had standard deviations less than 1% [45].
The principle of FTIR involves converting the received infrared spectra M s (λ) into spectral power Vs(λ) by using the internal photoelectric effect: where λ is the wavelength, r(λ) is the linear response of the FTIR, and M 0 (λ,T inst ) is radiation at the temperature of the FTIR.The FTIR was calibrated before measuring the radiation from the sample.M s (λ) is the sample radiation measured by the FTIR after calibration.
According to Kirchhoff's law, for opaque objects such as surfaces, the sum of the absorptivity and reflectivity (R s ) is 1, and the absorptivity is equal to ε; thus, The distance from the ground measurement sample to the sensor was less than 1 m, where the upward radiation could be treated as negligible at this distance and the atmospheric transmittance was considered to be 1.Thus, the radiation from the sample was expressed as: where is the radiation that is reflected by environmental influences such as the atmosphere.M dwτ (λ) is the fallout radiation comprising downward radiation in the atmosphere reflected by the sample detected with the FTIR.B(λ, T s ) is black body radiation.
ε s (λ) was obtained by transforming formula (3) [45,46], as follows, In Equation (4), M s (λ) and M dwτ (λ) were obtained by the FTIR.Thus, sample calibration was necessary before measuring M s (λ) and M dwτ (λ) to ensure the accuracy of the measurements.Therefore, ε s (λ) was converted into ε: where λ 1 and λ 2 are wavelength ranges with thermal infrared atmospheric window values of 8 and 14, respectively.In order to facilitate the computation, the integral equation was discretized as follows.
In fact, in order to improve the accuracy, the wavelength range interval of 8-14 µm was divided into 375 ∆λ.
Land use is the same in all directions in TD (large scale), but the sand dunes are large and undulating near to the station.There is a flat and bare ancient riverbed SE (NW) of the station's tower, which was consistent with the directions of 90% (70%) of the flux source area.The flux source area could be observed well in the windward direction (prevailing wind direction) [47].In this study, the bulk transfer coefficients (C d and C h ) were calculated based only on data measured in the NNE-ESE wind directions (local dominant wind direction) (Figure 1d-f).We then compared the results with those obtained in all wind directions.
The ground surface soil heat flux was calculated according to the soil temperature and moisture gradients [48].Due to errors in the land surface temperature during the sand dust season, the ground surface temperature was calculated using the observed ε and radiation values.

Theory and Methodology
According to many other observational studies [49,50], the energy balance closure is not achieved completely for different underlying types.The energy balance closure is usually formulated as: R n − G 0 = H + LE, where R n is the net radiation flux, G 0 is the ground surface soil heat flux, H is the sensible heat flux, and LE is the latent heat flux.The ratio of the energy balance closure is usually formulated as: EBR = ∑ (H + LE)/∑ (R n − G 0 ) * 100%.The energy balance is formulated relative to the residuals as: In the analysis, α was calculated from the observed surface solar radiation components [51] using Equations ( 7) and ( 8): where S ↑ is shortwave upward radiation, S ↓ is shortwave downward radiation, α is the mean of α (calculated using the weighted mean method), and the subscript i is the time index.
In this study, we focused on the measurements in sunny and dry weather.In order to calibrate the radiation, the FTIR had to be calibrated every 10-20 min to a black body.The temperature of the cold black body was 10 • C lower than the environmental temperature, whereas the temperature of the hot black body was 10 • C higher than the surface temperature.After setting the temperatures of the cold black body and the hot black body, radiation spectrum data were measured for the cold and hot black body and saved.The accuracy achieved for the black body emissivity was 0.994-0.998± 0.002 and the accuracy of the temperature was ±0.1 • C. Thus, the error caused by the black body was less than 0.004.The temperature fluctuations according to the interferometer remained within 0.1 • C and the calibration error for the black body was less than 0.002 [46].In order to reduce the interference in the instrument's noise signal, we set the spectrum stacking number at 10 times and the average values were obtained.
In this study, we performed three steps to minimize the error during operation in order to obtain the surface emissivity spectrum with high accuracy: (1) we measured the cold black body, hot black body, and diffuse reflection radiation from a gold plate; (2) we measured the surface radiation; and (3) we repeated step (1).These three steps were performed as quickly as possible, where we limited the time to 10 min for each group of measurements.Radiation correction was performed in steps ( 1) and (3) to evaluate the influence of the environment over time on the surface spectral radiation, thereby reducing the error.The surface temperature of the diffuse gold plate was measured using a thermoelectric coupling thermometer.In general, the average value was taken based on five measurements.
For the desert surface, Korb et al. [45] suggested that the maximum emission rate for the band fitted (black body radiation spectrum fitted to the surface radiation spectrum) at 7.45-7.65µm is 0.95.By using this method to obtain the surface radiation temperature, the surface emission spectra were acquired for a thermal infrared window at wavelengths of 8-14 µm in the TD with high efficiency and accuracy.
We used the physically-based and semi-empirical methods of Yang et al. [11] to calculate ε.The computed surface sensible heat (H sfc ) was compared with the observed surface sensible heat (H obs ) to fit the ε values.H sfc can be obtained from: is the heat capacity of air at constant pressure, ρ is the air density (kg m −3 ), T 0 is the aerodynamic surface temperature (K), T a is the air temperature (K), and r h is the aerodynamic resistance for heat (s m −1 ) , where Pr is the Prandtl number (=1 if z/L ≥ 0 and 0.95 if z/L < 0), L[≡ T a u 2 * /(kgT * )] is the Obukhov length, u * is the frictional velocity (m s −1 ); T * [≡ −H/(ρc p u * )] is the frictional temperature; ψ m (ζ) and ψ h (ζ) are the integrated stability correction function for momentum transfer and temperature profiles, respectively; k (=0.4) is the von Kármán constant, R ↑ lw and R ↓ lw are the upward longwave radiation and downward longwave radiation, respectively.σ (5.677 × 10 −8 W m −2 K −4 ) is the Stefan-Boltzmann constant near neutral (i.e., T 0 = T a ).The difference between H sfc and H obs is sensitive to the value of ε.Using an iterative algorithm (ε values from 0.8 to 1.0 with a step width of 0.001), the ε value was derived by minimizing the root mean square (RMS) between the calculated H sfc and observed H obs .
In the analysis, z 0m was estimated using Equations ( 9) to (11): where z is the observed height (m) and d is the displacement height (m), which was negligible in our study with no vegetation coverage.According to Dyer ([52], parameters 16 and 5 in Equations ( 10) and ( 11) are consistent with k = 0.40, while u is the average wind velocity (m s −1 ).Under unstable conditions [12,53], the following equation can be used: where, x = (1 − 15.2ζ) 1/4 .Under stable conditions [54,55], the equation is as follows.
The value of z 0h is difficult to determine because it cannot be measured directly, but it is possible to derive z 0h from the equations for H obs [23], where it can be obtained from the flux-gradient relationships in a surface layer based on Monin-Obukhov (M-O) similarity theory [56].The following equation can be used.
Under unstable conditions: Under stable conditions: where In the analysis, we employed the same procedure used by Yang et al. [11] for data quality control, as follows: (1) excluding periods when the absolute value of |H obs | < 10 W m −2 , which corresponds to the time with a low solar elevation angle when the observation error was large; (2) excluding periods when z 0h > z 0m , which is physically unrealistic; and (3) excluding data from sandy, rainy, and snowy periods for C d and C h .
In our analysis, C d and C h were calculated using two methods.First, we followed the Eddy correlation method using Equation ( 15): where τ is the surface stress (kg m −1 s −2 ).The other method for calculating C d and C h is based on M-O similarity theory where C d and C h are expressed as [57,58]: where Ψ m (z/L) and Ψ h (z/L) are the integration forms of the M-O similarity functions ϕ m and ϕ h at the station, respectively, and the equation is: where ∆ is the difference symbol, z 1 and z 2 are the height above ground (in this analysis, z 1 = 2 m, mboxemphz 2 = 4 m), θ * is the turbulent temperature scale, and θ is the potential temperature.
The lack of observations across the whole TD prevented us from capturing the spatial distributions of these surface parameters based on in situ observations alone, so we also compared our site observations with the values of α and ε from remote sensing products.This comparison allowed us to use in situ observations to calibrate the remote sensing products and to produce improved surface parameter estimates over a large spatial domain.In this analysis, Landsat 8 data from NASA were used to compare the observation data for α and ε.The empirical method described by Qin et al. [59] was used to calculate α and ε from the Landsat data.We only used data from 19 August 2013 at 11:00 (Local time), which was a sunny day.The solar elevation angle at the station was 59.05 • at 11:00 on 19 August 2013.

Surface Energy Balance Closure
The energy balance closure was not achieved completely at the station (Figure 2a).During the daytime (positive net radiation), EBR was about 40-60%, whereas it was less than 25% at night (negative net radiation) with no correction for the soil heat flux (Figure 2a). Figure 2b shows that during the daytime, EBR was about 60-70%, whereas it was 25-35% at night with G 0 .Table 3 shows during the dust season (March and April), EBR was less than 60%, whereas it was more than 60% during other seasons (it was not possible to obtain statistics for EBR due to the energy balance data synchronization matching problem in October, November, December, January, and February).Figure 2c shows the mean diurnal cycles of the energy fluxes.The monthly mean peak H, LE, R n , and 2.5-cm depth soil heat flux (G 2.5 ) reached 122.1, 27.7, 295.1, and 35.1 W m −2 , respectively; H and R n exhibited different patterns compared with LE and G 2.5 .Figure 2d,g shows that the monthly mean peak G 0 and 8-cm depth soil heat flux (G 8 ) reached 64.2 W m −2 and 22.3 W m −2 , respectively.The slope value of H + LE and R n − G 2.5 was around 0.44, and that of H + LE and R n − G 0 was around 0.61.In the regressions, G 0 increased the slope of the ordinary least squares regression by 39%, which is why it was necessary to determine the ground heat storage in the desert in this study.Over the desert, LE was a very small component of the surface energy balance with the weak diurnal variation where the maximum values for LE and G 0 (G 2.5 and G 8 ) occurred in the morning and noon (afternoon), respectively, because the soil heat capacity was larger than the air heat capacity.EBR was typically lower when the friction velocity was below 0.19 m s −1 , but greater when the friction velocity was above the 0.19 m s −1 (Figure 2h).EBR was typically lower when the Bowen ratio was above 4.6, but greater when the Bowen ratio was below the 4.6 (Figure 2h).Wilson et al. [49] suggested that the main reasons for the energy imbalance are systematic errors associated with sampling mismatch, systematic instrument bias, neglected energy sinks, low and high frequency loss of turbulent fluxes, and horizontal and/or vertical advection of heat and water vapor.Previous studies have suggested that the energy balance is almost closed in a flat desert region.However, the energy imbalance is large in the TD.At the station, the neglected surface soil heat flux was corrected based on the soil temperature and humidity; thus, it was mainly due to the neglected low frequency or high frequency turbulent flux.In spite of this, landscape heterogeneity at the station cannot be ignored as a contributor to incomplete energy balance closure [50].

Relationships between ϕ m , ϕ h , and z/L
Figure 3 shows the observed and calculated functions ϕ m and ϕ h against z/L in the NNE-ESE wind direction.In this analysis, data validation showed that data were in agreement with the theoretical relationship between the functions ϕ m , ϕ h , and z/L at the station under the same wind direction.
Atmosphere 2018, 9, x FOR PEER REVIEW 11 of 26 corrected based on the soil temperature and humidity; thus, it was mainly due to the neglected low frequency or high frequency turbulent flux.In spite of this, landscape heterogeneity at the station cannot be ignored as a contributor to incomplete energy balance closure [50].

Relationships between φm, φh, and z/L
Figure 3 shows the observed and calculated functions φm and φh against / in the NNE-ESE wind direction.In this analysis, data validation showed that data were in agreement with the theoretical relationship between the functions φm, φh, and z/L at the station under the same wind direction.The relationships between φm, φh, and / in the NNE-ESE wind direction were obtained as follows (Equations ( 18) and ( 19)), where these coefficients are based on the least square fits.

Surface Albedo α
To remove any possible contamination of clouds from α measurements, the data acquired during sunny periods were used for estimating this parameter.Figure 4a shows the relationship between α and the solar elevation angle (h) during sunny periods, which demonstrates that there was a significant exponential relationship between α and h (when taking it as a value) where α was fitted best by the following numerical value equation.
According to Figure 4a,b, α was 0.27 when h > 15°, which is in good agreement with the values obtained for other deserts, as well as the measurements from Tazhong station at the center of the TD [3,40,36,60,61].α (weighted mean) was 0.31 and 0.30 on the sunny days and cloudy days, The relationships between ϕ m , ϕ h , and z/L in the NNE-ESE wind direction were obtained as follows (Equations ( 18) and ( 19)), where these coefficients are based on the least square fits.

Surface Albedo α
To remove any possible contamination of clouds from α measurements, the data acquired during sunny periods were used for estimating this parameter.Figure 4a shows the relationship between α and the solar elevation angle (h) during sunny periods, which demonstrates that there was a significant exponential relationship between α and h (when taking it as a value) where α was fitted best by the following numerical value equation.
Figure 4b shows the relationship between α and h during cloudy periods by the following numerical value Equation (21).
According to Figure 4a,b, α was 0.27 when h > 15 • , which is in good agreement with the values obtained for other deserts, as well as the measurements from Tazhong station at the center of the TD [3,36,40,60,61].α (weighted mean) was 0.31 and 0.30 on the sunny days and cloudy days, respectively.According to the regression-like analysis of the outgoing versus incoming shortwave radiation during sunny and cloudy periods (Figure 4c,d), the results were similar to the weighted mean values, but the outgoing shortwave radiation was slightly lower when the assumptions underlying the weighted mean were met.The diurnal variations in α during sunny and cloudy periods are shown in Figure 4e.α was higher in the morning and evening, but lower in the daytime, where the values was 0.27 from 9:00 to 15:00 (LST) during sunny and cloudy periods.α was slightly lower during cloudy periods than sunny periods, especially in the morning.We also compared the observed values with the values of α derived from the Landsat 8 satellite over this area.The retrieved albedos (from Landsat 8 data) were close to the observed values at both the station and Tazhong, where they ranged from 0.24 to 0.30 near the station (Figure 4f).Therefore, our results suggest that the satellite-derived values of α from remote sensing products are highly reliable for use in the TD.
respectively.According to the regression-like analysis of the outgoing versus incoming shortwave radiation during sunny and cloudy periods (Figure 4c,d), the results were similar to the weighted mean values, but the outgoing shortwave radiation was slightly lower when the assumptions underlying the weighted mean were met.The diurnal variations in α during sunny and cloudy periods are shown in Figure 4e.α was higher in the morning and evening, but lower in the daytime, where the values was 0.27 from 9:00 to 15:00 (LST) during sunny and cloudy periods.α was slightly lower during cloudy periods than sunny periods, especially in the morning.We also compared the observed values with the values of α derived from the Landsat 8 satellite over this area.The retrieved albedos (from Landsat 8 data) were close to the observed values at both the station and Tazhong, where they ranged from 0.24 to 0.30 near the station (Figure 4f).Therefore, our results suggest that the satellite-derived values of α from remote sensing products are highly reliable for use in the TD.
Atmosphere 2018, 9, x FOR PEER REVIEW 13 of 26 variations in the surface albedo on sunny days and cloudy days.(f) Retrieval of the surface albedo from Landsat data.

Surface Emissivity ε
Figure 5a shows that the value of ε corresponding to the minimum RMS error value for Hsfc against Hobs was 0.85 under near neutral conditions.Thus, we determined the sensitivity of ε to Hobs.Table 4 shows the difference in ε with Hobs according to the method employed, which indicates that ε decreased (increased) as Hobs increased (decreased).For instance, ε increased by 2, 7.6, 9, 11.6, 12.8, 13.5, and 17.4% when Hobs increased by 5, 10, 15, 20, 25, 30, and 35%, respectively, whereas ε decreased by 5.2, 8.9, 16.    Figure 5b shows that the dependence of ε on the density according to the surface temperature method was 0.85.Table 5 shows that the values of ε averaged over a temperature difference of 0-0.1 • C in group 1 (0.839) and group 2 (0.818), estimated according to temperature gradients in near zero conditions, which was about 0.83, and the mean value of ε averaged over a temperature difference of 0.1-0.2• C in group 1 (0.849) and group 2 (0.825) estimated according to temperature gradients in near zero conditions, which was about 0.84, and so on.We also selected the same day to retrieve the value of ε for these deserts using the Landsat 8 data.Figure 5c shows that the retrieved ε values (from Landsat 8 data) were generally close to the observed value, with a range of 0.912 to 0.916 at the station.The observed value of ε was 0.91 according to the FTIR data, which is very close to that reported for other deserts, i.e., 0.91 to 0.97 [36,40,62].However, ε estimated based on sensible heat flux and surface temperature was both 0.85, which was lower than the observed value (0.91).
Figure 6a,b shows that the distributions of z 0m and z 0h were approximately log-normal.The peak frequency for ln(z 0m ) was −5.14, which was equivalent to z 0m = 5.858 × 10 −3 m in the normally distributed histogram.Similar values were obtained for deserts such as the Mojave Desert, Peruvian desert, Sonoran Desert, HEIFE (Heihe region) Desert, and Badain Jaran Desert, which were all in the order of 10 −3 m [40, [62][63][64].Figure 6b shows that the peak frequency was −8.54 for ln(z 0h ), which was equivalent to z 0h = 1.965 × 10 −4 m in the normally distributed histogram.Data with magnitudes of 10 −3 (z 0m ) and 10 −4 (z 0h ) were retained in the analysis.
Many previous studies have investigated the relationships between z 0m , the wind speed, the wind direction [12,65], and u * [66].z 0m does not change with the wind direction in areas with a flat and uniform underlying surface, whereas z 0m might be related to the wind direction in areas with a non-uniform underlying surface (rough elements are not evenly distributed in different directions).Figure 7a,b shows the changes in z0m and z0h with different wind directions at the station.z0m had a double peak when the wind was in the E direction and the SW-WSW wind direction.z0m was 3.2 × 10 −2 m when the wind was in the E wind direction and 2.8 × 10 −2 to 3 × 10 −2 m when it was in the SW-WSW wind direction at the station.In addition, z0h had a double peak when the wind was in the ENE wind direction and the NNW wind direction.z0h was 5.3 × 10 −3 m when the wind was in the ENE wind direction and 5.6 × 10 −3 m when it was in the NNW wind direction at the station.
Figure 7c,d shows the changes in z0m and z0h with different wind directions at the station, where the data magnitudes were 10 −3 and 10 −4 for z0m and z0h, respectively.z0m had two high values in areas when the wind was in the NE-E wind direction and the S-WSW wind direction.z0m was 5.6-6.3 × 10 −3 m when the wind was in the NE-E wind direction and 5.1-6.1 × 10 −3 m when it was in the S-WSW wind direction at the station.However, z0h only changed weakly with the wind direction.Figure 1d shows the characteristic rough surface elements at the station.
Therefore, the obvious changes in z0m with the wind direction at the station were affected mainly by the spatial inhomogeneity.The distributions of the surrounding sand dunes at the station are not uniform.In the NE-E and SSE-WSW directions from the tower (red dot in Figure 1), there are many undulating sand dunes, which were consistent with the directions of the peak z0m value.A flat and bare ancient riverbed is located in the SE direction from the tower, which is consistent with the directions of the 90% flux source area and the small values of z0m.Thus, the results were roughly consistent with the fact that the peak z0m values occurred with the prevailing wind and opposite to the prevailing wind, whereas they did not occur in the 90% flux source area at the station.Figure 7a,b shows the changes in z 0m and z 0h with different wind directions at the station.z 0m had a double peak when the wind was in the E direction and the SW-WSW wind direction.z 0m was 3.2 × 10 −2 m when the wind was in the E wind direction and 2.8 × 10 −2 to 3 × 10 −2 m when it was in the SW-WSW wind direction at the station.In addition, z 0h had a double peak when the wind was in the ENE wind direction and the NNW wind direction.z 0h was 5.3 × 10 −3 m when the wind was in the ENE wind direction and 5.6 × 10 −3 m when it was in the NNW wind direction at the station.
Figure 7c,d shows the changes in z 0m and z 0h with different wind directions at the station, where the data magnitudes were 10 −3 and 10 −4 for z 0m and z 0h , respectively.z 0m had two high values in areas when the wind was in the NE-E wind direction and the S-WSW wind direction.z 0m was 5.6-6.3 × 10 −3 m when the wind was in the NE-E wind direction and 5.1-6.1 × 10 −3 m when it was in the S-WSW wind direction at the station.However, z 0h only changed weakly with the wind direction.Figure 1d shows the characteristic rough surface elements at the station.
Therefore, the obvious changes in z 0m with the wind direction at the station were affected mainly by the spatial inhomogeneity.The distributions of the surrounding sand dunes at the station are not uniform.In the NE-E and SSE-WSW directions from the tower (red dot in Figure 1), there are many undulating sand dunes, which were consistent with the directions of the peak z 0m value.A flat and bare ancient riverbed is located in the SE direction from the tower, which is consistent with the directions of the 90% flux source area and the small values of z 0m .Thus, the results were roughly consistent with the fact that the peak z 0m values occurred with the prevailing wind and opposite to the prevailing wind, whereas they did not occur in the 90% flux source area at the station.The diurnal variations in the means of ln(z0m), ln(z0h), and kB −1 are shown in Figure 8.The range for z0m was 1.1 × 10 −3 to 8.4 × 10 −3 m, the range for z0h was 1.2 × 10 −4 to 6.1 × 10 −4 m, and kB −1 varied from 0.9 to 3.9.The ln(z0m), ln(z0h), and kB −1 terms were nearly constant during the daytime, but they fluctuated at sunrise and sunset (Figure 8).z0m depended on u* and z/L (Figures 2h and 10c).ln(z0m) was higher in July (5.5 × 10 −3 m), but lower in November (4.2 × 10 −3 m) (Figure 9a).Between December and June/July, the gradual increases in the wind velocity/surface-air temperature increased and LE (Figure 9b-d).Blyth and Dolman [67] showed that the vegetation structure, meteorological conditions, and soil surface resistance were the main factors that affected z0h, and thus ln(z0h) exhibited seasonal variations.We found that z0h was higher in April and September (4.8 × 10 −3 m and 4 × 10 −3 m, respectively), but lower in January and February (3.4 × 10 −3 m and 3.2 × 10 −3 m, respectively) (Figure 9a).The monthly variations in ln(z0h) correlated well with those in LE, but with a phase shift (Figure 9b).The diurnal variations in the means of ln(z 0m ), ln(z 0h ), and kB −1 are shown in Figure 8.The range for z 0m was 1.1 × 10 −3 to 8.4 × 10 −3 m, the range for z 0h was 1.2 × 10 −4 to 6.1 × 10 −4 m, and kB −1 varied from 0.9 to 3.9.The ln(z 0m ), ln(z 0h ), and kB −1 terms were nearly constant during the daytime, but they fluctuated at sunrise and sunset (Figure 8).z 0m depended on u * and z/L (Figure 2h and Figure 10c).ln(z 0m ) was higher in July (5.5 × 10 −3 m), but lower in November (4.2 × 10 −3 m) (Figure 9a).Between December and June/July, the gradual increases in the wind velocity/surface-air temperature increased H and LE (Figure 9b-d).Blyth and Dolman [67] showed that the vegetation structure, meteorological conditions, and soil surface resistance were the main factors that affected z 0h , and thus ln(z 0h ) exhibited seasonal variations.We found that z 0h was higher in April and September (4.8 × 10 −3 m and 4 × 10 −3 m, respectively), but lower in January and February (3.4 × 10 −3 m and 3.2 × 10 −3 m, respectively) (Figure 9a).The monthly variations in ln(z 0h ) correlated well with those in LE, but with a phase shift (Figure 9b).In addition, the monthly variations in z 0m and z 0h were related to the surface conditions.For example, the station often experiences frequent dusty weather events during the spring and summer, and the number of dusty days in the spring (March to May) accounts for about 46% of the annual dusty days.Thus, z 0h varied from 3.17 × 10 −4 m in February to 4.84 × 10 −4 m in April with a mean of 3.76 × 10 −4 m.In contrast to ln(z 0h ), the phase was the opposite for kB −1 .The kB −1 term was nearly constant in all seasons and the monthly mean (based on January, February, March, April, May, June, July, August, September, November, and December) was 2.5.These results are similar to those obtained by Brutsaert [68] who found that the mean kB −1 value was 2.3 for tall vegetation and it ranged between −1.6 to −0.16 for bare soil.Furthermore, our kB −1 results for smooth surfaces were considerably different from that reported by Kondo [69] (i.e., −1.1).On the Qinghai-Tibetan Plateau, the kB −1 values were 2.5 and 2.36 in North PAM (Portable Automated Mesonet) and Anduo (both plateau meadows) [25], while in a semiarid region with degraded grassland and farmland, the kB −1 values ranged among 7.1-9.2and 8.6-11.0,respectively [70].Verhoef et al. [23] showed that the kB −1 values were 6.9, 8.1, 12.4, and −0.9 in HAPEX-Sahel (Hydrologic-Atmospheric Pilot Experiment in the Sahel), vineyard, sparse vegetation, and bare soil areas, respectively.This is because z 0m increases with the vegetation height, so for a vegetated underlying surface, the vegetation height basically determines the magnitude of the roughness length.Therefore, with a tall vegetation underlying surface, the average value of kB −1 was large.z 0m was smaller in the desert than other areas, but the surface temperature difference and H were larger in the desert than other areas.Thus, kB −1 was different from that in other areas.

Bulk Transfer Coefficients C d and C h
According to Equations ( 15) and ( 16), the turbulent exchange is affected by C d and C h , which are the key factors in turbulent flux parameterization schemes over different underlying surfaces.We estimated C d and C h according to Equation (15) (Eddy correlation method) and Equation ( 16) (M-O similarity theory), respectively.Figure 10 and Figure 13a,b shows the results obtained using Equation ( 15), and Figure 13c,d shows the results produced with Equation ( 16).In addition, the monthly variations in z0m and z0h were related to the surface conditions.For example, the station often experiences frequent dusty weather events during the spring and summer, and the number of dusty days in the spring (March to May) accounts for about 46% of the annual dusty days.Thus, z0h varied from 3.17 × 10 −4 m in February to 4.84 × 10 −4 m in April with a mean of 3.76 × 10 −4 m.In contrast to ln(z0h), the phase was the opposite for kB −1 .The kB −1 term was nearly constant in all seasons and the monthly mean (based on January, February, March, April, May, June, July, August, September, November, and December) was 2.5.These results are similar to those obtained by Brutsaert [68] who found that the mean kB −1 value was 2.3 for tall vegetation and it ranged between −1.6 to −0.16 for bare soil.Furthermore, our kB −1 results for smooth surfaces were considerably different from that reported by Kondo [69] (i.e., −1.1).On the Qinghai-Tibetan Plateau, the kB −1 values were 2.5 and 2.36 in North PAM (Portable Automated Mesonet) and Anduo (both   8 and 10a,b shows that the changes in z 0m and z 0h were higher after sunrise and sunset, where they were consistent with the variations in C d and C h .Figure 10a,b shows that the values of C d and C h were higher (lower) in the daytime (at night).According to the strong diurnal variations in the surface temperature, z/L < 0 after sunrise, but it was stable (z/L > 0) after sunset.Thus, u/u * was lower (higher) in the daytime (at night) (Figure 10d).As a result, C d and C h also exhibited diurnal variations, where they responded to variations in the surface roughness length, atmospheric stratification and wind speed.
According to our observations, the mean values calculated over 24 h for C d and C h were 6.34 × 10 −3 and 5.96 × 10 −3 , respectively.The mean values calculated during the daytime for C d and C h were 8.30 × 10 −3 and 1.11 × 10 −2 , respectively.The mean values calculated during the nighttime for C d and C h were 2.44 × 10 −3 and 8.50 × 10 −4 , respectively.Compared with previous studies, the mean C d and C h values were larger at the station than in the Gobi area (e.g., Dunhuang and HEIFE), the hinterland of the TD (e.g., Tazhong), and semiarid areas (e.g., Tongyu), but similar to the C d and C h values for the HEIFE desert [8,[70][71][72][73][74].These differences may be explained by the combined effects of the surface roughness and overlying atmospheric stratification.With the same underlying surface, the atmospheric stratification was more unstable and the bulk transfer coefficients were larger.C d was higher in near zero conditions with greater values of z 0m [70].
Figure 11 shows the monthly variations in C d and C h .Figure 11a,b shows that the highest values for C d and C h occurred in May and March (0.99 × 10 −2 and 0.54 × 10 −2 , respectively) whereas the lowest were in February and January (0.30 × 10 −2 and 0.27 × 10 −2 , respectively) because the atmospheric stratification was more unstable in the summer (the air temperature increased in the spring and summer).The order was C d > C h (i.e., the momentum transfer was greater than the heat transfer) in all months except February (the standard deviation of C d was high in February) and the average ratio was 1.8.According to our observations, the mean values calculated over 24 h for Cd and Ch were 6.34 × 10 −3 and 5.96 × 10 −3 , respectively.The mean values calculated during the daytime for Cd and Ch were 8.30 × 10 −3 and 1.11 × 10 −2 , respectively.The mean values calculated during the nighttime for Cd and Ch were 2.44 × 10 −3 and 8.50 × 10 −4 , respectively.Compared with previous studies, the mean Cd and Ch values were larger at the station than in the Gobi area (e.g., Dunhuang and HEIFE), the hinterland of the TD (e.g., Tazhong), and semiarid areas (e.g., Tongyu), but similar to the Cd and Ch values for the HEIFE desert [8,[70][71][72][73][74].These differences may be explained by the combined effects of the surface roughness and overlying atmospheric stratification.With the same underlying surface, the atmospheric stratification was more unstable and the bulk transfer coefficients were larger.Cd was higher in near zero conditions with greater values of z0m [70].
Figure 11 shows the monthly variations in Cd and Ch. Figure 11a,b shows that the highest values for Cd and Ch occurred in May and March (0.99 × 10 −2 and 0.54 × 10 −2 , respectively) whereas the lowest were in February and January (0.30 × 10 −2 and 0.27 × 10 −2 , respectively) because the atmospheric stratification was more unstable in the summer (the air temperature increased in the spring and summer).The order was Cd > Ch (i.e., the momentum transfer was greater than the heat transfer) in all months except February (the standard deviation of Cd was high in February) and the average ratio was 1.8.The variation in z0h was low (from 5.11 × 10 −5 m to 7.90 × 10 −4 m), so both Equation (15) and Figure 9 indicate that the bulk transfer coefficients were mainly dependent on the wind speed and z/L, where the latter was dominated by the difference in temperature between the land surface and the atmosphere.According to our analysis, there was a good power function between the values of Cd and Ch and the wind speed (Figure 12).In general, Figure 12 suggests that Cd and Ch increased rapidly as the wind speed decreased under unstable conditions (z/L < 0) when the wind speed was less than 4.5 m s −1 .
By contrast, Cd and Ch exhibited nonlinear variations with the wind speed when the wind speed > 4.5 m s −1 , where Cd and Ch approached constant values regardless of the wind speed while the atmosphere stratification tended to be neutral.These findings are similar to those reported by Feng et al. [70] who found that Cd and Ch decreased as the wind speed increased in unstable conditions and they were close to constant values as the wind speed increased with degraded grassland and cropland surfaces.The variation in z 0h was low (from 5.11 × 10 −5 m to 7.90 × 10 −4 m), so both Equation ( 15) and Figure 9 indicate that the bulk transfer coefficients were mainly dependent on the wind speed and z/L, where the latter was dominated by the difference in temperature between the land surface and the atmosphere.According to our analysis, there was a good power function between the values of C d and C h and the wind speed (Figure 12).In general, Figure 12 suggests that C d and C h increased rapidly as the wind speed decreased under unstable conditions (z/L < 0) when the wind speed was less than 4.5 m s −1 .
By contrast, C d and C h exhibited nonlinear variations with the wind speed when the wind speed >4.5 m s −1 , where C d and C h approached constant values regardless of the wind speed while the atmosphere stratification tended to be neutral.These findings are similar to those reported by Feng et al. [70] who found that C d and C h decreased as the wind speed increased in unstable conditions and they were close to constant values as the wind speed increased with degraded grassland and cropland surfaces.less than 4.5 m s −1 .
By contrast, Cd and Ch exhibited nonlinear variations with the wind speed when the wind speed > 4.5 m s −1 , where Cd and Ch approached constant values regardless of the wind speed while the atmosphere stratification tended to be neutral.These findings are similar to those reported by Feng et al. [70] who found that Cd and Ch decreased as the wind speed increased in unstable conditions and they were close to constant values as the wind speed increased with degraded grassland and cropland surfaces.We found that the estimates of C d and C h varied with different wind directions.Under the NNE-ESE wind direction, the values of C d and C h decreased (increased) as the wind speed increased under stable conditions (unstable conditions).For these wind directions, the values of C d or C h had the same magnitude using the Eddy correlation method and the revised M-O similarity function method.These results were related mainly to the flux source at the station, i.e., the maximum flux source area in 90% of cases was in the E-SE direction from the 3-m tower.The SW direction flux was affected by the 3-m tower.By contrast, the areas to the SE and NW from the 3-m tower comprised a flat ancient river bed, and those in the NE and SW directions from the 3-m tower were high and rolling dunes.

Discussion
Compared with other underlying surfaces around the world, desert is quite distinctive in the physical process of the underlying surface, especially the heating process between the land surface and the atmosphere.Desert has higher α, lower ε, z 0m and z 0h , compared to humid areas.Due to the remarkable land heating effect on the atmosphere in the desert, the mixing layer develops much deeper than other areas, and the Bowen ratio shows significant variations [76]).Sandstorms occur frequently in desert, resulting in obvious changes in the surface absorption of solar radiation and the proportion of sensible heat flux, latent heat flux and soil heat flux by the direct and indirect radiation effects of dust aerosols.Also, parameters (α, ε, z 0m , z 0h , kB −1 , C d , and C h ) were influenced in sandstorms weather.Consequently, it will influence the land-atmosphere interaction in the desert greatly.In addition, because of the extremely low precipitation and strong evaporation in desert, there is few water vapor in the soil as well as atmosphere.Therefore, the energy transfer of the desert is dominated by sensible heat flux, while the latent heat flux is very small, which distinguishes the desert from other areas.What is more, land-atmosphere interaction in desert plays an important role in global and regional energy balance and climate change [77].

Conclusions and Further Work
The overland-surface processes in the desert region of northwest China play important roles in the regional weather and climates.α, ε, z 0m , z 0h , C d , and C h are key parameters for the land-air interactions.In this study, these parameters were investigated based on observations on the northern marginal zone of TD and compared with the values estimated from the observations in the central area of TD.The main findings are summarized as follows.
(1) In the northern marginal zone of the TD, α and ε were 0.27 and 0.91, respectively, which are consistent with the values obtained based on observations in the hinterland of the TD as well as being similar to the dry parts of the Great Basin desert in North American.Also, the retrieved α and ε from Landsat 8 data were close to the observed values at both the northern marginal zone and hinterland of the TD.The satellite-derived values of α and ε from remote sensing products are highly reliable for use in the Taklimakan Desert.The values of α varied from sunny to cloudy days, and were 0.31 and 0.30, respectively.Also, the relationship between α and the solar elevation angle (h) during sunny periods and cloudy periods were different, which were α = 0.2586 + 0.2415e −h/8.852 and α = 0.4424h −0.1404 , respectively.The values of ε were sensitive to the observed value of H using the iterative method, that ε decreased (increased) as H obs increased (decreased).
(2) In the marginal zone of the TD, z 0m was dependent on the wind direction.z 0m clearly changed with the wind direction, which could be affected by the spatial inhomogeneity (The distributions of the surrounding sand dunes at the station are not uniform, there were many undulating sand dunes in the prevailing wind and opposite to the prevailing wind).These are similar to the results obtained in vegetated areas where z 0m was highly dependent on the seasonal variations in vegetation conditions.In general, in the NE-E and SSE-WSW directions, there were many undulating sand dunes, which were consistent with the directions of the peak z 0m value, which z 0m was 5.6 × 10 −3 to 6.3 × 10 −3 m in the NE-E wind direction and 5.1 × 10 −3 to 6.1 × 10 −3 m in the S-WSW wind direction.z 0h changed with the wind direction where the magnitudes of the data were 10 −3 and 10 −4 for z 0m and z 0h , respectively.The optimal values for z 0m and z 0h were 5.858 × 10 −3 m and 1.965 × 10 −4 m according to the normally distributed histogram, and the magnitudes of z 0m and z 0h at the station were consistent with the results reported by Stull [63].
(3) z 0m and z 0h varied seasonally.Z 0m was higher in July (5.5 × 10 −3 m) and lower in November (4.2 × 10 −3 m); z 0h was higher in April and September (4.8 × 10 −3 m and 4 × 10 −3 m, respectively) and lower in January and February (3.4 × 10 −3 m and 3.2 × 10 −3 m, respectively).The monthly z 0m varied between 3.2 × 10 −3 (February) and 4.8 × 10 −3 m (April), and the average annual value was 4.78 × 10 −3 m.The monthly z 0h varied between 3.17 × 10 −4 (February) and 4.84 × 10 −4 m (April), and the average annual value was 3.76 × 10 −4 m.The values of ln(z 0m ), ln(z 0h ), and kB −1 exhibited obvious fluctuations at sunrise and sunset.The excess resistance to heat transfer, kB −1 , varied inversely with ln(z 0h ).The daily mean kB −1 was between 0.9 and 3.9, and the average annual value was 2.5.These results were different from those obtained in the vegetated areas with bare soil and smooth surfaces such as the HEIFE Gobi desert, but they were similar to those determined for the Qinghai-Tibetan Plateau and HAPEX-Sahel.
(4) Due to the high H value in this area, C d and C h exhibited large diurnal variations in the marginal zone of TD.The fluctuations in C d and C h were large before sunrise and after sunset, and they were consistent with z 0m .The order was C d > C h in all months except February.The daily mean values of C d and C h were 6.34 × 10 −3 and 5.96 × 10 −3 , respectively, and their neutral values were 4.30 × 10 −3 and 1.70 × 10 −3 .The highest and lowest values for C d were 0.99 × 10 −2 (May) and 0.54 × 10 −2 (March), respectively, while the highest and lowest values for C h were 0.30 × 10 −2 (February) and 0.27 × 10 −2 (January).
Our results also suggested that the parameters derived from remote sensing observations were good quality compared with the in situ observations.Therefore, remote sensing provides a better method for capturing the characteristics of the surface parameters and land-surface processes in a vast desert region where in situ observations are very sparse and difficult to maintain.In future research, we will use our observations to evaluate whether the uncertainty of parameters derived from observations can affect the capacity of land-surface models for modeling land-air interactions in a desert climate.We will also assess how much the uncertainty in surface parameter estimates can affect simulations of the desert weather and climate in this region.

Figure 1 .
Figure 1.(a) Map of Xiaotang station in China and Xinjiang Province.(b) High resolution map of Xiaotang station.(c) Photograph of Xiaotang station in the Taklimakan Desert.(d) Wind direction at Xiaotang station during 5 January to 31 December 2010.(e) Flux footprint (10%, 30%, 50%, and 70%) at the observation station (concentric center), where the labels around the circle show the direction in geographical coordinates (0-360°) and the radius length of the circle is the distance.(f) Same as (e) but for 90%.

Figure 1 .
Figure 1.(a) Map of Xiaotang station in China and Xinjiang Province.(b) High resolution map of Xiaotang station.(c) Photograph of Xiaotang station in the Taklimakan Desert.(d) Wind direction at Xiaotang station during 5 January to 31 December 2010.(e) Flux footprint (10%, 30%, 50%, and 70%) at the observation station (concentric center), where the labels around the circle show the direction in geographical coordinates (0-360 • ) and the radius length of the circle is the distance.(f) Same as (e) but for 90%.

Figure 2 .
Figure 2. (a) Mean daily variations in the energy closure from 5 January to 31 December 2010.(b) Same as (a) but corrected for the soil heat flux.(c) Mean daily variations in the flux (H is the surface sensible heat flux, LE is the latent heat flux, R n is the net radiation, and G 0 is the surface soil heat flux).(d) Same as (c) but corrected for the soil heat flux.(e) Scatter plot with regression for H + LE and R n − G 2.5 .(f) Scatter plot with regression for H + LE and R n − G 0 .(g) Mean daily variations in the soil heat flux at depths of 2.5 cm and 8 cm.(h) Mean daily variations in the friction velocity and Bowen ratio.

Figure 3 .
Figure 3. (a) Observed and calculated functions of φm against / .(b) Same as (a) but for φh.

Figure
Figure4bshows the relationship between α and h during cloudy periods by the following numerical value Equation(21).

Figure 3 .
Figure 3. (a) Observed and calculated functions of ϕ m against z/L.(b) Same as (a) but for ϕ h .

Figure 4 .
Figure 4. (a) Changes in the surface albedo with solar elevation angles on sunny days at Xiaotang station during 5 January to 31 December 2010.(b) Same as (a) but on cloudy days.(c) Scatter plot of outgoing against incoming radiation on sunny days.(d) Same as (c) but on cloudy days.(e) Diurnal

Figure 4 .
Figure 4. (a) Changes in the surface albedo with solar elevation angles on sunny days at Xiaotang station during 5 January to 31 December 2010.(b) Same as (a) but on cloudy days.(c) Scatter plot of outgoing against incoming radiation on sunny days.(d) Same as (c) but on cloudy days.(e) Diurnal variations in the surface albedo on sunny days and cloudy days.(f) Retrieval of the surface albedo from Landsat data.

Figure 5 .
Figure 5. (a) Dependence of the surface emissivity ( ε ) on the RMS in Hobs and Hsfc in near neutral conditions at Xiaotang from 5 January to 31 December 2010.(b) Dependence of the surface emissivity ( ε ) on density according to the surface temperature method.(c) Retrieval of the surface emissivity from Landsat data.

Figure 5 .
Figure 5. (a) Dependence of the surface emissivity (ε) on the RMS in H obs and H sfc in near neutral conditions at Xiaotang from 5 January to 31 December 2010.(b) Dependence of the surface emissivity (ε) on density according to the surface temperature method.(c) Retrieval of the surface emissivity from Landsat data.

Figure 6 .
Figure 6.(a) Frequency distribution of the estimated ln(z0m) parametric density at Xiaotang station from 5 January to 31 December 2010.The optimal value of ln(z0m) was −5.14, which has the highest frequency in the curve.(b) Same as (a) but for ln(z0h) where ln(z0h) was −8.54.

Figure 6 .
Figure 6.(a) Frequency distribution of the estimated ln(z 0m ) parametric density at Xiaotang station from 5 January to 31 December 2010.The optimal value of ln(z 0m ) was −5.14, which has the highest frequency in the curve.(b) Same as (a) but for ln(z 0h ) where ln(z 0h ) was −8.54.

Figure 7 .
Figure 7. (a) Changes in z0m with different wind directions at Xiaotang from 5 January to 31 December 2010.(b) Same as (a) but for z0h.(c) Same as (a) but for a magnitude of 10 −3 .(d) Same as (b) but for a magnitude of 10 −4 .

Figure 7 .
Figure 7. (a) Changes in z 0m with different wind directions at Xiaotang from 5 January to 31 December 2010.(b) Same as (a) but for z 0h .(c) Same as (a) but for a magnitude of 10 −3 .(d) Same as (b) but for a magnitude of 10 −4 .

Figure 9 .
Figure 9. (a) Monthly variations in the mean values of ln(z0m), ln(z0h), and kB −1 from 5 January to 31 December 2010.(b) Same as (a) but for the surface-air temperature.(c) Same as (a) but for the wind speed.(d) Same as (a) but for the sensible heat and latent heat fluxes.

Figure 9 .
Figure 9. (a) Monthly variations in the mean values of ln(z0m), ln(z0h), and kB −1 from 5 January to 31 December 2010.(b) Same as (a) but for the surface-air temperature.(c) Same as (a) but for the wind speed.(d) Same as (a) but for the sensible heat and latent heat fluxes.

Figure 9 .
Figure 9. (a) Monthly variations in the mean values of ln(z 0m ), ln(z 0h ), and kB −1 from 5 January to 31 December 2010.(b) Same as (a) but for the surface-air temperature.(c) Same as (a) but for the wind speed.(d) Same as (a) but for the sensible heat and latent heat fluxes.

Figure 9 . 26 Figure 10 .
Figure 9. (a) Monthly variations in the mean values of ln(z0m), ln(z0h), and kB −1 from 5 January to 31 December 2010.(b) Same as (a) but for the surface-air temperature.(c) Same as (a) but for the wind speed.(d) Same as (a) but for the sensible heat and latent heat fluxes.

Figure 10 .
Figure 10.(a) Mean daily variation in C d , with wind directions of NNE-ESE from 5 January to 31 December 2010.(b) Mean daily variations in C h , with wind directions of NNE-ESE from 5 January to 31 December 2010.(c) Mean daily variations in z/L with wind directions of NNE-ESE from 5 January to 31 December 2010.(d) Mean daily variations in u/u * with wind directions of NNE-ESE from 5 January to 31 December 2010.

Figure
Figure 10a-d shows the variations in C d and C h according to their relationships with z/L and u/u * .Figures8 and 10a,b shows that the changes in z 0m and z 0h were higher after sunrise and sunset, where they were consistent with the variations in C d and C h .Figure10a,bshows that the values of C d and C h were higher (lower) in the daytime (at night).According to the strong diurnal variations in the surface temperature, z/L < 0 after sunrise, but it was stable (z/L > 0) after sunset.Thus, u/u * was lower (higher) in the daytime (at night) (Figure10d).As a result, C d and C h also exhibited diurnal variations,

Figure 11 .
Figure 11.(a) Monthly variations in the mean values of Cd under NNE-ESE wind direction from 5 January to 31 December 2010.(b) Monthly variations in the mean values of Ch under NNE-ESE wind direction from 5 January to 31 December 2010.

Figure 11 .
Figure 11.(a) Monthly variations in the mean values of C d under NNE-ESE wind direction from 5 January to 31 December 2010.(b) Monthly variations in the mean values of C h under NNE-ESE wind direction from 5 January to 31 December 2010.

Figure 12 .
Figure 12.(a) Relationship between the C d and wind speed from 5 January to 31 December 2010.(b) Relationship between the C h and wind speed from 5 January to 31 December 2010.

Figure 26 Figure 12 .
Figure 13a,b also shows the relationship estimated between C d and C h using the Eddy correlation method and the atmospheric stability in the NNE-ESE wind direction.For comparison, Figure 13c,d shows the results based on the estimates obtained using the M-O similarity function (Figure 3c is unrevised and Figure 3d is revised).As shown in Figure 13, C d varied from 2 × 10 −3 to 10 × 10 −3 and C h from 0-9 × 10 −3 when z/L ranged from −3 to 0.5 at the station.C d and C h increased as z/L decreased (z/L ≤ 0), which is consistent with the results of previous studies in mid-latitudinal regions [4,70,75].The variations in C d and C h versus the wind speed under different values of z/L are shown in Figure 14.C d and C h decreased rapidly as the wind speed increased under weak wind conditions (≤3.0 m s −1 ) and then remained almost constant as the wind speed increased further, where the minimum value was around 2 m s −1 .By contrast, with the revised M-O similarity function, the difference was very small for C h but large for C d .

Figure
Figure 13a-b also shows the relationship estimated between Cd and Ch using the Eddy correlation method and the atmospheric stability in the NNE-ESE wind direction.For comparison, Figure 13c-d shows the results based on the estimates obtained using the M-O similarity function (Figure 3c is unrevised and Figure 3d is revised).As shown in Figure 13, Cd varied from 2 × 10 −3 to 10 × 10 −3 and Ch from 0-9 × 10 −3 when z/L ranged from −3 to 0.5 at the station.Cd and Ch increased as z/L decreased (z/L ≤ 0), which is consistent with the results of previous studies in mid-latitudinal regions [4,70,75].The variations in Cd and Ch versus the wind speed under different values of z/L are shown in Figure 14.Cd and Ch decreased rapidly as the wind speed increased under weak wind conditions (≤3.0 m s −1 ) and then remained almost constant as the wind speed increased further, where the minimum value was around 2 m s −1 .By contrast, with the revised M-O similarity function, the difference was very small for Ch but large for Cd.

Figure 13 .
Figure 13.(a) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the Eddy correlation method.(b) Same as (a) but with wind directions of NNE-ESE.(c) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the unrevised M-O similarity function.(d) Same as (c) but revised.

Figure 13 .
Figure 13.(a) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the Eddy correlation method.(b) Same as (a) but with wind directions of NNE-ESE.(c) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the unrevised M-O similarity function.(d) Same as (c) but revised.

Figure 13 .
Figure 13.(a) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the Eddy correlation method.(b) Same as (a) but with wind directions of NNE-ESE.(c) Relationship between the bulk transfer coefficients and atmospheric stability with wind directions of NNE-ESE from 5 January to 31 December 2010 according to the unrevised M-O similarity function.(d) Same as (c) but revised.

Figure 14 .Figure 14 .
Figure 14.(a) Relationship between the bulk transfer coefficient Cd and wind speed with wind directions of NNE-SES from 5 January to 31 December 2010.(b) Same as (a) but for Ch.(c) Same as (a) Figure 14.(a) Relationship between the bulk transfer coefficient C d and wind speed with wind directions of NNE-SES from 5 January to 31 December 2010.(b) Same as (a) but for C h .(c) Same as (a) but according to the revised M-O similarity function.(d) Same as (b) but according to the revised M-O similarity function.Table 6 shows the values of C d and C h calculated using the Eddy correlation method and the M-O similarity function (with NNE-ESE wind) method.In near natural conditions, the mean values of C d (C h ) were 20.72 × 10 −3 , 4.30 × 10 −3 , 4.23 × 10 −3 , and 4.44 × 10 −3 (6.62 × 10 −3 , 1.70 × 10 −3 , 2.74 × 10 −3 , and 1.55 × 10 −3 ) with the Eddy correlation method (all wind directions), Eddy correlation method (with NNE-ESE wind), unrevised M-O similarity function (with NNE-ESE wind), and revised M-O similarity function (NNE-ESE), respectively.For the results derived based on all wind conditions, the mean values of C d (C h ) were 20.31 × 10 −3 , 6.34 × 10 −3 , 12.14 × 10 −3 , and 7.94 × 10 −3 (67.13 × 10 −3 , 5.96 × 10 −3 , 3.81 × 10 −3 , and 3.33 × 10 −3 ), respectively.The values of C d (C h ) at the station calculated using the Eddy correlation method (all wind direction) were larger than those calculated with the other three methods, i.e., Eddy correlation method (with NNE-ESE wind direction), and unrevised and revised M-O similarity function (with NNE-ESE wind direction).The values of C d (C h ) at the station calculated with the unrevised M-O similarity function (NNE-ESE) were larger than those calculated with the revised M-O similarity function (NNE-ESE).

Table 1 .
Measurement levels and instruments employed for acquiring near-surface meteorological variables and turbulent fluxes at Xiaotang station.

Table 2 .
EddyPro parameters used in this study.

Table 3 .
EBR in every month.

Table 4 .
Difference in ε with Hobs according to the method employed.

Table 4 .
Difference in ε with H obs according to the method employed.

Table 5 .
Estimates of ε according to the temperature gradients in near zero conditions.

Table 6 .
Values of C d and C h calculated using different methods.