A method based on high‐frequency temperature measurements to estimate the sensible heat flux avoiding the height dependence

A new method to estimate the sensible heat flux H using high‐frequency temperature data has been proposed. The new method proposes to scale the mean eddy vertical velocity responsible for the air parcels renewal using similarity formulae. It is shown that the empirical coefficient resulting from this scaling (apparently height dependant) is relatively constant with height when comparing estimated versus measured H over half hour periods. This work suggests that the coefficient is universal. This allows to us to propose a new method that provides useful advantages for field applications.


Introduction
[2] Direct measurement of sensible heat flux density is possible with fast response turbulence instruments, which are able to sense the fluctuations of wind speed and temperature. This method, known as eddy covariance method, is limited by its stringent instrumental requirements and by the high cost in terms of money and maintenance. Hence many indirect methods have been developed to estimate the sensible heat flux from accurate measurement of several meteorological variables using more robust and economical instrumentation. The fundamentals and relationships of some indirect methods (e.g., Bowen ratio method, aerodynamic method, etc.) can be found in traditional textbooks like Brutsaert [1982], Stull [1988], and Monteith and Unsworth [1990].
[3] Paw U and Brunet [1991] pioneered the surface renewal method [Higbie, 1935] to estimate the sensible heat flux over natural surfaces. To apply this method, temperature data measured at high frequency (typically 2 -10 Hz) are needed. The method is attractive because it uses only a low-cost fine wire thermocouple, and it is extremely interesting because it bypasses the difficulties of the traditional methods. Methods based on similarity theory are difficult to apply in sparse, tall canopies, and they cannot be used within the canopy. The eddy covariance method cannot be used within very dense, short canopies like grass, mulch, etc. However, despite the advantages in terms of maintenance and low cost as mentioned above that permit and excellent spatial cover, to increase the number of experiments, etc., for long-term monitoring, this approach may not be convenient owing to the huge amount of data recorded.
[4] Organized structures have been described for both above and within the plant canopy [Gao et al., 1989;Paw U et al., 1992]. These structures are responsible for the majority of the vertical momentum, heat, and mass transfer . This fact provides the key to estimate fluxes and determines the basis of the surface renewal model to estimate the sensible heat flux over natural surfaces. The model assumes that the turbulent exchange of any scalar is driven by the regular replacement of the air parcel that is in contact with a surface, where the exchange occurs. An air parcel from above the surface sweeps to the surface and replaces a parcel that ejects from the canopy Katul et al., 1996;Snyder et al., 1996;Spano et al., 1997;Chen et al., 1997a].
[5] When high-frequency temperature data are plotted versus time, the organized coherent structures look like ramp events. A comprehensive, idealized scheme of the process (referred as scheme 1) is illustrated in the papers of Paw U et al. [1995], Snyder et al. [1996], Katul et al. [1996], and Spano et al. [1997]. Under unstable conditions (scheme 1 in Figure 1) the ramp is characterized by a nearly constant increase of temperature as the air parcel is heated by contact with the surface and a sharp drop associated with a nearly instantaneous replacement of a cooler air parcel from above. This is followed by a quiescent period when there is no heating or cooling. For stable conditions a warm air parcel sweeps in from above and replaces the parcel in contact with the surface (Figure 1).
[6] Chen et al. [1997a] proposed another scheme (scheme 2 in Figure 1), which accounts for a finite microfront instead of a near instantaneous drop in the temperature. They assumed that the quiescent time period is negligible (Figure 1). For a given period of time the sensible heat flux H at any height z is the net exchange of heat conducted by all the ramps formed during this period. From Paw U et al. [1995] the sensible heat flux at any z is obtained by where r and C p are the density and specific heat of air, respectively, and A and L are the mean ramp amplitude and duration between ramps over the sampling period. As shown in Figure 1, L is the sum of a quiescent time period and the ramp period when the parcel is being heated or cooled. Equation (1) represents the sensible heat flux density lost through the top of the air parcel.
[7] Details of the derivation of equation (1) are given by Paw U et al. [1995]. The measurement height z is equal to the mean volume per unit area. The variable a corrects this mean volume for unequal heating. If the air parcels are uniformly heated at all heights within the air volume under the measurement height, then a = 1.0. Following Snyder et al. [1996], this is realistic when the measurements are made well above the canopy top into the inertial sublayer. However, when the measurements are made close to the canopy top, the a parameter takes into account the vertical lapse rate of temperature in the air parcel. When the measurements are taken just above the canopy top h, Paw U et al. [1995] suggested 8z = 0.5h for tall, dense canopies (taller than 2 m). Using 8 = 0.5 over a uniform forest terrain and temperature measurements above the canopy top (13 m), Katul et al. [1996] obtained good results. Snyder et al. [1996] and Spano et al. [1997] found a % 1.0 when the data were collected at several times the canopy height over short, dense canopies. Chen et al. [1997b] found that a was close to 0.5 (0.45-0.7) either at or above the canopy height. In addition, Chen et al. [1997b] suggested that the variability in the value of a around 0.5 reported by Snyder et al. [1996] and Spano et al. [1997] is mainly due to the method used to determine the residence time of the parcel in contact with the heat sources (or sinks) since it is overestimated by a factor $2.5.
[8] To determine A and L in equation (1), Snyder et al. [1996] and Spano et al. [1997] used a method proposed by Van Atta [1977] that will be referred as VA. Then, small errors in the determination of the amplitude of the ramp result in large errors when determining the duration of the ramp because it depends on the cube of the amplitude (see equation (A5)). Other authors have used different methods. However, as mentioned above, Chen et al. [1997a] used scheme 2 rather than scheme 1 to identify ramps and used a higher measurement frequency (80 Hz) than the other studies. Consequently, the two methods are not directly comparable.
[9] When the transfer of momentum to the canopy is high, near the canopy top, ramp events may not be detectable if the sampling frequency is too low. This is a crucial difference between the papers presented by Chen et al. [1997aChen et al. [ , 1997b, Snyder et al. [1996], and Spano et al. [1997]. They found ramps when measuring close to the surface, but their sampling frequency was $10 times as fast as the others. However, the surface renewal method becomes impractical if we need to sample more frequently than a reasonable capability. Hence, in some situations, it is desirable to set up the instrumentation well above the canopy top. From Snyder et al. [1996] and Spano et al. [1997] it can be inferred that the assumption a z = z does not always lead to accurate estimations. Snyder et al. [1996] and Spano et al. [1997] found good results after calibration of the a parameter, but a varied with the height above the crop and into the inertial sublayer. This suggests that there was unequal heating even well above the canopy. No rule was given to estimate the a parameter for a given measurement height. Calibration using a sonic anemometer was recommended.
[10] The objective of this work is to propose a rather empirical model to estimate the sensible heat flux as an alternative to the surface renewal model (equation (1), scheme 1 using 8 Hz frequency) that bypasses the problem associated with the measurement height. Different methods have been reported in the literature to determine the ramp parameters; however, the method of VA is used in this paper.

Theory
[11] Over a homogeneous surface, turbulent transfer of heat is predominantly vertical, and the turbulent diffusion approach, commonly known as the K theory, has been successful in describing the sensible heat flux.
where T is the time averaged temperature of the air, K h is the turbulent transfer coefficient or eddy diffusivity for heat, and z is the height from the zero plane displacement. The justification of equation (2) rests on empirical grounds and the analogy between the mean free path of molecular motion and the scale of turbulence. The idea of relating fluxes to gradients is also implicit in similarity theory for a uniform boundary layer. In turn, the variable (z a) in equation (1) physically represents the ''effective eddy size'' responsible for the air parcel renewal. Then, since the ramplike structure (characterized by the amplitude A over period L) contributes to vertical transport, the local vertical temperature gradient must scale with A/(z a). It can be derived that This is based on the surface renewal concept that the air descending to the plant canopies comes from some typical height with some scalar value. Then, when the parcel of air descends to the canopy and begins to be enriched or depleted of the scalar, the time rate of change of the scalar exhibits the slow change part of the ramp. The magnitude of the sudden drop or rise, at the peak of valley of the slow ramp change, defines amplitude A. If one takes the baseline scalar value as representing that of the scalar at the typical height from where the air parcel originated, then the ramp amplitude A should be directly proportional to the mean temperature difference between the parcel origination height and the canopy height. Therefore the mean air temperature gradient should be directly proportional to A, as expressed in equation (3). Pereira and Paw U [1994] briefly discussed this concept as applied to temperature differences between the effective surface and the air, which represents a bulk temperature gradient.
[12] When measurement height is well above the canopy top, in practice z > h + 2(h À d) [Sellers et al., 1986], where d is the zero plane displacement, the Monin-Obukov similarity theory can be used. Therefore a suitable expression for K h is the following: This parameterization involves large-scale eddies, and for practical calculations, it is usually assumed constant over sampling periods about half an hour or an hour.
[13] From equations (2), (3), and (4) a new relationship to determine the sensible heat flux into the inertial sublayer can be expressed as We have assumed here that the height at which parcels are brought into the canopy, by coherent structures, is at least as high as the height at which similarity theory may be assumed. In principle, b 1 is a height-dependent unknown parameter, as is a. In equation (4), z is the elevation above the zero plane displacement, k is Von Karman's constant (taken as 0.4), u * is the friction velocity, and f h (z) is the stability function of the surface layer on the heat transfer processes. The widely accepted formulation of f h (z) is [Businger et al., 1971] Unstable where z is a dimensionless buoyancy parameter defined as z = z r /L MO , z r is a reference height, and L MO is the Monin-Obukov length that can be expressed as where T is the air mean absolute temperature and g is the gravitational acceleration. Friction velocity was calculated using the expression [Dyer, 1974] where u r is the wind speed at the reference height z r , z o is the surface roughness height, which can be estimated as 0.1h [Wieringa, 1993], and C m is the diabatic profile function for momentum where x = (1 À 16z) 1/4 .
[14] When measurement height is in the roughness sublayer, where thickness is assumed to be between z = h and z < h + 2(h À d) following Sellers et al. [1986], the turbulence is mainly dominated by mechanical factors or by canopy roughness properties. It is not the focus of this work to study the best parameterizations of K h since it still is a subject of intense research [Lee, 1996]. A linear relationship between K h and ðu * zÞ is proposed. Then the new formula to estimate the sensible heat flux in the roughness sublayer can be expressed as where, for the same reasons stated for equation (5), in principle, b 2 is a height-dependent unknown parameter. [15] To estimate the sensible heat flux using equations (5) and (10), determination of the ramp amplitude is needed. Chen et al. [1997a] inferred that VA is a suitable method to determine A for short canopies. In addition, possible inaccuracy in determining the ramp period is eliminated. Therefore the VA method was used to determine A. The A value used was obtained as the mean from three A values using different time lags as suggested by Chen et al. [1997a]. The time lags selected were 0.25, 0.5, and 0.75 s. They were inferred from their results. For bare soil and straw mulch (6.5 cm thick) the optimum time lags were in a narrow range around 0.1 s, and they were in a wide range around 1 s for a Douglas fir forest (16.5 m tall).
[16] To determine u * and z, iteration of equations (7), (8), and (9) was done setting z = 0 to start the process. Initially, u Ã was calculated from equation (8), then z was calculated after calculation of L MO from equation (7), where the measured sensible heat flux with a sonic anemometer was used. The procedure is repeated until convergence is achieved. Three iterations were, in general, needed for unstable conditions.

Data
[17] Three different crops were selected to test the method proposed. Two of them, grass and wheat, were with a homogeneous canopy. The other one was grape vineyard. The grass experiments were conducted over 0.1 m tall grass at Davis (California). The measurement heights were at 0.6, 0.9, and 1.2 m above the ground during the springtime period (from day of year 86 to 88) and at 0.7, 1.0, and 1.3 m above the ground during the summertime period (from day of year 213 to 214). Data were collected on day 86 from 1400 to 1730 LT, on day 87 from 1100 to 1700 LT, and on day 88 from 1000 1530 LT. Data were collected on day 213 from 0900 to 1500 LT and on day 214 from 0900 to 1300 LT. The wheat experiments were conducted over 0.7 m tall wheat also at Davis. The measurement heights were at 0.7, 1, and 1.3 m. Data were collected during two days in spring season (day 148 from 1000 to 1800 LT and day 149 from 0930 to 1930 LT). Details about the location and the experiment are given by Snyder et al. [1996] for grass and by Spano et al. [1997] for wheat.
[18] The grapevine experiment was conducted at the Oakville Field Station in Napa Valley, California (38°26 0 N, 122°24 0 W, 58 m above sea level). The Cabernet Sauvignon grapevines were oriented in north-south rows and were $2.0 m tall, and there was $100 m of fetch in the upwind direction (south) during the experiment. The vines were irrigated with a drip irrigation system and trained in a conventional curtain system. Ground cover was $60%. Temperature data were collected at the canopy top (2.0 m) and at 2.3, 2.6, and 2.9 m.
[19] In all three experiments, temperature data were measured using a 7.6 Â 10 À5 m diameter fine wire thermocouple and recorded at 8 Hz. In addition, standard meteorological data were measured nearby using automated weather stations. Half hourly data sets were created in order to compare estimated versus measured sensible heat flux using a one-dimensional sonic anemometer.

Results
[20] The performance of the different estimates using the proposed method, equation (5) or equation (10), and the original method, equation (1), is shown in Table 1 for the different measurement heights. Simple linear fitting analysis of the estimates versus the actual values was used. A linear fit was done to determine the intercept to see more clearly any bias of the methods. Another fit was forced through the origin to find a ''practical field value'' for each method, level, and crop. The root-mean-square error (RMSE) was determined using the corresponding practical field value to determine the reliability of each method. From the practical field values determined at each level a value was assigned as ''a rule of the thumb'' for each crop and will be referred as the crop value.

Grass
[21] The top of the roughness sublayer is estimated to be at 0.16 cm, so clearly all the data gathered were measured in the inertial sublayer. Then equation (5) was used to determine the sensible heat flux at each measurement height. The results obtained from the whole data set are shown in Figure 2. Table 1 shows the statistics of the linear fit of the estimates using both methods in two data sets, spring and summer.
[22] Overall, the performance of equation (5) is excellent. Whatever the level of measurement, the slope is nearly constant, and the bias can be neglected; it is statistically not different from zero at 5% level of significance. The regression coefficients are, in general, high. Similar results were obtained for the linear fits through the origin. A crop value for this crop was well defined since it was not height dependent (see Table 1). The results obtained from equation (1) were not good in general since the statistics determined depend on the day and height. The best results for the spring season data set are obtained at 0.9 m level, while for summer season data set they are obtained at 0.7 m level. It is shown that the practical field value or a parameter is far from constant with height.
[23] The proposed method was superior to equation (1) since at each level and, consequently, for all the data sets the root mean square error was lower. The crop value using the equation (5) was 0.1 and 1.0 using equation (1), as suggested by Snyder et al. [1996] for measurements well above the canopy top.

Wheat
[24] The top of the roughness sublayer is estimated to be at 1 m, so the lowest level was located in the roughness sublayer, while the upper level was located at the bottom of the inertial sublayer. The measurements were made close to the transition zone, and a question arises concerning which equation should be used to test the proposed method. In order to see the different performance of equations (5) and (10) when they are not applied in their corresponding sublayer, both equations were used to determine the sensible heat flux at each measurement height. The results obtained are shown in Figures 3 and 4 for equations (5) and (10), respectively.
[25] Despite some measurements made in the roughness sublayer, equation (5) gave, in general, the best performance (see Table 1). At all levels, equation (5) was superior to equation (10), and the variation of b 1 was lower than b 2 . Noteworthy is the low variation of coefficient b 1 with a value similar that the obtained for the grass data set. Despite equation (1) giving the best performance of the sensible heat flux at the 1.0 m level, it gave the worst performance for the others. The variation with height of the a parameter was  high for this canopy with a variation coefficient of 19%. For all data sets, the crop values assigned for this canopy could be the same as that assigned for grass for the equations (5) and (1).

Grape vineyard
[26] It is difficult to estimate the top of the roughness sublayer in this case because the canopy is not homogeneous. In this case the sensible heat flux comes from the soil surface and from the canopy. Assuming the bottom of the inertial sublayer at a height twice the height of the canopy [Raupach et al., 1989], all measurement levels were located in the roughness sublayer. As was done in section 3.2.2 for wheat, in order to see the difference in performance of equations (5) and (10) when they are not used in their corresponding sublayer, both equations were used to determine the sensible heat flux at each measurement height. The results obtained are shown in Figures 5 and 6 for equations (5) and (10), respectively.
[27] Similar to what happened with the wheat data set and despite the measurements being made in the roughness sublayer, equation (5) gave, in general, the best performance (see Table 1). Whatever the level, the standard error of the estimate was greater, but the bias was low enough to introduce the lowest root mean square error rather than using equations (10) and (1). The variation of the coefficients b 1 and b 2 can be assumed constant with height. Noteworthy is that the b 2 values are close to that obtained for wheat at the level of 0.7 m. Recall for wheat that the lowest level was in the roughness sublayer and some hesitation may arise from the 1 m level. The a parameter of equation (1) was nearly constant in the three lowest levels (2, 2.3, and 2.6 m), although not close to 0.5 or to unity, $0.85. For the highest height (2.9 m) it was notably lower, 0.74. These results disagree with the assumption that the higher the measurement level, the higher the a parameter, with mean boundaries values of 0.5 and 1 for close and well  above the canopy top, respectively. For all data sets, the crop values for this crop were different to that assigned for the other canopies, but for the proposed method the variation was smaller than using the original method. Noteworthy is the good performance of the two methods for this canopy since it assumes the energy advection term negligible.

Summary and Concluding Remarks
[28] The proposed new method to estimate the sensible heat flux empirically links the standard surface renewal method (Lagrangian in nature) proposed by Paw U et al. [1995] with a quasi-stationary one-dimensional diffusion process with a constant heat eddy diffusivity over a half hour time period.
[29] Here, using the Van Atta [1977] analysis method, it is shown that a suitable vertical velocity scale is proportional to u * f(z). Despite the fact that the proposed scale is well defined in the inertial sublayer, good performance was also obtained in the roughness sublayer, which is difficult to explain. This is similar to the results obtained by Chen et al. [1997b], who proposed to scale the vertical velocity (az/L) with the friction velocity. It is not the objective of this work to find the best parameterization of the eddy diffusivity. However, the stability function takes into account the increase or decrease of the vertical Figure 6. Comparison between the sensible heat fluxes H 2 (W m À2 ) estimated using equation (10) for grape vineyard and the measured sensible heat flux H m (W m À2 ). The 1:1 line is represented for comparison. Figure 5. Comparison between the sensible heat fluxes H 1 (W m À2 ) estimated using equation (5) for grape vineyard and the measured sensible heat flux H m (W m À2 ). The 1:1 line is represented for comparison.
velocity under unstable or stable conditions respectively. In other words, the air parcels are renewed by eddies whose length scale is comparable to z for near-neutral conditions but larger or lower than z for unstable or stable conditions, respectively. This supports some results obtained by Katul et al. [1996].
[30] However, the proposed method is not as attractive as the original from Paw U et al. [1995] since it is air velocity dependent. Nevertheless, the proposed method achieves the aim of this work with the scope on practical field practices since a cup anemometer is not expensive. For practical field applications, after calibration of the crop value, the iterative procedure to determine the stability function can be done using the sign of the ramp amplitude that determines the sign of z (Figure 1).
[31] This work suggests that the coefficient b 1 in equation (5) is universal for homogeneous short canopies (with a value of 0.1); thus the temptation arises to suggest that it may be stable for a given homogeneous canopy. In addition, it is noteworthy that the crop value obtained for grapevines was 0.15 despite the large differences of the two canopies and turbulent nature properties of the sublayer where the measurements were carried out. Also, the general good performance of the equation (10) in the roughness sublayer should be pointed out since coefficient b 2 was also stable, ranging from 0.23 to 0.33, for the different canopies.
[32] To conclude, a new empirical method based on highfrequency temperature measurements has been proposed which after calibration and despite being based on somewhat rough physical assumptions, is able to estimate the sensible heat flux with reasonable errors (see Table 1). For all data sets, the errors were lower than the original surface renewal method proposed by Paw U et al. [1995]. The main advantage of the new method is based on the fact that vertical velocity of the mean eddies responsible of the renewal process has been properly scaled for the corresponding amplitude of the temperature of the mean ramp events. Although it cannot be confirmed with the data sets used in these experiments, the results suggest that temperature sensors can be mounted at any height above the canopy. In the roughness sublayer, coefficient b 2 is nearly constant as well, taking into account that the canopies of wheat and grapevines are very different. This supports the results obtained by Chen et al. [1997b], which used a different ramp-like structure scheme, frequency of measurement temperature data, and a variety of canopies very different from the used in this work. The horizontal wind speed is needed, so the proposed method is less attractive than the original one. However, it is noteworthy that the ramp period is not needed, avoiding a source of error in the estimates of the sensible heat flux.

Appendix A
[33] Structure functions, equation (A1), and the analysis technique, equations (A2)-(A5), from Van Atta [1977] were used to determine A and L: where m is the number of data points in the 30 min interval measured at frequency f (in Hz), n is the power of the function, j is a sample lag between data points corresponding to a time lag (r = j/f ), and T i is the ith temperature sample. According to Van Atta [1977], the time lag r must be much less than l + s. An estimate of the mean value for amplitude A during the time interval is determined by solving equation (A2) for the real roots where p ¼ 10 S 2 ðrÞ À S 5 ðrÞ S 3 ðrÞ ðA3Þ q ¼ 10 S 3 ðrÞ: The methodology to solve equation (A2) for real roots of A is given in most mathematical handbooks. The inverse ramp frequency L was calculated using [34] Acknowledgments. This work was supported by the Centre Tecnològic Forestal de Catalunya under project InterReg II (2.3.GC. code 82), by the CICYT projects HID96-1295-C04-03 and HID97-397, and a grant from the CIRIT. The authors want to gratefully acknowledge K. T. Paw U for his valuable comments and R. L. Snyder and D. Spano for providing the data used in this work. In addition, the suggestions made by the anonymous reviewers improved substantially the readability of the manuscript.