The surface energy balance model used in this study was established based on the principles of energy balance equations (Bornstein, 1986; Oke, 1987; Terjung, 1978) which were developed based on the assumption that the surface in the modeling area is homogeneous in radiative, thermal, and geometrical properties (Carlson et al., 1981; Outcalt, 1971; 1972). In case of a nonhomogeneous area, such as a hilly area, or a partially vegetated area, the surface has to be partitioned into small, homogeneous facets so that the equations can be applied to each uniform facet (Terjung, 1978). Surface partitioning has been a difficult task in surface energy balance model establishment. This difficulty became one of the obstacles for model development and application.
Since GRID and TIN structure in a GIS decompose surfaces into relatively homogeneous facets, the integration of a surface energy balance model with GRID or TIN will enhance the model's ability to simulate the spatial variation of surface temperature and heat fluxes in heterogeneous areas. The results from the model in a GIS can be easily classified or indexed with the spatial analysis functions of a GIS, and displayed spatially, which significantly improve the model's applicability.
This research is an attempt to embed an energy balance climate model into a GIS based on the existing GIS analysis capabilities and macro language functions (Goodchild, 1993).
where Q* is the net all-wave radiation which is the balance between the net shortwave radiation (K*) and net longwave radiation (L*). The net radiation input is balanced by three heat fluxes: sensible heat flux (QH), latent heat flux (QE), and subsurface heat flux (QG). The equation is solved by searching a unique equilibrium surface temperature that results in a balance.
On an opaque surface, there is no radiative transmission through the surface, therefore, with given albedo (a ) or absorptivity ( zetashort), the reflected shortwave radiation (Kup ) in Equation (2) can be derived based on incoming shortwave radiation (Kdown),
Kdown involves two major components, direct solar radiation (S) and diffuse solar radiation (D),
Direct solar radiation (S) can be represented by,
where S0 is direct radiation at normal incidences, and is a parameter determined by slope angle (phi), slope azimuth (OMEGA sl), zenith angle (Z), and solar azimuth (OMEGA s):
Between sunrise and sunset, the zenith angle is given by:
where PHI , delta , and h are latitude, declination angle, and hour angle. Declination angle (delta ) is a dependent variable of the date of the year and is derived from an array of values. The hour angle h defines the time of a day. By setting the zenith angle to 90, the sunrise and sunset hour angles (hrise, hset) can be derived with Equation (5b).
The solar azimuth angle OMEGA s is a function of latitude, declination angle, and zenith angle at a given time:
S0 in Equation (5) is determined by,
where I0 is the solar irradiance (I0 = 1370 Wm-2) (Kale, 1985), V represents radius vector which adjusts the distance between the earth and the sun, and is the transmission coefficient which is a function of atmospheric pressure (p), precipitable water vapor (w), optical air mass (m), and dust-haze factor (d) (Brooks 1959),
where p0 is sea level pressure (p0 = 101.3 kPa). The dust-haze factor d has a value ranging from 0.6 to 1.0 for "clear" conditions, and from 1.4 to 2.0 for polluted conditions (Terjung, 1978). The optical air mass m is adjusted by the elevation z, p is the pressure at elevation z. The precipitable water vapor w varies with vapor pressure, ea.
Diffuse solar radiation (D) comes from two sources: the diffuse radiation from the sky, and the reflected radiation from the environment. Because obtaining reflected radiation from the environment is a very complicated process, a separate study of the integration between this part of the model and a GIS is required. In this study, the model is simplified to assume that the diffuse radiation comes from the sky only. Since the equations in this study are applied to relatively small facets (cells in GRID and triangles or small irregular facets in TIN), when there is no drastic change of the terrain, this assumption could be true (the facets don't "see" their neighboring facets, therefore they don't receive the reflection from the neighbors). Under such an assumption, the diffuse radiation can be calculated by,
where Ss and Sh are direct radiation on the top of the atmosphere and on a horizontal surface at the bottom of the atmosphere (both can be derived from equation (5) and (6)), m is the optical air mass, and w is the precipitable water vapor. Equation (7) is based on an empirical approach developed by Terjung and Louie (1973) from observations in the Los Angeles area. Because this equation was developed for Los Angeles, adjustments may be necessary for other locations. An isotropic condition is assumed for the diffuse radiation from the sky.
All the variables in Equation (2) are either given or be solved yielding the shortwave radiation budget components for Equation (1).
where long is the longwave absorptivity of the surface.
The outgoing longwave radiation from a surface (Lup ) can be determined by the Stefan-Boltzmann law:
where is the Stefan-Boltzmann constant (5.67x10-8 Wm-2K-4), s is the emissivity of the surface given as a surface property, and Ts is the surface temperature in Kelvin.
The incoming longwave radiation comes from two sources: the sky and the environment. Similar to the diffuse radiation, it is assumed in this study that the incoming longwave radiation comes from the sky only:
where Ta is the shelter height air temperature and ea is the vapor pressure. This formula was developed by Brunt (1932), which, compared to other equations, was found to be more stable with temperature change than other methods (Jiménez et al, 1987).
All the variables required to solve longwave radiation are known except surface temperature (Ts) in Equation (9) which will be discussed later.
Sensible heat flux QH is expressed as a function of the temperature difference between the surface and the adjacent air and a heat transfer coefficient (hc)(Terjung and O'Rourke, 1980):
where u is the wind speed measured at shelter height (m s-1). Since Ta in Equation (11) is given by input data, the only unknown variable in Equation (11) is the surface temperature (Ts).
To derive latent heat flux (QE), a Bowen Ratio approach (Oke, 1987) is used. Bowen Ratio (beta ) is the ratio between sensible heat flux and latent heat flux. Typical values of beta range from less than 1 for very humid vegetated environments to greater than 10 for very dry desert environments (Oke, 1987). Based on Equations (11), and with given beta , QE can be presented as:
All the variables in Equation (12) are known except Ts .
The last term in Equation (1) is the subsurface heat flux (QG), the energy exchange between the surface and the subsurface. QG is determined by the subsurface temperature gradient and thermal conductivity of the substrate (k):
Tg is the subsurface temperature at a distance (d) below the surface and is given as input data. k is given based on surface material and moisture availability (Oke, 1987). The only unknown variable in Equation (13) is, again, the surface temperature Ts.
Surface temperature (Ts) is obtained by an interactive numerical approach (Gerald and Wheatley, 1994). The processes begin with an initial Ts value. Each term is solved based on the initial value. If the result exceeds zero, Ts is adjusted by a given interval and each function of Ts is recalculated. After every calculation, the interval is reduced by half. The process is repeated until the equation balances when the difference between both sides of the equation falls within a given tolerance level (for example, less than 1.0 Wm-2). At this point, the surface temperature is solved and so are L*, QH, QE, and QG.
Although the specific procedures in integrating the model with GRID and TIN are different (the former is based on raster procedures and the latter vector procedures), the principles of the integration are identical in both cases. Both GRID and TIN divide a continuous surface into uniform units, square cells in GRID and irregular facets in TIN (Burrough, 1986; ESRI, 1991a, 1991b). It is assumed that no variation occurs across a single unit surface, whether it is a square cell or an irregular triangle or polygon. This provides a foundation for the surface energy balance model because the energy balance can be sought for each individual unit which meets the assumption of homogeneity for the energy balance equations. When simulations are carried out for each and all the cells or facets that compose a nonhomogeneous surface, the spatial variation over the surface is derived.
The integrated model with GRID is illustrated by Figure 1.
The input data of the model are given by three different formats. The astronomical data (time, date, and latitude, etc.) and the non-surface dependent meteorological data (input air temperature, pressure, etc.) are given by assigning respective values to system wide variables in macro language blocks. The spatially varied input data (elevation and surface type) are given by GRID format. The gird for the surface type can be created based on classified surface types (eg. vegetated surface, bare soil surface, paved surface, etc.) or based on classified aerial photographs. The radiative and thermal properties associated with each surface type are given by attribute tables related to the surface grid. The initialization stage of the simulation involves further processing of the input data. The grid layers representing slope and aspect are created from the elevation layer with existing GRID functions. Incoming longwave radiation, diffuse radiation and the direct radiation at normal incidences are calculated based on input data by macro language blocks. An initial surface temperature grid is created based on either the subsurface temperature or the air temperature (the initial value will be replaced by simulated Ts values). After the initialization, simulations of the radiation terms and the heat fluxes terms are carried out based on the initial values with existing grid functions. The energy balance equations are solved for every grid in the study area. The energy balance for each of the cells is checked. The cells that reached the balance are recorded to an output layer and a mask is created to block out these cells for further processing. For the unbalanced cells, the initial surface temperatures are changed based on a given interval and the simulations are repeated until all the cells reach their balance. The results are then written to several output grid layers. From there, further classification and conversion to vector format are possible.
The integration of the model with TIN is illustrated by Figure 2.
The input component of the integrated model with TIN is identical to the GRID except that the spatially varied input data now are entered with polygon layers. At the initial stage, a TIN is generated from the elevation layer. Because TIN is not directly associated with any attribute tables, it has to be converted back to polygon layers (all the polygons are triangles) so that the slope and aspect information can be retained in attribute tables associated with the polygon layers. The slope and the aspect layers are then combined with the surface type layer by the existing vector analysis function of the GIS to create a polygon layer composed of uniform facets (no longer triangles because the combination with the surface type) in terms of its slope, aspect, and surface types. From this point on, almost all the simulations are carried out in attribute format using the facets as the basic units.
a. The model is supplied with spatially varied input data which enhance the model's capability of providing spatial patterns;
b. The cell or small polygon-based simulation unit ensures that the assumption of homogeneity for the energy balance equations is met more precisely; and
c. The spatial analysis function of the GIS enables the model's output to be displayed spatially and to be classified in various ways which greatly enhances the results analyses process and application.
The integrated model at this stage still has limitations. As discussed earlier, the model is simplified in some aspect to make the integration possible. One of the simplifications is the elimination of the radiative input from the environment - the neighboring cells or facets. Such a simplification will certainly affect the model's accuracy under certain conditions when longwave radiation and diffuse radiation from the environment are important factors. Since simulation of the environmental effects is a complicated process (Johnson and Watson, 1984; Oke, 1987; Steyn, 1980) existing GIS functions at current level cannot provide easy solutions. Further research is needed to study the possibility of this part of the integration.
Bornstein, R. D., (1986) Urban climate models: nature, limitations and applications. WMO No.652, pp 237-276.
Brooks, F. A., (1959) An Introduction to Physical Microclimatology Davis: Univ. of Calif.
Brunt, D., (1932) Notes on radiation in the atmosphere. Quart. J. Roy. Meteor. Soc., 58: 389-418.
Burrough, P.A., (1986) Principles of Geographical Information Systems for Land Resources Assessment. Oxford Univ. Press.
Carlson, T. N., Dodd, J. K., Benjamin, S. G. and Copper, J. N., (1981) Satellite estimation of the surface energy balance, moisture availability and thermal inertia. Jour. of Climate and Appl. Meteor., 20: 67-87.
Gerald, C. F. and Wheatley, P. O. (1994) Applied numerical analysis., 5th ed. Addison-Wesley Pub.
Goodchild, M. F., (1993) The state of GIS for environmental problem-solving. Environmental Modeling with GIS. Oxford Univ. Press. Pp 8-15.
ESRI, (1991a) Cell based modeling with GRID. ESRI
ESRI, (1991b) Surface modeling with TIN. ESRI
Jiménez, J. I. , Alados-Arboledas, L., Castro-Diez, Y. and Ballester, G., (1987) On the estimation of long-wave radiation flux from clear skies. Theor. Appl. Climatol., 38: 37-42.
Johnson G. T. and Watson I. D., (1984) The determination of view-factor in an urban canyon. Jour. of Climate and Appl. Meteor., 23: 329-335.
Kyle, H. L., Ardanuy, P. E. and Hurley, E. J., (1985) The status of the Nimbus-7 earth-radiation budget data set. Bull. Am. Meteor. Soc., 66: 1378-1387.
Oke, T. R., (1986) Urban climatology and the tropical city: an introduction. in Proc. of the technical conference of urban climatology and its applications with special regard to tropical areas, WMO No.652, pp1-25.
Oke, T. R., (1987) Boundary layer climates. 2nd ed., Methuen, N.Y.
Outcalt, S. I., (1971) A numerical surface climate simulator. Geog. Analysis, 3: 380-303.
Outcalt, S. I., (1972) The development and application of a simple digital surface-climate simulator. Jour. Appl. Meteor., 11: 629-636.
Steyn, D. G. (1980) The calculation of view factors from fish-eye lens photographs. Atmosphere-Ocean, 18: 254-258.
Terjung, W. H. and Louie, S. S., (1974) A climatic model of urban energy budgets. Geog. Analysis, 6: 341-367.
Terjung, W. H. and O'Rourke, P. A., (1978) An outline of boundary-layer climatology methods and analysis. Academic Publishing Service ASUCLA, Los Angeles.
Terjung, W. H. and O'Rourke, P. A., (1980) Energy exchanges in urban landscapes: selected climatic models. Publ. in Clim., Vol. XXXIII, C. W. Thornthwaite associates..
Lin Wu
Assistant Professor
Department of Geography and Anthropology
California State Polytechnic University, Pomona
3801 W Temple Ave
Pomona, CA91768
Telephone: (909)869-3578
Fax (909)869-3586
E-mail: LWU@CSUPomona.edu