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/m^{2})
of a vegetated surface can be estimated from 3 components: **G =
G**_{a}** + G**_{veg}** + G**_{s
}[** Eq. 1**], where G

The first step is to evaluate K_{a} and
K_{veg}. The storage coefficient **K**_{x}**
= (****r**** **_{x}**C**_{px}**)(h**_{Z}**/****D**** t)** [** Eq.
3**], where r is density (kg/m

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

K_{a} and K_{veg} can now be
applied to the observed previous and subsequent index
temperatures in Eq. 2 to calculate present values of G_{a}
and G_{veg}. The values are combined with observed G_{s}
in Eq. 1 to yield G (W/m^{2}).

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 (Q_{n}),
thermal storage flux (G), sensible heat flux (H) and latent heat
flux (LE). These fluxes are expressed in flux density units of
W/m^{2}. 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:

Q_{n} + 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 = Q_{n} + 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, Q_{n}
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/m^{2}) can be separated
into three components:

G (W/m^{2}) = G_{a} + G_{veg}
+ G_{s} (1)

where G_{a} is the rate of energy
storage within the layer of air between the soil surface and the
sensors above the canopy, G_{veg} is the rate of energy
storage within the mass of vegetation, and G_{s} 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:

G_{x} = K_{x}(T_{p} - T_{s}),
(2)

where K_{x} [W/m^{2}K] 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 (P_{o} =
101.33 kPa):

**K**_{x}
[W/m^{2}K] = (r _{x}C_{px})(h_{z}/D t) (3)

where subscript "x" is either
"a", "veg", or "s"; the volumetric
heat capacity [J/m^{3} K] is the product of the
appropriate density r [kg/m^{3}] and heat capacity C_{p}
[J/kg K]; h_{z} [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 K_{x} incorporate
values of h_{z} and D t and so the coefficients are site specific.

__Measuring Temperature Changes.__
Temperature changes (T_{p} - T_{s}, 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 G_{a} and G_{s},
and this is usually obtained from temperature profiles measured
with a number of sensors. With G_{veg}, 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 (G_{a}) and in the
vegetation (G_{veg}). This approximation method also
measured G_{s} 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 (Q_{n} + 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 K_{a}
can be obtained from Eq. 3. The specific heat of dry air is
essentially a constant at C_{p} = 1010 J/kg K with
ambient temperature and pressure. The variables h_{Z} 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/m^{3}] = 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 (P_{Z})can
be approximated from the International Standard Atmosphere (List,
1949, p. 268):

P_{Z}/P_{o} = ((288 -
0.0065Z)/288)^{5.256} (4b)

where Z [m] is the elevation of the measurement
site. For example, the ratio P_{Z}/P_{o} = 0.976
at 200 m, and 0.823 at 1600 m elevation. Since sea level pressure
is 101.33 kPa, Eq. 4b estimates P_{200} = 98.898 kPa and
P_{1600} = 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/m^{3}.

There are several possibilities for determining
the air temperature to be used in evaluating K_{a} 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 K_{a} and
G_{a} at that time step. A new K_{a} 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 T_{p} and T_{s.}

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/m^{3} at mean air
temperature of 20 ° C; C_{p} = 1010 J/kgK; h_{z} = 17 m;
and P/P_{o} = 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 K_{a}
= 5.63 W/m^{2}K for the HARTX elevation of 200 m. At a
higher site elevation, say Z = 1600 m on the Rio Grande of
central New Mexico, P_{Z}/P_{o} = 0.823 and K_{a}
= 4.75 W/m^{2}K at a mean air temperature of 20 C. The
value of K_{a} 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 K_{veg} 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/m^{2}. With an average
tree height of 12 m, this is equivalent to vegetation density of
0.84 kg/m^{3} over a layer thickness h_{Z} = 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 (C_{p} = 4190 J/kgK). The time
period for the central difference analysis remains at D t = 3600 s. These
values in Eq. 3 yields K_{veg} = 0.6(0.842)4190(12)/3600
= 7.06 W/m^{2}K at HARTX.

The vegetation storage coefficient K_{veg}
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 G_{s}. Use of soil heat flux
discs also eliminated the necessity for measuring soil
temperatures

CALCULATING THE STORAGE FLUX G

The storage flux G [W/m^{2}] is
calculated by summing the 3 components as given in Eq. 1. The
polarity of the storage fluxes is defined by the term T_{p}
- T_{s} 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/m^{2} in spring, and 500 W/m^{2} 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 G_{s}, for G_{s}+G_{a}
and for G = G_{s}+G_{a}+G_{veg}. 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/m^{2}
in spring, and -50 W/m^{2} in autumn. When averaged over
the entire day, G shows a gain of -14 W/m^{2} in the
spring and only -3 W/m^{2} 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 (G_{s}) is the
largest of the 3 storage components at the HARTX site, with
30-min peaks reaching -60 W/m^{2} in the spring and -26
W/m^{2} 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 G_{s }component. Daily sums
in the air layer (G_{a}) and in the vegetation (G_{veg})
tend toward zero. This provides a rationale for mean daily
studies to concentrate upon G_{s}, and place a lower
priority upon measurements of G_{a} and G_{veg}.,
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 (G_{s}) and
indexed estimates (G_{a} and G_{veg}) 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 (R^{2} = 0.87), s _{y.x} = ± 11 W/m^{2}),
and G' = 0.84G -1 in the fall (R^{2} = 0.70, s _{y.x}
= ± 11
W/m^{2}). 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/m^{2}, and R^{2}
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/m^{2}) 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/m^{2} (1000 J/m^{2} = 1 MJ/m^{2}).
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/m^{2} 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/m^{2}, which is
a little less than 3% of the 1744 MJ/m^{2} 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.

CONCLUDING COMMENTS

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.