THERMAL STORAGE RATES IN THE SURFACE ENERGY BALANCE

L.W. Gay, School of Renewable Natural Resources, UAZ, Tucson, AZ 85721. Telephone 502/621-7268; email: lgay@ag.arizona.edu. 7/9/98.

ABSTRACT

The thermal storage rate (G, W/m2) of a vegetated surface can be estimated from 3 components: G = Ga + Gveg + Gs [Eq. 1], where Ga is air layer storage rate, Gveg is biomass storage rate, and Gs is the soil storage rate. We ignore energy storage associated with changes in vapor content of the air layer. Fluxes are positive coming out of storage (out of air layer, out of biomass, out of soil), as commonly occurs at night. G can be estimated by measuring Gs directly with soil heat flux plates, and by estimating Ga and Gveg with respect to a single air temperature (T, C) above the canopy. Ga and Gveg are linear functions of the time change of the index temperature T. Then Gx = Kx(Tp - Ts) [Eq. 2],where Kx is a thermal storage coefficient, subscript "x" represents "a" (air) or "veg" (biomass), Tp is air temperature for the "previous" averaging period and Ts for the "subsequent" averaging period.

The first step is to evaluate Ka and Kveg. The storage coefficient Kx = (r xCpx)(hZ/D t) [Eq. 3], where r is density (kg/m3) and Cp is heat capacity (J/kg K), hZ (m) is thickness of the air layer or vegetation layer, and D t (s) is time interval between midpoints of previous and subsequent data averaging periods. We consider first Ka since air density varies with temperature and with pressure (i.e., elevation). Cp for air is essentially constant at 1010 J/kg at ambient temperatures and pressures. Density of dry air (kg/m3) can be estimated as a function of pressure (P, in kPa) and temperature (T, K) from the equation (List, 1949, p. 290) r (kg/m3) = 3.4838*P (kPa)/T (K) [Eq. 4a]. Temperature can be taken as the mean daily air temperature. Pressure PZ at elevation Z (m) of the measurement site can be estimated from the International Standard Atmosphere (List, 1949, p. 268) and sea level air pressure (Po = 101.33 kPa), as PZ/Po = ((288-0.0065Z)/288)5.256 [Eq. 4b]. Air density determined by Eqs. 4 A&B is used in Eq. 3 to determine the air storage coefficient for that site. Calculations will be demonstrated for the HARTX Scots pine forest site in southwestern Germany. On DOY 135 mean Ta = 20 C = 293.15 K, Z = 200 m, hZ = 17 m, and D t = 1 hr = 3600 s. From Eqs. 4 A&B we calculate r a = 1.176 kg/m3. Finally, Eq. 3 yields Ka = 5.61 W/m2 K. Ka follows mean air temperature, and can change daily.

We now calculate the storage coefficient in vegetation from Eq. 3. Biomass (typically in kg/m2) must be expressed as density (kg/m3) over the thickness of the vegetation layer. The density is then adjusted for an "active fraction" of about 0.6 for forests, and Cp is assumed equal to that of water (4190 J/kg K). If biomass at HARTX is 10.1 kg/m2 and hZ = 12 m, the average density for the vegetated layer is 0.842 kg/m3 and the active density becomes 0.505 kg/m3. D t remains unchanged. Thus Kveg = 7.053 W/m2 K from Eq. 3, and this remains constant until the density changes.

Ka and Kveg can now be applied to the observed previous and subsequent index temperatures in Eq. 2 to calculate present values of Ga and Gveg. The values are combined with observed Gs in Eq. 1 to yield G (W/m2).

INTRODUCTION

Only four thermal energy fluxes are needed to define a surface energy balance with reference to an infinitely thin, horizontal plane positioned near the earth's surface. The four fluxes are net exchange of allwave radiation (Qn), thermal storage flux (G), sensible heat flux (H) and latent heat flux (LE). These fluxes are expressed in flux density units of W/m2. Since the thin plane can store no energy, flux polarity in an energy balance analysis is defined as positive for transfer to the plane, and negative for transfer away. The definition is straight forward for radiation, sensible heat and latent heat, but application to the thermal storage fluxes warrants further explanation. The flux of thermal energy in and out of the soil (or air below the plane) is the rate of change of thermal storage in the soil (or air layer). During daylight this energy flow moves into the soil (or air layer) and away from the reference plane so it has negative polarity. At night, the flow of thermal energy is from the soil (or air layer) up to the plane and it then is positive. This defines an increase in thermal energy in the soil (or air layer) as negative with respect to the plane, while a decrease in the layer (i.e., energy coming out of storage) is positive.

The polarity is important because the fluxes obey the law of conservation of energy in that gains of energy must equal losses at the surface for any time period. Consequently, the surface energy balance is written (with appropriate polarity) as follows:

Qn + G + H + LE = 0.

This note describes evaluation of the thermal storage flux. This is normally the smaller of the four fluxes, and can approach zero when averaged over a daily cycle of 24 hours, and over an annual cycle of 365 days, and so the mean daily storage flux is often assumed to be zero. However, for short periods in daytime (say periods of an hour or so) and for certain surface conditions (bare soil, water, high forest, etc.) storage fluxes can become a significant component of the surface energy balance.

Precise measurements of "available energy" (A = Qn + G) are important for energy balance measurements of LE. The Bowen ratio energy balance model is such an energy balance method; it divides up the available energy by a ratio of gradients of air temperature and humidity (Tanner, 1960). In a different energy balance application, Qn and G are measured directly, sensible heat H is measured by eddy correlation or by an aerodynamic model, and LE is estimated as the residual term in the energy balance equation. The estimates of thermal storage fluxes illustrated in this report were prepared for use with such an "eddy correlation/energy balance" measurement system, labeled the One Propeller Eddy Correlation (OPEC) system by Blanford & Gay (1992)). With such systems, accurate measurements of the storage fluxes contribute directly to accuracy of latent energy fluxes.

THE THERMAL STORAGE FLUX

This note defines and demonstrates practical methods to estimate thermal storage fluxes for use in an energy balance analysis of a forest site. The primary reference is to the development of a 5-month energy budget of a Scots pine forest (Gay, et al., 1996) in the Hartheim Experiment (HARTX) where thermal storage fluxes were evaluated separately by two independent sets of measurements.

Storage Flux Components. The rate of change in thermal storage (G, W/m2) can be separated into three components:

G (W/m2) = Ga + Gveg + Gs (1)

where Ga is the rate of energy storage within the layer of air between the soil surface and the sensors above the canopy, Gveg is the rate of energy storage within the mass of vegetation, and Gs is the rate of rate of energy storage in the soil. The fluxes become larger at forest sites where both the biomass and the thickness of the air layer increase. The three storage flux components in Eq. 1 (air, vegetation and soil) are functions of temperature and specific heat capacity of the three media (air, vegetation and soil). Equation 1 excludes change in thermal energy associated with changes in vapor content of the air layer for simplification, since the OPEC system measures air temperature, but not water vapor content. Effects of this simplification are small.

A Practical Storage-Flux Model. This simplified model estimates the magnitude of the individual storage fluxes in Eq. 1 as linear functions of environmental temperature. The data used in this analysis are means for a data averaging period that typically ranges from 0.2 to 1.0 hr. The form of the model is given in Eq. 2:

Gx = Kx(Tp - Ts), (2)

where Kx [W/m2K] is a thermal storage coefficient, subscript "x" is either "a", "veg", or "s", T [K] is temperature, and subscripts "p" and "s" identify temperatures in the previous and the subsequent data-averaging periods. The use of previous and subsequent data-averaging periods in Eq. 2 defines the "central difference" method described by Stewart (1988).

The storage coefficient is essentially "volumetric heat capacity" for the layer thickness of the storage medium as averaged over the time increment of the analysis. The thermal storage coefficient for air, vegetation or soil is defined in Eq. 3 for sea level pressure (Po = 101.33 kPa):

Kx [W/m2K] = (r xCpx)(hz/D t) (3)

where subscript "x" is either "a", "veg", or "s"; the volumetric heat capacity [J/m3 K] is the product of the appropriate density r [kg/m3] and heat capacity Cp [J/kg K]; hz [m] is the thickness of the storage layer; and D t [s] is the time interval between midpoints of the previous and subsequent data-averaging periods. It is also feasible to use "backward differences" (previous to present period) and "forward difference" (present to subsequent period). These two alternatives have a time period half as large as the central difference method.

Note that density in the air storage coefficient varies with temperature and with air pressure. These adjustments will be explained in the discussions that follow. Also, all three storage coefficients Kx incorporate values of hz and D t and so the coefficients are site specific.

Measuring Temperature Changes. Temperature changes (Tp - Ts, units K) in storage media of the general model (Eq. 2) define the magnitudes of the air, vegetation and soil storage flux components. The average layer temperature is desired for Ga and Gs, and this is usually obtained from temperature profiles measured with a number of sensors. With Gveg, however, the vegetation temperatures are not measured directly, but are indexed by temperature of the air layer. Vogt, et al. (1996) discusses application of profile measurements of temperature and moisture in air and soil to evaluate thermal storage fluxes during HARTX.

In contrast to the profile measurements of Vogt, et al., this note describes a simple approximation that used a single air temperature sensor at 17 m at HARTX to index thermal storage changes in the air (Ga) and in the vegetation (Gveg). This approximation method also measured Gs directly with soil heat flux discs rather than by inferring changes in volumetric heat content from changes in soil termperature. Also, storage changes associated with fluctuations in water vapor in the air layer are ignored, since only air temperature data are available from the OPEC system. We shall see that results from the simplified index method compare well with those from the more elaborate system of profile measurements.

APPROXIMATING THERMAL STORAGE COMPONENTS

Storage flux measurements described here seek to accurately define available energy (Qn + G) in the surface energy balance equation in order to improve eddy correlation/energy balance estimates of latent energy. Storage flux calculations with the central difference method are normally done retrospectively after observations are complete, since this method uses a subsequent data point in each calculation. There is seldom a compelling reason to obtain final energy balance results in real time.

Thermal Storage in the Air. In principle, the value of the air storage coefficient Ka can be obtained from Eq. 3. The specific heat of dry air is essentially a constant at Cp = 1010 J/kg K with ambient temperature and pressure. The variables hZ and D t are fixed by conditions of measurement. However, air density (and hence volumetric heat capacity) varies with temperature and pressure and so correct values of air density must be determined for conditions at the measurement site. Density of dry air can be estimated as a function of pressure (P, in kPa) and temperature (T, K) from an equation of List (1949, p. 290):

r a [kg/m3] = 3.4838*P (kPa)/T [K]. (4a)

Note that T (K) = T (C) + 273.15, and kelvin and celsius units are equal for temperature differences D T. Temperature in Eq. 4a can be taken as mean daily temperature observations. Density from Eq. 4a is evaluated at a pressure representative of the elevation Z of the site. Pressure at the site (PZ)can be approximated from the International Standard Atmosphere (List, 1949, p. 268):

PZ/Po = ((288 - 0.0065Z)/288)5.256 (4b)

where Z [m] is the elevation of the measurement site. For example, the ratio PZ/Po = 0.976 at 200 m, and 0.823 at 1600 m elevation. Since sea level pressure is 101.33 kPa, Eq. 4b estimates P200 = 98.898 kPa and P1600 = 83.395 kPa. The density of dry air predicted by Eq. 4a at the ISA pressure for 200 m and at a mean temperature of 293.15 K (20 C) is 1.175 kg/m3.

There are several possibilities for determining the air temperature to be used in evaluating Ka in Eq. 3. First, one can use a spatially averaged temperature to represent the air layer temperature instead of the single air temperature used as an index in this paper. This is normally obtained by calculating a mean for temperature profiles defined through the air layer. Either the daily mean profile or the daily mean index temperature will require only one evaluation of Ka each day, but it may prove awkward to obtain the needed daily mean. .

A satisfactory alternative is to average the previous and the subsequent temperature at each time step to obtain the mean temperature needed to evaluate Ka and Ga at that time step. A new Ka would be determined and used at each time step, but no disadvantage is associated with this. The HARTX air storage calculations using one mean temperature each day gave the same results as those from use of Tp and Ts.

Calculation of an air thermal storage coefficient is illustrated for the HARTX site on May 14 (DOY 135), 1992, using mean daily temperature of 20 C. Dry air density at sea level from Eq. 4a is r a = 1.209 kg/m3 at mean air temperature of 20 ° C; Cp = 1010 J/kgK; hz = 17 m; and P/Po = 0.9765 at the 200 m elevation of the site. The HARTX data-averaging periods are 0.5 hr, so the time difference between previous and subsequent periods (means at the midpoint of each period) is D t = 1.0 hr = 3600 s. Evaluating Eq. 3 yields Ka = 5.63 W/m2K for the HARTX elevation of 200 m. At a higher site elevation, say Z = 1600 m on the Rio Grande of central New Mexico, PZ/Po = 0.823 and Ka = 4.75 W/m2K at a mean air temperature of 20 C. The value of Ka is then applied in Eq. 3 to the observed previous and subsequent air temperatures to calculate the mean air storage rates for each time period throughout the day.

Thermal Storage in Vegetation. The vegetation storage coefficient Kveg can be obtained from Eq. 3, although the volumetric heat capacity of the biomass is evaluated differently from that of the air layer. For the Scots pine forest at Hartheim, Jaeger and Kessler (1996) report biomass measurements of 10.1 kg/m2. With an average tree height of 12 m, this is equivalent to vegetation density of 0.84 kg/m3 over a layer thickness hZ = 12 m. Stewart (1988) simplified his analysis of Scots pine storage by estimating the "active fraction" of the forest entering into thermal exchange to be 0.6 (i.e., the larger elements such as tree trunks do not fully participate in the thermal storage). Heat capacity of the green biomass is assumed equal to that of water (Cp = 4190 J/kgK). The time period for the central difference analysis remains at D t = 3600 s. These values in Eq. 3 yields Kveg = 0.6(0.842)4190(12)/3600 = 7.06 W/m2K at HARTX.

The vegetation storage coefficient Kveg is then applied to the previous and subsequent air temperatures using Eq. 2 to estimate the mean storage rate in the vegetation in the present period. The HARTX air temperatures measured at 17 m height become the index to thermal storage in the 12 m layer of vegetation

Thermal Storage in the Soil. In principle, the soil storage coefficient can be determined with Eq. 3 after estimating the volumetric heat capacity of the soil layer involved in change of storage. To apply Eq. 3 in practice it is necessary to measure profiles of temperature and moisture in the soil layer to determine density and heat capacity. It is particularly difficult to account for changes in soil moisture which can have large effects on volumetric heat capacity. Vogt et al. (1996) describes application of this method at the Hartheim forest site. In the alternative chosen for this report, Gay et al. (1996) buried two pairs of soil heat flux discs at a depth of 1 cm, and recorded Gs directly. The results can be affected by calibration and placement of the soil heat flux discs (Fritschen and Gay, 1979). Use of four discs provides some spatial sampling. Changes in heat storage in the upper 1 cm layer of soil was assumed to be negligible, and mean soil heat flux measured by the four discs was accepted as Gs. Use of soil heat flux discs also eliminated the necessity for measuring soil temperatures

CALCULATING THE STORAGE FLUX G

The storage flux G [W/m2] is calculated by summing the 3 components as given in Eq. 1. The polarity of the storage fluxes is defined by the term Tp - Ts in Eq. 2. In accord with energy balance usage, polarity of storage fluxes in this paper are negative when their sense is away from the exchange surface (i.e., into the air, into the vegetation, and into the soil as at midday) and positive when moving to the exchange surface (out of the air, out of the vegetation, and out of the soil, as at night). As a note of caution, some authors assign the polarities of fluxes G, H and LE by direction where up is positive and down is negative. Flux polarities should always be viewed within the context of their intended use.

Results of Sample Calculations. This section demonstrates the evaluation of the thermal storage components at the HARTX forest site. The focus is upon short-term 30-min means within individual days. Two days of thermal storage rates calculated at HARTX (30-min mean data) are presented as examples in Tab. 1. The only information required to recreate these calculations is the sequence of 30-min mean air temperatures in column "T2", and the parameter values that were given in the earlier HARTX examples. The 30-min mean net radiation (Q) is included to confirm that skies were clear on these two days and thus the magnitudes of thermal storage should be near maximal values for the dates. The peak 30-min means of Q exceed 700 W/m2 in spring, and 500 W/m2 in fall. The first day shown in Tab. 1A, May 14 (DOY 135), 1992, occurs during spring warming; the second day in Tab. 1B is Sep. 13 (DOY 257), occurs during the period of near-neutral thermal exchange at the time of the autumnal equinox.

The thermal storage components in Tab. 1 are plotted as time series in Fig. 1 to illustrate magnitudes at HARTX on these two days. The thermal storage fluxes are plotted as "stacked curves", piggy-backed upon each other so Fig. 1 contains curves for Gs, for Gs+Ga and for G = Gs+Ga+Gveg. Looking first at the total thermal storage rate G, the peak 30-min mean gain near the middle of the day reaches about -80 W/m2 in spring, and -50 W/m2 in autumn. When averaged over the entire day, G shows a gain of -14 W/m2 in the spring and only -3 W/m2 in the fall. The values are small, but consistent with the cyclic supply of solar energy which increases rapidly in late spring, and diminishes substantially by the autumnal equinox.

Figure 1 confirms that soil heat flux (Gs) is the largest of the 3 storage components at the HARTX site, with 30-min peaks reaching -60 W/m2 in the spring and -26 W/m2 in the fall. Also, if we integrate the heat storage rate G

This page contains Table 1 A&B, tabulating 30-min means of thermal storage components at HARTX for a spring day and a fall day. The table "TMAYGD.doc" was printed as a text file on DOS-based printer at 16.67 cpi (132 characters on 8.5 inch paper). The file can be opened in W95 and fits on a page with 0.5 inch margins all the way around when using font size 10. The columns are wavy in W95, but there should be a way around this.

This page contains Figs 1A & 1B, showing HARTX thermal storage components throughout a clear day in spring, and again in fall.
over the 24-hour period to determine daily totals, the largest contribution comes from the Gs component. Daily sums in the air layer (Ga) and in the vegetation (Gveg) tend toward zero. This provides a rationale for mean daily studies to concentrate upon Gs, and place a lower priority upon measurements of Ga and Gveg., especially at sites with continuous, uniform canopies as at HARTX.

Comparison of Methods at HARTX. This section compares results obtained at HARTX with two independent measurements of the thermal storage rate. The first method is the simplified approach of direct measurement (Gs) and indexed estimates (Ga and Gveg) from a single air temperature well above the surface. These results are compared to those from Vogt's, et al. (1996) more elaborate analysis of measured profiles of temperature and moisture through the soil, plant and air layers. The HARTX results are compared over two periods. The spring comparison extends from DOY 135.67 through 143.46 (length of 8.8 days). This contained 438 30-min periods, of which concurrent data were obtained on 405 periods. The fall comparison extended from DOY 254.00 through 274.0 (length of 21 days) over 1008 30-min periods. The two systems operated concurrently for 960 30-min periods in the fall.

The first comparison is illustrated in Fig. 2A with a time series plot of the total thermal storage from the simplified method (G, displayed as a single line) with that from the profile method (G'). A similar comparison is shown in Fig. 2B for the fall period. The fall period extended 21 days, but only 10 days are plotted in the time series in order to show both periods at the same time scale. The two estimates of G are surprisingly close together in both spring and fall, given that they are based on somewhat different measurements approaches at two locations separated by 50 m within the Hartheim forest.

Statistical analysis are illustrated in Fig. 3A for comparisons in spring, and in Fig. 3B for fall. The well defined linear relationship is G' = 0.97G +1 in spring (R2 = 0.87), s y.x = ± 11 W/m2), and G' = 0.84G -1 in the fall (R2 = 0.70, s y.x = ± 11 W/m2). A phase shift between G and G' is evident in the fall regression in Fig. 3B. If the data for G' is lagged (retarded) by one hour, the fall regression improves to G' = 0.91G - 1, with s improving to ± 8 W/m2, and R2 increasing to 0.82.

Magnitude of Thermal Storage. The importance of the thermal storage fluxes can be inferred from their magnitude with respect to the other fluxes that are included in the surface energy balance equation (Eq. 1). The first example is for a single spring day, May 14 (DOY 135), which is plotted in Fig. 4 as rates (W/m2) averaged over 30-min periods. Figure 4 shows that the sensible and latent fluxes at this uniform forest site are approximately equal, and the thermal storage rate is relatively minor throughout the day. This can be seen more clearly by examining daily totals obtained by integrating the flux density rates over the day.

Energy totals for each 30-min period can be obtained in joules per square meter by multiplying each data point by the length (s) of its averaging period (1800 seconds in a 30-min periods), and then for convenience expressing the totals as MJ/m2 (1000 J/m2 = 1 MJ/m2). Totals are derived for DOY 135 and 257 from Tab. 1 A&B for thermal storage rate (G) and for net radiation (Q). Totals are presented in Tab. 2 for approximate daylight (Q>0), for approximate night (Q<0) and for the entire 24-hr period. The Table shows that the thermal storage total was about 10 % of net radiation for daylight in spring, and 7 % for the 24-hr day. The fall totals were about 8 % in daylight and 3 % for the full day. The thermal storage fluxes are illustrated for clear days in late spring and in late fall. The figures show that thermal storage is small, but not negligible.

Table 2. Net radiation and thermal storage totals at Hartheim for a clear day in spring and in fall. Units are MJ/m2 for the indicated period.

 Period Spring DOY 135Q Spring DOY 135G Fall DOY 257Q Fall DOY 257G daylight 19.9 -2.1 12.2 -1.0 night -2.0 0.8 -2.0 0.7 24-hr 17.9 -1.2 10.3 -0.3

The magnitude of thermal storage over seasonal periods is shown in Fig. 5. Flux totals are determined for each day, and then accumulated for the 157-day period May 11 through Oct 14 (DOY 133 through 288). This period is roughly equivalent to the annual growing season, and it contains a mix of clear, cloudy and rainy days. Total thermal storage at the HARTX site was G = -50 MJ/m2, which is a little less than 3% of the 1744 MJ/m2 of net radiation accumulated over the same period. The analysis stops just as winter is beginning, but we can expect that essentially all of the energy gain stored at the site on DOY 288 will be gone by midwinter.

The rate of thermal storage can be expressed as a linear function of environmental index temperatures through a thermal storage coefficient. The magnitude of the flux depends upon characteristics of the site. The storage component is affected by the volumetric heat capacity and thickness of the storage media, and length of the data averaging period. Thermal storage is less important when averaged over diurnal or annual cycles. Storage for surfaces covered with low, uniformly dense canopies is quite small, even for short periods. Magnitude of thermal storage can be large for certain surfaces (especially water and bare soil), and for tall forests with large amounts of biomass.

The computation of thermal storage rates was demonstrated with data obtained in a 12-m Scots pine plantation over the summer growing season. The thermal storage rates occasionally ranged up to 10% of net radiation averaged over 30-min periods, but were typically only a few percent over 24-hr periods. Essentially the same results were obtained with several different methods for defining the thermal storage coefficient and/or developing average temperatures in the storage media. Two independent sets of measurements were in good agreement at the site of the example forest. Energy gained and lost to the soil was the most important thermal storage component at the forest site Direct measurements of the soil thermal storage rate with heat flux plates at this site agreed well with a more elaborate method based upon definition of profiles of soil temperature and moisture.

LITERATURE CITED

Blanford, J.H., and L.W. Gay. 1992. Tests of a robust eddy correlation system for sensible heat flux. Theor. Appl. Climatol., 46:53-60.

Fritschen, L.J., and L.W. Gay. 1979. Environmental Instrumentation. Springer Verlag, New York. 216 pp.

Gay, L.W., Vogt, R., and A. Kessler. 1996. The May-October energy budget of a Scots pine plantation at Hartheim, Germany. Theor. Appl. Climatol., 46:79-94.

Jaeger, L., and A. Kessler. 1996. The HartX period May 1992, seen against the background of twenty years of energy balance climatology at the Hartheim pine plantation. Theor. Appl. Climatol., 46:9-21.

List, R.J. 1949. Smithsonian Meteorological Tables. 6th Rev.Ed., 6th reprint. Smithsonian Misc. Collections, Vol. 114. Smithsonian Inst. Press, Washington, DC. 527 pp.

Stewart, J.B. 1988. Modeling surface conductance of pine forest. Agr., For. Meteor. 43:19-35.

Tanner, C.B. 1960. Energy balance approach to evapotranspiration from crops. Soil Sci. Soc. Amer. Proc. 24:1-9.

Vogt, R., Ch. Bernhofer, L.W. Gay, L. Jaeger, and E. Parlow. 1996. The available energy over a Scots pine plantation: what's up for partitioning? Theor. Appl. Climatol. 46:23-32.