## Abstract

Numerical weather and climate prediction systems necessitate accurate land surface–atmosphere fluxes, whose determination typically replies on a suite of parameterization schemes. The authors present a field investigation over tall grass in a Beijing suburb, where the aerodynamic roughness length (*z*_{0}) and zero-plane displacement height (*d*) are found to be 0.02 and 0.44 m, respectively (the value of *d* is close to two-thirds the average grass height during this field experiment). Both *z*_{0} and *d* values are then used as input parameters of an analytic model of flux footprint; the footprint model reveals that eddy-covariance flux measurements are mainly representative of the tall grass surface concerned herein, potential influences from anthropogenic sources in this suburban area notwithstanding. Based on the “fair weather” data (with an energy balance ratio of 0.89), the authors evaluate four parameterizations of turbulent surface fluxes, namely, a total of three traditional “iterative” schemes and one “noniterative” scheme developed recently to reduce computational time. Their performances are intercompared in terms of the estimations of the sensible heat flux and two turbulence components (the friction velocity and temperature scale). In weakly stable to unstable conditions, two schemes are recommended here for their good performance overall; the first scheme stems jointly from the work of Högström and Beljaars and Holtslag, and the second stems from that of Li et al.. For this tall grass surface, the choice of *z*_{0}/*z*_{0h} = 100 (*z*_{0h} is the thermal roughness length) is more appropriate than another choice of 10, because the former produces comparatively accurate sensible heat flux estimations.

## 1. Introduction

Numerical weather and climate prediction systems depend on the accurate partitioning of available energy at the surface between turbulent sensible heat and latent heat fluxes (Sellers et al. 1996; Chen and Dudhia 2001; Holt et al. 2006; Ban et al. 2010). The partitioning patterns may differ at individual grassland sites located in different climate zones. For example, Bi et al. (2007) measured the energy partitioning over grassland in the tropical monsoon environment of southern China and found that the latent heat flux was a primary consumer of available energy in summer. Bian et al. (2003) reported that sensible heat accounted for about 80% of the available energy during a dry period at the Qamdo site on the Tibetan Plateau, while the latent heat flux was almost double the sensible heat flux over a prairie in core regions of the Tibetan Plateau during the monsoon season (Gao et al. 2004). Lately, Gao et al. (2009) found that the sensible heat flux was a primary consumer of available energy for the entire year by examining energy partitioning over a steppe prairie in Inner Mongolia, northern China. The different patterns in grassland surface energy partitioning are caused by surface properties such as soil moisture and vegetation cover. So far, to our knowledge, there has been no reporting of energy partitioning over the grasslands around the city of Beijing, although strong evaporation in this area tends to trigger or enhance local convection systems that usually bring precipitation to Beijing. To understand the turbulent transfer processes in a Beijing suburb, a field experiment was conducted at a suburb within the Beijing–Tianjin–Tanggu metropolitan region in northern China from 19 July through 6 August 2010. The first objective of this work was to characterize the radiation budget and surface energy balance by using micrometeorological measurements over a tall grass site in this suburban area.

Another issue addressed in our work concerned the evaluation of several often adopted schemes of surface flux parameterization, which are constructed within the framework of bulk parameterization. It is of interest to compare their performance in terms of turbulent flux estimations, based on the available measurements over tall grass. Specifically, turbulent fluxes on the land surface–atmosphere interface are modulated by a variety of meteorological parameters, such as the mean wind speed, near-surface gradients in air temperature and other scalar quantities concerned, atmospheric thermal stratification (static stability), and aerodynamic and thermal roughness lengths. Extensive efforts have been devoted to develop parameterization schemes of surface turbulent fluxes for numerical weather and climate prediction systems. In early versions of related numerical models, the turbulent transfer coefficients were derived from the Monin–Obukhov similarity theory (Monin and Obukhov 1954), where an iterative algorithm was usually necessary to solve surface wind stress and heat/moisture fluxes by consulting the bulk parameterization approach (e.g., Businger et al. 1971; Dyer 1974; Holtslag and de Bruin 1988; Beljaars and Holtslag 1991; Högström 1996). Because computationally intensive iterations were needed for all grid nodes at the lowest level of the atmospheric model for all time steps, a huge amount of computer power was consumed, especially for a high-resolution simulation over regional or global scales (see Li et al. 2010 for reference). Louis (1979) and Louis et al. (1982) proposed a compact scheme where the bulk Richardson number (Ri_{B}) instead of the Obukhov length (*l*) was used for the parameterization of the momentum and sensible heat fluxes; the scheme is more efficient because values of Ri_{B} are readily available from routine weather station observations. Since then, several different noniterative parameterization schemes have been developed. Recently, Li et al. (2010) proposed new formulations to calculate turbulent transfer coefficients, which were constructed on the basis of classic iterative computational results given by Launiainen (1995).

As far as practical atmospheric modeling is concerned, two traditional schemes have been most commonly adopted, which are developed by Businger et al. (1971) and Dyer (1974) and are hereafter referred to as the Businger scheme and the Dyer scheme, respectively. Another set of schemes has also received much attention, namely, the schemes proposed by Beljaars and Holtslag (1991) for stable conditions and by Högström (1996) for unstable conditions (see Guo and Zhang 2007; Guo et al. 2007 for related evaluations); their schemes are collectively referred to herein as the BHH scheme. The above three schemes are computationally intensive because of the involved numerical iterations. Moreover, these traditional schemes (namely the Businger, Dyer, and BHH schemes) were selected as references for comparison against a recently developed noniterative algorithm (Li et al. 2010; referred to as the new scheme herein). Therefore, in order to provide data support for turbulent flux parameterization concerned with climate and weather models, the second objective of this work is to verify whether the three traditional “iterative” schemes and one “noniterative” scheme lead to reliable flux estimations over the tall grass surface, through comparing the estimated results against eddy-covariance flux data collected during the field experiment.

Given the fact that various schemes are originally proposed for use over relatively flat and homogeneous surfaces (e.g., with uniform vegetation), it has long been of interest to evaluate surface flux parameterizations over heterogeneous surfaces characterized, for example, by a mixture of land covers. One primary component of our investigation, therefore, is to evaluate several of the most commonly adopted ones (the Businger, Dyer, BHH, and new schemes), based on micrometeorological measurements over a tall grass site. An interesting aspect lies in the fact that related evaluations previously made elsewhere are often targeted at bare land or sparsely vegetated surfaces, while our tall grass site also involves the determination of zero-plane displacement height. Conclusions derived from this particular vegetated site are supposed not only to guide the choice of a comparatively accurate scheme, but also to ascertain the general usefulness of the candidate schemes over fairly high, complex vegetation (though their accuracy in the flux estimations is found to vary by a significant degree).

## 2. Materials and methods

### a. Site

Micrometeorological measurements were conducted at a tall grass site (40.1627°N, 116.1121°E, 60 m above sea level) about one kilometer west of Beijing from 19 July through 6 August 2010. The grass was about 0.65 m in height. The prevailing wind was from the southeast and northwest in the summer season (see Fig. 1). The upwind fetch is comparatively uniform to the southeast of the site, as shown in Fig. 1.

### b. Micrometeorological measurements and meteorological conditions

Two towers were set up to mount the meteorological instruments. One was 12 m high and the other was 3 m. On the 12-m tower, meteorological instruments included an eddy-covariance system using a CSAT3 three-dimensional sonic anemometer (Campbell Scientific Inc., Logan, Utah, United States) and a LI-7500 open-path gas analyzer (Li-Cor Inc., Lincoln, Nebraska, United States). They measured three components of turbulent wind velocity, virtual temperature, and water vapor and CO_{2} concentrations; the turbulence raw data were collected at 10 Hz. Also available were CNR1 pyronometers and pyrgeometers (Kipp & Zonen, Delft, the Netherlands), which measured four radiation components. The eddy-covariance (CSAT3 and LI-7500) and CNR1 instruments were mounted at a height of 10.0 m relative to the surface. Additionally, a total of three HMP45C temperature and relative humidity sensors (housed in radiation shields; Campbell Scientific Inc.) and three 010C-1 wind speed sensors (Met One Instruments Inc., Grants Pass, Oregon, United States) made auxiliary measurements of the ambient air, namely, at the heights of 4.0, 7.5, and 11.0 m above the surface. Installed respectively at the heights of 12.0 and 2.0 m, 020C-1 wind direction (Met One Instruments Inc.) and CS106 barometric pressure (Campbell Scientific Inc.) sensors were used to measure the mean wind direction and air pressure. On the 3-m tower, the CSAT3 and LI-7500 instruments were both installed at a height of 2.5 m above the surface, and CNR1 pyronometers and pyrgeometers were installed at a height of 2.0 m. These two towers were separated by about 5 m in the horizontal direction, so that the instruments mounted on them could well represent the tall grass concerned herein; notice that there was no discernible surface slope at this site.

Soil temperature was measured using four TCAV averaging soil thermocouple probes (Campbell Scientific Inc.) at 0.05-, 0.10-, 0.20-, and 0.40-m depths. The volumetric water content of the soil was measured at 0.10- and 0.20-m depths using two soil moisture reflectometers (CS616, Campbell Scientific Inc.). The soil heat flux was measured at 0.05- and 0.10-m depths using two heat flux plates (HFP01SC, Campbell Scientific Inc.,). The soil measurements were made at a distance of 2 m away from the 10-m tower.

Among the aforementioned measurements, various sensors were used to provide half-hourly average values (except for the eddy-covariance instruments). Besides, half-hourly rainfall amount was monitored using a TE525MM tipping-bucket rain gauge (Campbell Scientific Inc., Logan).

Contemporary flux measurements made over various ecosystems often adopt 10 or 20 Hz as the sampling frequency of the eddy-covariance system. As a common practice among the flux community, in order to determine whether the adopted sampling frequency is appropriate for deriving eddy-covariance turbulent fluxes, ogive functions are often calculated to get relevant information. We, accordingly, have calculated the ogive functions for both the momentum and sensible heat fluxes. The results enable us to confirm that 10 Hz is an appropriate sampling frequency over the tall grass site concerned herein, for both measurement heights of 2.5 and 10.0 m (figure omitted).

Figure 2 shows the time series of meteorological variables: wind speed, wind direction, air temperature, relative humidity, atmospheric pressure, and precipitation during the entire experiment. The first four variables were measured at the 11-m height on the 12-m tower. It can be seen that the wind speed was generally weak and varied between about 1.0 m s^{−1} northwesterly at night and 2–5 m s^{−1} southeasterly during daytime. A clear diurnal signal was present in the wind speed, with the maximum wind speed occurring around noon. The minimum and maximum air temperatures were about 19° and 37°C, respectively. Highest humidity was concurrent with low-pressure and precipitation episodes. The diurnal cycles of the radiative components, except for downward longwave, were also clear during a large portion of the study period; downwelling shortwave radiation (DSR), upwelling shortwave radiation (USR), downwelling incoming longwave radiation (DLR), and upwelling longwave radiation (ULR) are shown in Fig. 3. The albedo of the underlying surface, defined as the ratio of the daily peak values of upward and downward shortwave radiation, is about 0.13 averaged over the entire observational period.

### c. Theory

The surface energy balance over a grass canopy can be approximated by

where is the net radiation (W m^{−2}); and are the sensible heat and latent heat fluxes, respectively; is the surface soil heat flux (W m^{−2}); and is the residual energy involved in various processes, such as photosynthesis and respiration. The ratio of the sum of sensible heat and latent heat fluxes () to available energy [the difference of net radiation and surface soil heat flux ()] is used to determine the surface energy balance ratio. The net radiation is given by .

The heat flux at the soil surface is estimated by using a combination of soil calorimetry and measurement of soil heat flux () at a certain depth (0.05 m in the present study) in the soil using heat flow transducers. The heat storage [*C _{g}δz* (

*δT*/

*δt*)] of the soil layer above the heat flux plate is included as

where is the volumetric heat capacity of soil (J m^{−3} K^{−1}), *δz* is the soil thickness (0.05 m in this paper), is the soil temperature (K), and is the time (s). The volumetric heat capacity was estimated as follows (De Vries 1963): , where is the volumetric water content; is the saturated value of ; and and are the volumetric heat capacities of dry soil and water, respectively.

The eddy-covariance method was used to calculate fluxes of momentum (*τ*), sensible heat, latent heat, and carbon dioxide (*F*_{CO2}) using values of instantaneous fluctuations observed by the instrument. These fluxes can be expressed as

where is the mean air density (kg m^{−3}); is the specific heat capacity of air (1005.0 J kg^{−1} K^{−1}); is the latent heat of vaporization (2.5 × 10^{6} J kg^{−1}); and *u*′ (*υ*′ as well), , , , and *c*′ are the fluctuations in the horizontal and vertical wind components (m s^{−1}), potential temperature (K), specific humidity (kg kg^{−1}), and concentration of CO_{2} (mg m^{−3}), respectively.

Based on single-level sonic anemometer data, Martano (2000) proposed a method for estimating the zero-plane displacement height and the aerodynamic roughness length by using the least squares method. Footprint analysis is commonly used to quantify the contributing source areas to the scalar flux measurement (Schmid 1994); it is worth recalling that analytic footprint models are based on ideal conditions, such as steady airflow and optimal surface homogeneity, though they are extensively adopted in the literature, with less regard to the surface conditions encountered. To estimate the average footprint for the entire experiment, the cumulative function of flux density is expressed as

where is the upwind distance (m) from the measurement site, the von Kármán constant (Andreas et al. 2006), is the zero-plane displacement height (m), *u*_{*} is the friction velocity (m s^{−1}), and is the average wind speed (m s^{−1}) between the displacement height and the observation height (m) above the surface. Assuming a logarithmic profile (without effects of stratification) for horizontal wind speed , the average wind speed can be computed by

where *z*_{0} is the aerodynamic roughness length (m) and for brevity hereafter, *z*′ stands for *z* − *d*.

Following the Monin–Obukhov similarity theory (Monin and Obukhov 1954), the friction velocity () and temperature scale () at the lowest level (, as concerned herein) of atmospheric models can be obtained by integrating the flux-profile equations:

where (K) is the potential temperature; *l* (m) is the Obukhov length; (K) is the surface potential temperature defined at the height of (m), which stands for the thermal roughness length; and and are the integrated stability correction functions for momentum and heat, respectively. The momentum () and heat () fluxes are calculated from the bulk formulae

where is the drag coefficient and is the heat transfer coefficient. It is obvious that

and, for instance, Businger et al. (1971), Dyer (1974), and Högström (1996) gave the functional expressions of both and . Holtslag and de Bruin (1988) and Beljaars and Holtslag (1991) were concerned specifically with their formulations under stable stratification conditions. For brevity, their equations are not given here. Following the ideas of Louis's compact scheme (Louis 1979; Louis et al. 1982), Li et al. (2010) regressed the relationship between the nondimensional stability parameter and the bulk Richardson number (*g* is the gravitational acceleration, m s^{−2}) within 5% error and obtained the following.

When ,

where , , , , , , , , , and .

When ,

where , , , , , , , and .

When ,

where , , , , and .

The statistical bias (BIAS), mean absolute error (MAE), and root-mean-square error (RMSE) that will be used in this study are given by

where is the result modeled using the parameterization schemes, is the result measured using the eddy-covariance method, and is the total number of data points.

## 3. Results and discussion

### a. Determination of the aerodynamic roughness length and zero-plane displacement height

The aerodynamic roughness length is practically sensitive to friction velocity (Johnson et al. 1998). To reduce the uncertainty and make the computations using the Martano (2000) method converge faster, we selected data with the ratio of *u*_{*} to *U* ranging from 0.053 to 0.083, and under moderate wind conditions (1.5 m s^{−1 } 3.5 m s^{−1}). The standard deviation of the quantity , , exhibited a clear minimum with m, as shown in Fig. 4a. Accordingly, was readily found by applying the minimum deviation principle. Figure 4b shows (the standard deviation of wind speed *U*) was minimized with m. The determined values of and fall into their respective ranges as expected; explicitly, is proposed to be about two thirds the plant canopy height over vegetated surfaces, and values typically vary on the order of centimeters (i.e., one order of magnitude higher than those over sparse vegetation or bare land) (cf. Garratt 1992; Stewart and Oke 2012).

### b. Eddy-covariance flux corrections and footprint analysis

To ensure the accuracy of eddy-covariance turbulent fluxes, several steps are followed for the sake of necessary corrections and data screening (cf. Mauder and Foken 2006; Aubinet et al. 2012). We use 30-min bulk averaging throughout this study. To calculate sensible heat fluxes, we follow Schotanus et al. (1983) to convert the sonic temperature into actual temperature. Eddy-covariance latent heat and CO_{2} fluxes involve a Webb–Pearman–Leuning (WPL) correction for density effects (cf. Webb et al. 1980; see Lee and Massman 2011 for a historical revisit). For corrections of potential flux loss at high frequencies (possibly caused by path averaging and limited sampling frequency; see Kaimal and Finnigan 1994), raw flux data are corrected using the transfer-function technique with the cospectral models of Kaimal et al. (1972) (see Moore 1986 for algorithm). The planar fit coordinate is then applied by following the recommendation of Wilczak et al. (2001). The above corrections are implemented by using the EdiRe software, developed by the University of Edinburgh (see http://www.geos.ed.ac.uk/abs/research/micromet/EdiRe/). Following Foken and Wichura (1996) and Foken et al. (2004), a stationarity test is made for the corrected fluxes. A total of 275 half-hourly eddy-covariance turbulence measurements are used for further analysis; notice that the data during rainfall periods are rejected on the basis of the LI-7500 automatic gain control (AGC) values.

Defensible efforts to evaluate surface flux parameterizations are demanding, not least because of the importance to ensure sufficient spatial representativeness of the eddy-covariance fluxes (used in our work as “benchmark” data against which the parameterized fluxes are compared). With m, m, m, and the mean m s^{−1}, the average wind speed (*U*) was calculated to be 3.2 m s^{−1} using Eq. (8); then, the cumulative function of flux density was computed from Eq. (7). Related results are presented in Fig. 5, where upwind distances with the highest crosswind contributions are found to be about 157.0 and 23.5 m for measurements at the 2.5- and 10.0-m heights, respectively. Given the geographic information in Fig. 1, we confirm that primary source areas of the fluxes reside well within the tall grass surface that had a fetch longer than 300 m in all directions. Nevertheless, we are unable to exclude potential influences exerted by land covers other than the targeted grassland (such as several adjacent villages) because cumulative footprint contributions in Fig. 5 reveal that 90% of the measured fluxes stem from upwind contributions that extend for as far as 3000 and 450 m or so for the 10.0- and 2.5-m heights, respectively. Therefore, our eddy-covariance flux measurements at the 2.5-m height are deemed well representative of the targeted grassland, but those at the 10.0-m height probably involve nonnegligible influence that should arise from a mix of land covers beyond the grassland boundaries (surface heterogeneity is characteristic of urbanized landscapes in the Beijing suburb). It is worth noting that our footprint calculations are based merely on data in neutral conditions, but flux source areas vary with the atmospheric stability, which extend in spatial coverage particularly in stable conditions (cf. Schmid 1997; Vesala et al. 2008); stable conditions occur for about 36% of the time within our dataset, when the flux data may be less representative of the targeted tall grass surface.

### c. Diurnal variations of energy components and CO_{2} flux

Figure 6 shows the diurnal variations of , , , and ; and were measured by the fast response instruments and calculated respectively using Eqs. (4) and (5), and was estimated by using Eq. (2). Both and , measured at two different towers (10-m level for the 12-m tower and 2.5-m level for the 3-m tower), demonstrate that these two instruments give similar variation of surface heat exchange and consistent values. The coefficient of determination *r*^{2}, RMSE, and BIAS between the fluxes measured by these two instruments are shown in Table 1. Based on the comparative statistics given in Table 1, there is an overall agreement between the flux measurements at the 2.5-m and 10.0-m heights, particularly for the sensible heat () and latent heat () fluxes, as evidenced by the relatively high *r*^{2} values derived, that is, 0.95 and 0.83 for the *H* and *LE* comparisons, respectively; by contrast, the counterpart for the *F*_{CO2} comparison (0.68) is somewhat lower because of more data scatter (figure omitted). Such an agreement is indicative of the likely existence of a “constant flux” surface layer over the present tall grass site. Despite the surface patchwork seen on the kilometer scale (see Fig. 1), our flux dataset should be regarded to involve insignificant influences by scalar sources beyond the targeted grassland. To some extent, the agreement, as noted above, should also indicate the appropriate adoption of related flux correction approaches over this vegetated site (see section 3b).

Temporal evolution in the flux measurements was also evident during this field experiment. For instance, relatively high evaporation following 17–19 July precipitation reduced the soil volumetric water content gradually (not shown) and the value of was often larger than the value of during the drying period from 20 to 21 July. After that, daytime magnitude of gradually decreased and reached about 100 W m^{−2}, and daytime magnitude of tended to increase and reached about 200 W m^{−2} during the clear days from 22 to 25 July. In terms of the magnitude of both heat flux components, they are largely similar to the measurements previously carried out over grassland ecosystems for a wide range of elevations (e.g., Laubach and Teichmann 1999; Gao et al. 2006; Bian et al. 2012; Liu and Feng 2012). Moreover, the sensible heat and latent heat fluxes lead to Bowen ratio values that mostly vary between 0.1 and 2.0, which suggest that the flux measurements represent intermediate and dry hydrological conditions (e.g., Choi et al. 2004; Moene and Schüttemeyer 2008; Guo et al. 2009), probably as a result of the lack of precipitation preceding the experiment.

The maximum value of was about 660 W m^{−2} on 20 July, and consumed 360 W m^{−2}, that is, about 55% of , at midday peak. The second largest consumer was , which accounted for about 27% (180 W m^{−2}) of the net radiation (), while represented only 15% of (95 W m^{−2}), and the remaining Re was about 25 W m^{−2}. During dry periods, and dominated daytime energy partitioning. At nighttime, was overwhelmingly dominant.

Similar to Fig. 6, the daily cycles of CO_{2} fluxes (*F*_{CO2}; mg m^{−2} s^{−1}) measured at the two heights are shown in Fig. 7. The CO_{2} flux was in phase with with midday uptake rates peaking at over 0.6 mg m^{−2} s^{−1} (absolute values) for these two levels. The daytime CO_{2} absorption after rainy events (17–19 July) was significant because of the strong photosynthesis of the plants (concurrent with relatively high values of the net radiation, ). The daytime CO_{2} flux (absolute values) then dropped to less than 0.3 mg m^{−2} s^{−1}, with slight increases following a precipitation event at midnight on 4 August. A noteworthy comparison concerns the phenomenology of CO_{2} exchange during the daytime. Over urbanized landscapes, the surfaces are known to be a net source of CO_{2}, mainly as a result of intensive anthropogenic emissions (cf. Oke 1987; Hiller et al. 2011; Nordbo et al. 2012). The urban core area of Beijing also appears as a strong source of CO_{2} emissions, based on long-term flux measurements (see Liu et al. 2012). By contrast, our CO_{2} flux data during the daytime hours indicate a sink of carbon over the tall grass ecosystem. During the nighttime, the wind directions were often from the northwest-west, and the flux measurements represented a larger area corresponding to a more extensive footprint. The average *F*_{CO2} was about 0.18 mg m^{−2} s^{−1}, which indicated a surface source of CO_{2} at night. This positive flux could result from either the soil respiration over the grassland or the anthropogenic influences from the adjacent villages.

### d. Energy balance closure

The intercomparison of and is given in Fig. 8. In this analysis, we focus on the data collected on six clear-sky days (i.e., from 21 through 24 July and 31 July through 1 August) because the sonic anemometer would be in malfunction during and after rain events. Then, 285 points were chosen for the analysis of the surface energy balance ratio (). The value of was found to be 0.89, and the correlation coefficient between and was 0.95 (see Fig. 8). Notice that a lower value of was derived from the entire dataset (figure omitted); therefore, we regard it better to select the “fair weather” measurements for further analysis.

In the literature, the lack of a perfect surface energy balance closure has been attributed to various factors, with primary ones including the omission of energy storage components, instrumental problems and/or flux measurement errors caused, for example, by flow distortion, limited high-frequency response, sonic anemometer attack angle, and flux averaging period; the omitted flux contributions in the form of horizontal/vertical advections are often quoted as possible factors, especially over complex terrain (e.g., Wilson et al. 2002; Mauder et al. 2007; Foken 2008; Mauder et al. 2010; Kidston et al. 2010; Kilinc et al. 2012). Following the regression analysis in Fig. 8, the energy balance ratio is estimated as 89% for this tall grass site. This value falls within the wide range reported in previous investigations, and, compared to the mean imbalance ratio on the order of 20% (see Wilson et al. 2002), the outcome of our energy balance assessment herein seems to be acceptable. Because commonly used flux corrections have been implemented in our dataset, the most possible cause of the observed imbalance may be the omission of potential horizontal advection on the spatial scale of order one kilometer over the present suburban site.

### e. Turbulent flux parameterization

Based on Eqs. (9), (10), and (12), , , and were modeled using the four schemes and compared against the measurements (Fig. 9). The linear regressions with 95% confidence interval were carried out for these three variables. In Fig. 9 (left column), it is observed that all four candidate schemes tended to underestimate the friction velocity , especially the Businger scheme. The negative value of BIAS (see Table 2) also indicates this underestimation, and the error metrics of these four schemes showed that the differences produced by the Businger scheme were more evident than those produced by the other schemes. The results also revealed that these four schemes generated unacceptable errors in very stable conditions ( > 0.2), which may partly stem from the fact that most parameterization schemes were developed mostly for unstable conditions. Figure 9 (middle column) compares the temperature scale modeled by the schemes against measurements. In contrast to the results of , the regressions show that the traditional Businger and Dyer schemes obviously tended to overestimate , while the BHH and new noniterative schemes reproduced well, and the new one gave a very good agreement with the slope approaching 1.0. The BHH and the new schemes also resulted in smaller error metrics, as shown in Table 2. Similar to , under strongly stable conditions, all the schemes failed to simulate . Finally, the modeled [calculated using Eq. (12)] showed that the new scheme gave large correlation coefficients (0.89), the best agreement (slope = 1), and a BIAS of only 2.8 W m^{−2} (see Fig. 9, right column).

The drag coefficient and the heat transfer coefficient were calculated respectively using Eqs. (13) and (14). Most of the values of and scattered around 10^{−3}. The value of was determined as 0.02 m, so 500 for the 10-m height. However, it is very difficult to calculate the thermal roughness length accurately. As is well known, the drag coefficient for momentum and the bulk transfer coefficient for heat are quite dependent upon the value of the ratio mostly in the range of 1–100 (Garratt 1992; Mascart et al. 1995; Li et al. 2010). Therefore, over the grassland, and were used in Eqs. (13) and (14), and ranged from −2 to 1. We chose the BHH scheme and the newly developed noniterative scheme to estimate the bulk coefficients. The relationships between and calculated by different schemes and the bulk Richardson number are shown in Fig. 10. It is obvious that these two schemes gave very similar results and could capture the relationship well between and for both and . For modeling the heat transfer coefficient , assuming can lead to significant errors in unstable conditions, and when , results were much more realistic.

A multitude of previous investigations has dealt with the evaluation of surface flux parameterization schemes in the literature. As far as the frequently adopted approach of bulk parameterization is concerned, most previous efforts are devoted to specific formulations within the underlying theoretical framework, namely, the Monin–Obukhov similarity theory (cf. Cassano et al. 2001; Andreas 2002); for example, several works tackled the formulations of flux-profile relationships (e.g., Guo and Zhang 2007; Guo et al. 2007), and others were concerned with those of aerodynamic or thermal surface roughness lengths (e.g., Andreas 1987; Su et al. 2001; Sun 1999; Yang et al. 2008; Chen et al. 2010; Guo et al. 2011). Although a variety of land surface types have been covered in previous endeavors (from smooth snow/ice to rough urban surfaces), a complete evaluation that encompasses the three traditional iterative schemes (including Businger, Dyer, and BHH) and the recently developed noniterative scheme (see Li et al. 2010), however, is still not available. Although the new scheme by Li et al. (2010) has been incorporated in efforts to propose new flux parameterizations (see Wouters et al. 2012), its performance has not been verified against field data. Our work thus provides the first data for evaluating the new scheme, while being equally concerned with traditional schemes that are still being used in numerical models.

## 4. Conclusions

Based on micrometeorological measurements over a tall grass surface in a Beijing suburb, we investigated the radiation budget and the available energy partitioning near the surface. The aerodynamic roughness length (*z*_{0}) over this grass surface was determined to be 0.02 m and the zero-plane displacement height (*d*) was determined to be 0.44 m. Both *z*_{0} and *d* values were then used as input parameters of an analytic model of flux footprint; we confirmed that eddy-covariance flux measurements are mainly representative of the tall grass surface herein, potential influences from anthropogenic sources notwithstanding. The human-induced influences, nonetheless, were regarded as relatively insignificant, because the measured CO_{2} fluxes during daytime hours were directed downward herein (the surface was a sink of CO_{2}), leading to a phenomenological contrast with other (sub)urban areas that are commonly found to be a source of CO_{2}.

The selected fair-weather data give a relatively satisfactory value of energy balance ratio (0.89), and these data are thus used to evaluate four parameterizations of turbulent surface fluxes, namely, a total of three traditional iterative schemes and one noniterative scheme developed to reduce computational time. Their performance is intercompared in terms of the estimations of the sensible heat flux and two turbulence components (the friction velocity and temperature scale). The following conclusions are drawn.

All the candidate schemes lead to rather inaccurate estimations in strongly stable conditions (with high bulk Richardson number higher than 0.2), primarily owing to the heat fluxes of relatively low magnitude.

In weakly stable to unstable conditions, two schemes are recommended herein because of their accuracy; the first scheme stems jointly from the work of Högström (1996) and Beljaars and Holtslag (1991), and the second stems from that of Li et al. (2010).

A sensitivity analysis concerns the ratio of

*z*_{0}/*z*_{0h}(*z*_{0h}is the thermal roughness length); for this tall grass surface, the choice of*z*_{0}/*z*_{0h}= 100 is more appropriate than another choice of*z*_{0}/*z*_{0h}= 10 because the former choice leads to more realistic relationships between the bulk transfer coefficients (*C*and_{M}*C*) and the bulk Richardson number (Ri_{H}_{B}).

To conclude, given the fact that the recently developed scheme by Li et al. (2010) has the advantage of being explicit (no numerical iterations are needed), we therefore recommend this scheme for future use in related atmospheric simulations.

## Acknowledgments

This study was supported by the National Key Basic Research Program under Grants 2010CB428502 and 2012CB417203, by the China Meteorological Administration under Grant GYHY201006024, by the CAS Strategic Priority Research Program Grant XDA05110101, and by the National Natural Science Foundation of China under Grant 40975009. Xiaofeng Guo acknowledges the support of State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences (Grant LAPC-KF-2009-02). Elie Bou-Zeid would like to an anonymous donor for financial support of research efforts related to water challenges in China. We are very grateful to three anonymous reviewers for their careful review and valuable comments, which led to substantial improvement of the manuscript.

## REFERENCES

_{2}fluxes over grassland in the tropical monsoon region of southern China. J. Geophys. Res., 112, D10106, doi:10.1029/2006JD007889

_{2}fluxes over typical steppe prairie in Inner Mongolia, China

_{2}fluxes over a suburban lawn: The influence of traffic emissions

_{2}flux measurements

_{2}flux over an urban area in Beijing

_{2}emissions from cities