SEBAL Surface Energy Balance Algorithms for Land Idaho Implementation
Advanced Training and Users Manual August, 2002 Version 1.0
Waters Consulting
University of Idaho
WaterWatch, Inc.
Nelson, British Columbia
Kimberly, Idaho
Wageningen, the Netherlands
Ralf Waters
Richard Allen Masahiro Tasumi Ricardo Trezza
Wim Bastiaanssen
Funded by a NASA EOSDIS/Synergy grant fromthe Raytheon Company through The Idaho Department of Water Resources
1
TABLE OF CONTENTS Page 4
List of Symbols List of Figures
6
Introduction
7
1. Background 2. Concept 3. Objectives 4. Who can use this manual?
The Theoretical Basis of SEBAL 1. Overview – How SEBAL computes evapotranspiration (ET) 2. Data requirements and data preparation A. Cloudy conditions B. Header file information C. Land-use map D. Weather data
9 11 11 12 12 12
The Operation of the SEBAL Model 1. Setting up the image 2. Acquiring header file information, weather data, and reference ET values 3. Surface radiation balance equation A. Surface albedo B. Incoming shortwave radiation C. Outgoing longwave radiation D. Choosing the “cold” and “hot” anchor pixels E. Incoming longwave radiation F. Solving the surface energy balance equation 4. Surface energy budget equation A. Soil heat flux B. Sensible heat flux C. Latent Heat Flux, λET, Instantaneous ET, and Fraction ETr 34
14 14 14 16 19 19 23 23 24
5. 24-hour evapotranspiration 6. Seasonal evapotranspiration
35 38
Appendices 1. 2. 3. 4. 5.
24 25
40
Acquiring header file information Land-use map and momentum roughness length Reference evapotranspiration calculation Preparing the srcinal satellite image for use in SEBAL Time issues around weather data and reference evapotranspiration 2
6. Tables for surface albedo computations 7. Methods for selecting “cold” and “hot” pixels 8. Iteration spreadsheet used for the atmospheric stability correction 9. SAVI and LAI for southern Idaho 10. G/Rn computations for water surfaces and snow 11. The stability of the atmosphere 12. The SEBAL Mountain Model
3
LIST OF SYMBOLS
cp air specific heat J/kg/K d zero plane displacement m de-s relative earth-sun distance dr inverse squared relative earth-sun distance ET ETr G H k L Lλ rah Rn RS ↓
actual evapotranspiration rate ference re evapotranspiration soil heat flux sensible heat flux von Karman’s constant Monin-Obukhuv length spectral radiance for bandλ aerodynamic resistance to heat transport net radiation flux incoming shortwave radiation
mm/hr mm/hr W/m2 W/m2 0.41 m W/m2/sr/µm s/m W/m2 W/m2
RL ↓
incoming longwave radiation
W/m2
RL ↑
outgoing longwave radiation
W/m2
corrected thermal radiance from the surface Rp path radiance in the 10.4 – 12.5 µm band Rsky narrow band downward thermal radiation for a clear sky Ta near surface air temperature K Tsurface temperature K s u wind velocity u* friction velocity z height zom momentum roughness length
W/m2/sr/µm
α αpath-
surface albedo albedo path radiance
-
albedo at top of atmosphere solar elevation angle aspect angle of the pixel
degrees radians
Rc
W/m2/sr/µm W/m2/sr/µm
m/s m/s m m
radiance
αtoa β γ
4
δ εο εΝΒ εa θ λ λET ρ ρλ σ τsw τΝΒ φ Ψh ψm ϖ ωλ
declination of the earth broad band surface emissivity narrow band surface emissivity atmospheric emissivity solar incidence angle latent heat of vaporization latent heat flux air density reflectivity for bandλ Stefan-Boltzmann constant shortwave transmissivity of air narrow band transmissivity of air latitude of the pixel stability correction for heat transport stability correction for momentum transport mountain wind speed weighting coefficient weighting coefficient for bandλ
5
radians degrees J/kg W/m2 kg/m3 5.67 × 10-8 W/m2/K4 radians -
LIST OF FIGURES Page Figure 1. Surface Energy Balance
9
Figure 2. Surface Radiation Balance
10
Figure 3. Flow Chart of the Net Surface Radiation Calculation
15
Figure 4. Flow Chart of the Iterative Process for the Calculation of Sensible Heat (H)
29
Figure 5. Plot of ET and ETrF vs. time
36
Figure 6. A Map of 24-hour Evapotranspiration
37
Figure 7. A Map of Seasonal Evapotranspiration
39
6
INTRODUCTION 1. BACKGROUND Land managers today are faced with multiple challenges. They must, as a result of expanding economic pressures, become advocates for the sustainable use of our natural resources while ensuring that the quality of the resource is maintained. The public and the courts are increasingly demanding objective and effective strategies for managing these exhaustible resources. Fresh water has become our most precious natural resource and the wise management of this resource is one of our greatest challenges. An understanding of the natural systems and the physical laws that govern each component of the hydrologic cycle is very important for the water resource manager. In particular, the evaporation processes from the various surfaces of the earth need to be understood in order to achieve a sustainable development of our water resources. Cropland irrigation is a major consumer of water in semiarid and arid regions and an efficient and reliable method for determining the consumptive use of water by crops is crucial for adequate water management. As fresh water becomes scarcer, competition for fresh water intensifies and better irrigation management will be required to achieve greater efficiency in the use of this valuable resource. The combination of the two separate processes whereby water is lost to the atmosphere due to evaporation from the soil and transpiration from vegetation is referred to as evapotranspiration (ET). The capability to predict levels of ET would be a valuable asset for water resource managers. ET is a good indicator of irrigation effectiveness and total water consumption from vegetation. Evapotranspiration information is useful for irrigation supply planning, water rights regulation, and river basin hydrologic studies.
2. CONCEPT Experimentally, the calculation of evapotranspiration can be made quite accurately using weighing lysimeters, eddy correlation techniques, and the Bowen ratio technique. These methods are limited, however, because they provide point values of ET for a specific location and fail to provide the ET on a regional scale. This limitation has motivated the development of using remotely sensed data from satellites to evaluate ET over vast areas. The major advantage of remote sensing is that ET can be computed without quantifying other complex hydrological processes. ET is highly variable in both space and time. It is variable in space due to the wide spatial variability of precipitation, hydraulic characteristics of soils, and vegetation types and densities. It is variable in time due to variability of climate. Satellite images provide an excellent means for determining and mapping the spatial and temporal structure of ET.
7
Remote sensing has great potential for improving irrigation management, along with other types of water management by providing ET estimations for large land surface areas using a minimal amount of ground data.
3. OBJECTIVES This manual explains a remote image-processing model for predicting ET termed SEBAL (Surface Energy Balance Algorithm for Land). SEBAL calculates ET through a series of computations that generate: net surface radiation, soil heat flux, and sensible heat flux to the air. By subtracting the soil heat flux and sensible heat flux from the net radiation at the surface we are left with a “residual” energy flux that is used for evapotranspiration (i.e. energy that is used to convert the liquid water into water vapor). This manual describes the theoretical basis of SEBAL using images from Landsat 5 and 7 satellites. However, the theory is independent of the satellite type and this manual could be applied to other satellite images if used with appropriate coefficients.
4. WHO CAN USE THIS MANUAL? The user of this manual should have background knowledge in hydrologic science or engineering and environmental physics (theories of mass and momentum transfer), as well as solar radiation physics. SEBAL was developed for use with ERDAS IMAGINE’s Model Maker tool and experience with this software is a prerequisite. As with any computer model, only knowledgeable users who understand the theory and who can consistently recognize the reliability of the computations should use SEBAL. As will be shown, it is also important that the user be familiar with the area being studied and has some “on the ground” knowledge of land use and topography for the area of interest. Following is a list of useful references, which provide a good background of the scientific methods used in SEBAL: 1. 2. 3.
4.
Monteith, J. L. and Unsworth, M. H. (1990). Principles of Environmental Physics, Second Edition, Butterworth Heinemann. ISBN 0-7131-2931- X Campbell, G. S. and Norman, J. M. (1998). An Introduction To Environmental Biophysics, Second Edition, Springer. ISBN 0-387-94937-2 Allen, R. G., Pereira, L., Raes, D., and Smith, M. (1998). Crop Evapotranspiration, Food and Agriculture Organization of the United Nations, Rome, Italy. ISBN 92-5-104219-5. 290 p. Allen, R.G. (2000). REF-ET: Reference Evapotranspiration Calculation Software for FAO and ASCE Standardized Equations, University of Idaho. www.kimberly.uidaho.edu/ref-et/
8
THE THEORETICAL BASIS OF SEBAL 1. OVERVIEW – HOW SEBAL COMPUTES EVAPOTRANSPIRATION (ET) In the SEBAL model, ET is computed from satellite images and weather data using the surface energy balance as illustrated in Figure 1. Since the satellite image provides information for the overpass time only, SEBAL computes an instantaneous ET flux for the image time. The ET flux is calculated for each pixel of the image as a “residual” of the surface energy budget equation:
λET = –RnG – H
(1)
where; λET is the latent heat flux (W/m2), Rn is the net radiation flux at the surface (W/m2), G is the soil heat flux (W/m2), and H is the sensible heat flux to the air (W/m2). The surface energy budget equation is further explained in part 4 of this section.
Figure 1. Surface Energy Balance
9
The net radiation flux at the surface (Rn) represents the actual radiant energy available at the surface. It is computed by subtracting all outgoing radiant fluxes from all incoming radiant fluxes (Figure 2). This is given in the surface radiation balance equation: Rn = RS↓ - α RS↓ + RL↓ - RL↑ - (1-εo)RL↓
(2)
where; RS↓ is the incoming shortwave radiation (W/m2), α is the surface albedo (dimensionless), RL↓ is the incoming longwave radiation (W/m2), RL↑ is the outgoing longwave radiation (W/m2), and εo is the surface thermal emissivity (dimensionless).
shortwave radiation
RS↓
α RS↓
(incident shortwave)
(reflected shortwave)
longwave radiation
RL↓
(1-εo)RL↓
(incident longwave)
(reflected longwave)
RL↑ (emitted longwave)
vegetation surface
Net surface radiation =gains – losses Rn = (1 - α) RS↓ + RL↓ - RL↑ - (1-εo)RL↓
Figure 2. Surface Radiation Balance
In Equation (2), the amount of shortwave radiation (R S↓) that remains available at the surface is a function of the surface albedo α ( ). Surface albedo is a reflection coefficient defined as the ratio of the reflected radiant flux to the incident radiant flux over the solar spectrum. It is calculated using satellite image information on spectral radiance for each satellite band. The incoming shortwave radiation (R S↓) is computed using the solar constant, the solar incidence angle, a relative earth-sun distance, and a computed atmospheric transmissivity. The incoming longwave radiation (R L↓) is computed using a modified Stefan-Boltzmann equation with atmospheric transmissivity and a selected L↑) is computed using the surface reference temperature. Outgoing longwave radiation (R Stefan-Boltzmann equation with a calculated surface emissivity and surface temperature. Surface temperatures are computed from satellite image information on thermal radiance.
10
The surface emissivity is the ratio of the actual radiation emitted by a surface to that emitted by a black body at the same surface temperature. In SEBAL, emissivity is computed as a function of a vegetation index. The final term in Equation (2), (1εo)RL↓, represents the fraction of incoming longwave radiation that is lost from the surface due to reflection. The radiation balance equation is explained in detail in part 3 of this section. In Equation (1), the soil heat flux (G) and sensible heat flux (H) are subtracted from the net radiation flux at the surface (Rn) to compute the “residual” energy available for evapotranspiration (λET). Soil heat flux is empirically calculated using vegetation indices, surface temperature, and surface albedo. Sensible heat flux is computed using wind speed observations, estimated surface roughness, and surface to air temperature differences. SEBAL uses an iterative process to correct for atmospheric instability due to the buoyancy effects of surface heating. Once the latent heat flux (λET) is computed for each pixel, an equivalent amount of instantaneous ET (mm/hr) is readily calculated by dividing by the latent heat of vaporization (λ). These values are then extrapolated using a ratio of ET to reference crop ET to obtain daily or seasonal levels of ET. Reference crop ET, termed ET r, is the ET rate expected from a well-defined surface of full-cover alfalfa or clipped grass and is computed in the SEBAL process using ground weather data (see Appendix 3). SEBAL can compute ET for flat, agricultural areas with the most accuracy and confidence. Appendix 12 describes additional SEBAL applications for mountain areas where elevation, slope, and aspect need to be considered.
2. DATA REQUIREMENTS AND DATA PREPARATION SEBAL requires a satellite image and some weather data. A land-use map for the area of interest is also helpful. This manual applies specifically to images produced from the Landsat satellites. Landsat Thematic Mapper (TM) bands 1-5 and 7 provide data for the visible and near infrared bands. The pixel size for these bands is 30m by 30m. TM band 6 provides data for longwave (thermal) radiation. The pixel size for this band is 60m by 60m for Landsat 7 and 120m by 120m for Landsat 5. We recommend using images from Landsat 7 since there is more error possible with recent Landsat 5 images due to sensor degradation. Some Landsat 5 satellite constants that are used in the computations are updated to 1986 and may be in error for 2002 without adjustment. A. Cloudy Conditions It is important that the image used is for a totally clear sky. ET cannot be computed for cloud covered land surfaces, because even a thin layer of cloud can considerably drop the thermal band readings and cause large errors in calculation of sensible heat fluxes. Therefore, all satellite images should be thoroughly checked for the occurrence of cloud cover and if found, these areas should be masked out and dealt with individually. In order to detect clouds in an image, view the thermal band using a range of colors to differentiate temperature values. Clouds will then show up as uniquely colored areas in the image and
11
can be flagged. Masking can be done following the processing of ET, but is an important final step. B. Header File Information TM images are generally created with an associated “header” file. The header file for the satellite image is a relatively small file that contains important information for the SEBAL process. The following information must be obtained from the header file for entry into SEBAL: • The satellite overpass date and time • The latitude and longitude of the center of the image
sun elevation angle (β) at the overpass time •• The Gain and Bias levels for each band The satellite overpass time is expressed as Greenwich Mean Time (GMT) and must be converted to local time (see Appendix 5). For Landsat 7, the header file includes the gains and biases for bands 1-5 and 7. These parameters are used to convert digital numbers (DN’s) in the srcinal files into energy units. For band 6, the thermal band, both high and low gain images are supplied and the user must decide which one to use for SEBAL. We recommend using the low gain image, which yields slightly lower resolution, but is less likely to suffer from saturation . When the low gain image is selected, use the gain and bias information for low gain. If the header file is missing, the user must use the internet to find gain and bias data for the image. Appendix 1 gives information about acquiring header file information. C. Land-Use Map
A land-use map is not a requirement for SEBAL but is highly recommended since it is useful in estimating the surface roughness parameter (z om). The land-use map divides the area of interest into various general classes of land-use such as agriculture, city, water, desert, forest, grassland, etc. Appendix 2 lists the general land-use types used in the Idaho application of SEBAL. D. Weather Data The following weather data are recommended for processing SEBAL: • Wind speed (required, hourly data is preferred) • Precipitation (daily data are recommended for a several week period prior to the image)
•
In addition to wind speed data, the following weather data are recommended for calculating hourly and daily reference ET. If the reference ET data are available for the image area, then the following data are not needed. Humidity (hourly data such as vapor pressure or dew point temperature)
• Solar radiation (hourly is preferred) • Air Temperature (hourly is preferred)
12
• Reference evapotranspiration, ETr (hourly is preferred, see Appendix 3 for calculation). The alfalfa ETr is preferred over the clipped grass ETr since it more nearly represents the upper limit of ET from well-watered vegetation.
The wind speed (u) at the time of the satellite overpass is required for the computation of sensible heat flux (H) and for the ETr calculation. Precipitation data is used to evaluate the general “wetness” of areas that have received rain within four or five days before the image date. Humidity data are necessary for the ET r calculation. Solar radiation data are useful for the estimation of the cloudiness of the image and for adjusting the atmospheric transmissivity (τsw). Reference evapotranspiration (ETr) is the estimated ET for a wellwatered reference crop, usually alfalfa. It is used to compute H at the “cold” anchor pixel and to compute the ETr fraction (ETrF) that is used to predict 24-hour and seasonal ET. Weather conditions might not be uniform within a satellite image. If the area of interest has a variety of terrains and land-uses such as mountains, valleys, agriculture, desert, and grassland, the weather conditions may vary significantly. If the primary interest is ET from agricultural fields, one should use a weather station located in an agricultural area. In cases of widely varying terrain or land-use, one should consider using data from two or more weather stations. It is important that the weather data come froma location that is near (within 50 km) the locations of the “hot” and “cold” anchorpixels (described later).
13
THE OPERATION OF THE SEBAL MODEL
1. SETTING UP THE IMAGE a nd ACQUIRING HEADER FILE INFORMATION First, a satellite image for the area of interest is required. This manual uses Landsat 5 TM images or Landsat 7 ETM+ images. Next, the srcinal satellite image must be georectified so that the final can ET “map” is accurately to earth coordinates. Theparty such as georectification be done by the userreferenced or the image can be sent to a third Earth Satellite Corporation. It is recommended that the georectified image besaved in the GeoTIF format and that nearest neighbor resampling be used during georectification to preserve spectral information. The geo-rectified disk is used with the ERDAS IMAGINE software to produce a layered spectral band image necessary for the SEBAL procedure. Once a layered image is produced, a subset image can be made if a smaller area of interest is desired (Appendix 4 contains the ERDAS steps required for this image set up). Then, header file information is obtained as explained in Appendix 1. Original images for which the header file information is most accessible are Landsat 7 ETM+ in the NLAPS or FAST format.
2. WEATHER DATA and REFERENCE ET VALUES Gather weather data from a relevant weather station (see section 2D above) and organize it in a table. Now, compute the reference ET as described in Appendix 3 or using other means. The results are hourly estimates of ET r for the day of the image. Now, time issues regarding how the hourly weather data were gathered and what time period the data represent need to be examined (see Appendix 5).
3. SURFACE RADIATION BALANCE EQUATION The first step in the SEBAL procedure is to compute the net surface radiation flux (R n) using the surface radiation balance equation: Rn = (1-α)RS↓ + RL↓ - RL↑ - (1-εo)RL↓
(2)
This is accomplished in a series of steps using the ERDAS Model Maker tool to compute the terms in the above equation. A flow chart of the process is shown in Figure 3:
14
Figure 3. Flow Chart of the Net Surface Radiation Computation
Rn = (1-α)RS↓ + RL↓ - RL↑ - (1-ε0)RL↓ model F09
Surface albedo
Incoming shortwave
Outgoing longwave
Incoming longwave
α
RS↓
RL↑
RL↓
model F04
s readsheet
model F08
Albedo-top of atmosphere
s readsheet
Surface temperature
toa
Surface emissivities
model F03
&
TS model F07
model F06
Reflectivity
ρλ model F02
NDVI SAVI LAI
Spectral radiance
model F05
L model F01
The computer model number used for each computation is given along with the variable name. The computation steps begin at the bottom of the figure with model F01 and continue upward to model F09 for the computation of R n. The two terms RS↓ and RL↓ are computed with a spreadsheet or a calculator rather than the Model Maker tool. 15
The terms of Equation (2) and how they are computed are now explained. A. Surface Albedo ( ) Surface albedo is defined as the ratio of the reflected radiation to the incident shortwave radiation. It is computed in SEBAL through the following steps: 1. The spectral radiance for each band (Lλ) is computed (model F01). This is the outgoing radiation energy of the band observed at the top of the atmosphere by the satellite. It is calculated using the following equation given for Landsat 5 and 7:
LMAX − LMIN × ( DN − QCALMIN ) + LMIN Lλ = QCALMAX − QCALMIN
(3)
where; DN is the digital number of each pixel, LMAXand LMIN are calibration constants, QCALMAX and QCALMIN are the highest and lowest range of values for rescaled radiance in DN. The units for Lλ are W/m2/sr/µm. For Landsat 5, QCALMAX = 255 and QCALMIN = 0 so that Equation (3) becomes:
LMAX − LMIN × DN + LMIN 255
Lλ =
(4)
Table 6.1 in Appendix 6 gives the values for LMAX and LMIN. For Landsat 7 with header file data on gains and biases, a simpler equation for λLis given:
Lλ = (Gain × DN) + Bias
(5)
where; Gain and Bias refer to the values given in the header file. The low gain data for band 6 is recommended. If header file data is not available for Landsat 7, Equation (3) must be used, where QCALMAX is 255 and QCALMIN is 0 or 1. For the NLAPS generating system, QCALMIN is 0. For the LPGS generating system, QCALMIN is 1. Table 6.2 in Appendix 6 gives LMIN and LMAX values for low and high gain. For Landsat 7, the gain for a band can shift from low to high or vice-versa for any given row. Therefore, it is important to confirm, for each scene (path/row), the specific gain of each band. The header file on the srcinal image CD or obtained from the internet provides information on whether to use low or high gain for a given band (see Appendix 1). To complete this first step, model F01_radiance, is loaded and run. Values for LMIN and LMAX, or Gains and Biases are entered into the tables of the model along with the subset image. The model is saved along with the output file for radiance (L λ).
16
2. The reflectivity for each band (ρλ) is computed (model F02). The reflectivity of a surface is defined as the ratio of the reflected radiation flux to the incident radiation flux. It is computed using the following equation given for Landsat images:
ρλ =
π ⋅ Lλ
(6)
ESUN λ ⋅ cos θ ⋅ d r
where; Lλ is the spectral radiance for each band computed in model F01, ESUN λ is the 2 mean solar exo-atmospheric irradiance for each band (W/m /µm), cos θ is the cosine of the solar incidence angle (from nadir), and d r is the inverse squared relative earth-sun distance. Values for ESUNλ are given in Table 6.3, Appendix 6. Cosineθ is computed using the header file data on sun elevation angle β( ) where θ = (900 - β). The term dr is defined as 1/de-s2 where de-s is the relative distance between the earth and the sun in astronomical units. dr is computed using the following equation by Duffie and Beckman (1980)1, also given in FAO 56 paper:Crop Evapotranspiration (Allen et al., 1998): dr
= 1 + 0.033 cos DOY
2π
365
(7)
where; DOY is the sequential day of the year (given in Table 6.5, Appendix 6), and the angle (DOY × 2 /365)is in radians. Values for dr range from 0.97 to 1.03 and are dimensionless. To complete step 2, model F02_reflectivity, is loaded and run. Values for ESUN λ are entered into the table of the model. (Since band 6 is not shortwave, it is not used in this computation and is given a dummy value of 1 for ESUN). The constants rdand cos θ are also entered along with the radiance file computed in model F01. Again, the output file for reflectivity (ρλ) and the model are saved. 3. The albedo at the top of the atmosphere α ( toa) is computed (model F03). This is the albedo unadjusted for atmospheric transmissivity and is computed as follows:
αtoa = Σ (ωλ × ρλ)
(8)
where; ρλ is the reflectivity computed in step 2 andωλ is a weighting coefficient for each band computed as follows:
ωλ =
ESUN λ
∑ ESUN
1
(9) λ
Duffie, J.A. and W.A. Beckman. 1980. “Solar Engineering of Thermal Processes.” John Wiley and Sons, New York, p 1-109.
17
Values for the weighting coefficient,ωλ, are given in Table 6.4, Appendix 6, for Landsat images. To complete step 3, model F03_albedo_toa, is loaded and run. Values for ωλ are entered into the table of the model with a dummy value of zero for band 6. The file for reflectivity (ρλ) computed in model F02 is also input. The model and output file are saved. 4. The final step is to compute the surface albedo (model F04). Surface albedo is computed by correcting theαtoa for atmospheric transmissivity:
α =
α toa − α path _ radiance τ sw
2
(10)
where; αpath_radiance is the average portion of the incoming solar radiation across all bands that is back-scattered to the satellite before it reaches the earths surface, and τsw is the atmospheric transmissivity. Values for αpath_radiance range between 0.025 and 0.04 and for SEBAL we recommend a value of 0.03 based on Bastiaanssen (2000). Atmospheric transmissivity is defined as the fraction of incident radiation that is transmitted by the atmosphere and it represents the effects of absorption and reflection occurring within the atmosphere. This effect occurs to incoming radiation and to outgoing radiation and is thus squared in Equation (10). τsw includes transmissivity of both direct solar beam radiation and diffuse (scattered) radiation to the surface. We calculateτsw assuming clear sky and relatively dry conditions using an elevation-based relationship from FAO-56:
τsw = 0.75 + 2 × z10-5 ×
(11)
where; z is the elevation above sea level (m). This elevation should best represent the area of interest, such as the elevation of the relevant weather station. To complete step 4, model F04_surface_albedo, is loaded and run. The file for αtoa from model F03 is entered along with the constants forαpath_radiance and τsw. The model and output file are saved. The surface albedo for all pixels has now been computed and these values can be checked to see if they are realistic. Use the ERDAS Viewer tool to view the surface albedo image and values. Compare the values for various known surfaces with typical albedo values given in Table 6.6, Appendix 6. For conditions having high levels of atmospheric aerosols or dust, the user may wish to predict α separately for each individual band using an atmospheric radiation transfer model such as Modtran along withradiosonde data for the location andimage date. In most situations, however, because of the design of SEBAL, where the energy balance
18
and ET are specified for the two extreme conditions (maximum ET and minimum ET), some error in estimated abledo caused by uncertainty in values for transmissivity and αpath_radiance, cause little error in the final ET computation. Therefore, for general application, Equation (10) is adequate. B. Incoming Shortwave Radiation (R S↓) Incoming shortwave radiation is the direct and diffuse solar radiation flux that actually reaches the earth’s surface (W/m2). It is calculated, assuming clear sky conditions, as a constant for the image time using:
Rs↓ = Gsc × cos θ × dr ×
τsw
(12)
2
where; Gsc is the solar constant (1367 W/m ), cos θ is the cosine of the solar incidence angle as in step 2 above, dr is the inverse squared relative earth-sun distance, andτsw is the atmospheric transmissivity. This calculation is done with a spreadsheet or calculator. Values for Rs↓ can range from 200 – 1000 W/m2 depending on the time and location of the image. C. Outgoing Longwave Radiation (RL↑) The outgoing longwave radiation is the thermal radiation flux emitted from the earth’s surface to the atmosphere (W/m2). It is computed in SEBAL through the following steps (see the flow chart, Figure 3): 1. Three commonly used vegetation indices are computed (model F05). Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), and Leaf Area Index (LAI) are computed using the reflectivity values found in model F02. Any
one of these indices can be used to predict various characteristics of vegetation, according to preferences of the user. The NDVI is the ratio of the differences in reflectivities for the near-infrared bandρ4() and the red band (ρ3) to their sum: NDVI = (ρ4 − ρ3) / (ρ4 + ρ3)
(13)
where; ρ4 and ρ3 are reflectivities for bands 4 and 3 and are found in the output file of model F02. The NDVI is a sensitive indicator of the amount and condition of green vegetation. Values for NDVI range between -1 and +1. Green surfaces have a NDVI between 0 and 1 and water and cloud are usually less than zero. The SAVI is an index that attempts to “subtract” the effects of background soil from NDVI so that impacts of soil wetness are reduced in the index. It is computed as: SAVI = (1 + L) (ρ4 − ρ3) / (L + ρ4 + ρ3)
(14)
where; L is a constant for SAVI. If L is zero, SAVI becomes equal to NDVI. A value of 0.5 frequently appears in the literature for L. However, a value of 0.1 is used to better 19
represent soils of southern Idaho. The value for L can be derived from analysis of multiple images where vegetation does not change, but surface soil moisture does (Appendix 9). The LAI is the ratio of the total area of all leaves on a plant to the ground area represented by the plant. It is an indicator of biomass and canopy resistance. LAI is computed for southern Idaho using the following empirical equation:
0.69 − SAVI ID 0.59
ln
LAI = −
0.91
(15)
where; SAVIID is the SAVI for southern Idaho calculated from Equation (14) using a value of 0.1 for L. The maximum value for LAI is 6.0, which corresponds to a maximum SAVI of 0.687. Beyond SAVI = 0.687, thevalue for SAVI “saturates” with increasing LAI and does not change significantly. Users must realize that the relationship between SAVI and LAI can vary with location and crop type. Equation (15) is for southern Idaho. See Appendix 9 for more details of how these equations worked in southern Idaho. To complete this step, model F05_vegetation_indicies, is loaded and run. The file for reflectivity (ρλ) is entered and the three output files saved along with the model. 2. Surface emissivity (ε) is the ratio of the thermal energy radiated by the surface to the thermal energy radiated by a blackbody at the same temperature. Two surface emissivities are used in SEBAL. The first is an emissivity representing surface behavior for thermal emission in the relatively narrow band 6 of Landsat (10.4 to 12.5 µm), expressed as εΝΒ. The second is an emissivity representing surface behavior for
thermal emission in the broad thermal spectrum (6 to 14 µm), expressed as εο. εΝΒ is used in calculation of surface temperature (Ts) and εο is used later on to calculate total long wave radiation emission from the surface. The surface emissivities are computed in model F06 using the following empirical equations, where NDVI > 0:
εNB = 0.97 + 0.0033 LAI;
for LAI < 3
(16a)
= ε00.95 + 0.01 LAI;
for LAI < 3
(16b)
and εΝΒ = 0.98 and ε0 = 0.98 when LAI ≥ 3. For water and snow we use “filters” in the model to set the value ofε
• For water; NDVI < 0 and α < 0.47, • For snow; NDVI < 0 and α ≥ 0.47,
20
ΝΒ εΝΒ = 0.99 and εo = 0.985 εΝΒ = 0.99 and εo = 0.985
and εo:
To complete step 2, model F06_surface_emissivity, is loaded and run. The files for NDVI, LAI and surface albedo (α) are input and the output file and model are saved. 3. The corrected thermal radiance from the surface (R c) is calculated following Wukelic et al. (1989)2 as: L6 − R p Rc = − (1 − ε NB ) R sky (17) τ NB 2
where; L6 is the spectral radiance of band 6 from model F01 (W/m/sr/µm), Rp is the path radiance in the 10.4 – 12.5 µm band (W/m2/sr/µm), Rsky is the narrow band 2 downward thermal radiation for a clear sky (W/m /sr/µm), and τNB is the narrow band 2 transmissivity of air (10.4 – 12.5 µm). Units for R c are W/m /sr/µm.
The corrected thermal radiance (Rc) is the actual radiance emitted from the surface whereas L6 is the radiance that the satellite “sees”. In between the surface and the satellite two things occur. First, some of the emitted radiation is intercepted by the atmosphere (transmissivity). Second, some thermal radiation is emitted by the atmosphere in the direction of the satellite (path radiance) and the satellite “thinks” that this is from the surface is confused as being emitted from the surface. Values for Rp and τNB require the use of an atmospheric radiation transfer simulation model such as Modtran and radiosonde profiles representing the image and date. In the absence of values for these terms, they can be ignored by setting R p = 0 and τNB = 1 and the Rsky term can also be ignoredby setting Rsky = 0. This converts Rc into an uncorrected radiance (L6). Fortunately, the effects of the three parameters on R c are largely self-canceling. However, the result of no correction to 6L will be a general underestimation of surface temperature (Ts) by up to about 5o C for warmer portions of an image. However, because SEBAL tailors the “dT” function, described later, around the Ts data calculated for an image, the impact on final ET values is generally small, especially for areas having low and high values of ET. Errors may be larger for midrange ET values, but generally less than a few percent. The R sky term is small enough that it can probably be ignored all of the time. If calculated, it can be based on an Idso-Jackson style of empirical formula as applied by Wukelic et al. (1989): R sky
2 Ta ] = (1.807 x 10 −10 ) Ta4 1 − 0.26 exp − 7.77 x 10 − 4 [273.15 −(18)
where; Ta is the near surface air temperature at the time of the image (K).
2
Wukelic, G.E., D.E. Gibbons, et al., (1989). Radiometric calibration of Landsat thematic mapper thermal band. Remote Sensing of Environment 28:339-347.
21
The corrected thermal radiance (Rc) is computed in model F07 along with the surface temperature as shown below. 4. The surface temperature (Ts) is computed in model F07 using the following modified Plank equation: Ts
=
K2
(19)
ε K ln NB 1 + 1 R c
where; Rc is the corrected thermal radiance from the surface using 6L from model F01. K1 and K2 are constants for Landsat images (Table 1). Units for R c must be the same as those for K1. Table 1. Constants for Equation 19 for Landsat 5 TM in mW/cm2/sr/µm (Markham and Barker, 1986), and for Landsat 7 ETM+ in W/m2/sr/µm (Landsat 7 Science User Data Handbook Chap.11, 2002) K1 Landsat5 TM Band6 Landsat7 ETM+ Band6
607.76 666.09
K2 1260.56 1282.71
To complete step 4, modelF07_surface_temp, is loaded and run. Values for K1 and K2 and Rp, τNB, and Rsky are input along with the files for spectral radiance (L λ) and narrow band emissivity (εΝΒ). The output file and model are saved. (0, 1, and 0 are default values for Rp, τNB, and Rsky). 5. The outgoing longwave radiation (RL↑) is computed (model F08). This is computed using the Stefan-Boltzmann equation:
RL↑ = εo × σ × Ts4
(20)
where; εο is the “broad-band” surface emissivity (dimensionless),σ is the StefanBoltzmann constant (5.67 × 10-8 W/m2/K4), and Ts is the surface temperature (K). To complete this step, model F08_longwave_out, is loaded and run. Files for the surface emissivity (εo) and surface temperature (TS) are entered and the output file and model saved. Values for RL↑ can range from 200 – 700 W/m2 depending on the location and time of the image.
22
D. Choosing the “Hot” and “Cold” Pixels The SEBAL process utilizes two “anchor” pixels to fix boundary conditions for the energy balance. These are the “hot” and “cold” pixels and are located in the area of interest. The “cold” pixel is selected as a wet, well-irrigated crop surface having full ground cover by vegetation. The surface temperature and near-surface air temperature are assumed to be similar at this pixel. The “hot” pixel is selected as a dry, bare agricultural field where ET is assumed to be zero. Both of these “anchor” pixels should be located in large and ×60m for Landsat 7 homogeneous areas that contain more than one band 6 pixel (i.e. 60m and 120m×120m for Landsat 4&5).
The selection in of SEBAL these “anchor” pixels skill and of practice. Thepixels. qualityAppendix of the ET7 computations depends on a requires careful selection these two contains a guide for choosing the “hot” and “cold” pixels. Once selected, the surface temperature and x and y coordinates for these two pixels are recorded for future computations. E. Incoming Longwave Radiation (R L↓) The incoming longwave radiation is the downward thermal radiation flux from the atmosphere (W/m2). It is computed using the Stefan-Boltzmann equation:
RL↓ = εa × σ × Ta4
(21)
where; εa is the atmospheric emissivity (dimensionless),σ is the Stefan-Boltzmann constant (5.67 × 10-8 W/m2/K4), and Ta is the near surface air temperature (K). The following empirical equation for εa by Bastiaanssen (1995)3 is applied using data from alfalfa fields in Idaho:
εa = 0.85 × (-ln τsw).09
(22)
where; τsw is the atmospheric transmissivity calculated fromEquation (11). The srcinal coefficients by Bastiaanssen (1995), derived for western Egypt, wereεa = 1.08 × (-ln τsw).265. Substituting Equation (21) into Equation (20) and using T cold from our “cold” pixel for Ta yields the following equation: RL↓ = 0.85 × (-ln
τsw).09 × σ × Tcold4
(23)
This computation is done using a spreadsheet or calculator. Values for LR↓ can range from 200 – 500 W/m2, depending on the location and time of image. Some automatic weather stations include a net radiometer that measures all incoming and outgoing short and longwave radiation. Measured RL↓ data can be used to check computed values. 3
Bastiaanssen, W.G.M. 1995. Regionalization of surface flux densities and moisture indicators in composite terrain: A remote sensing approach under clear skies in Mediterranean climates. Ph.D. dissertation, CIP Data Koninklijke Bibliotheek, Den Haag, The Netherlands. 273 p.
23
F. Solving the Surface Radiation Balance Equation for R n The net surface radiation flux (Rn) is now computed using Equation (2). Load and run model F09_net_radiation. Files for surface albedo α( ), outgoing longwave radiation (RL↑), and surface emissivity (εo) are input along with the incoming shortwave radiation (R S↓) and the incoming longwave radiation (RL↓). The output file for Rn and the model are saved. Values for Rn can range from 100 – 700 W/m2, depending on the surface. This completes the first step of the SEBAL procedure.
4. SURFACE ENERGY BUDGET EQUATION The second step of the SEBAL procedure is to compute the terms G and H of the surface energy budget equation. The net radiation flux (R n) is the net amount of radiant energy that is available at the surface for warming the soil, warming the air, or evaporating soil moisture. This is written as the surface energy budget equation: Rn = G + H + λET
(23)
where; Rn is the net radiation at the surface (W/m2), G is the soil heat flux (W/m2), H is the sensible heat flux to the air (W/m2), and λET is the latent heat flux (W/m2). SEBAL computes λET as a “residual” of net radiant energy after soil heat flux and sensible heat flux are subtracted:
λET = –RnG – H
(1)
This equation is solved through the following steps using the ERDAS Model Maker tool. A. Soil Heat Flux (G) Soil heat flux is the rate of heat storage into the soil and vegetation due to conduction. SEBAL first computes the ratio G/Rn using the following empirical equation developed by Bastiaanssen(2000) representing values near midday:
G/Rn = Ts/α (0.0038
+ 0.0074 2)(1 - .98NDVI4)
(24)
where; Ts is the surface temperature (o C), α is the surface albedo, and NDVI is the Normalized Difference Vegetation Index. G is then readily calculated by multiplying G/R n by the value for Rn computed in model F09. Soil heat flux is a difficult term to evaluate and care should be used in its computation. One must understand the area of interest in order to evaluate the accuracy of Equation (24). Values of G should be checked against actual measurements on the ground. Land classification and soil type will affect the value of G and a land-use map is valuable for identifying the various surface types (Appendix 2). Equation (24) predicts G measured for
24
irrigated crops near Kimberly, Idaho quite accurately (M.Tasumi and R.Allen, 2002, pers. commun.). In the Idaho applications of SEBAL, values of G/R n for water and snow are assigned as follows (Appendix 10), representing values near midday: • If NDVI < 0; assume surface is water; G/Rn = 0.5 • If Ts < 4 oC and α > 0.45; assume surface is snow; G/R n = 0.5 The estimation of G/Rn for deep, clear lakes is complex. It may be large in early summer when the lake is colder than the air and smaller in the fall when the lake is warmer. Any n carefully significant water bodies in the area of interest should be flagged and G/R analyzed (Appendix 10). G/Rn for turbid or shallow water bodies will be less than 0.5 due to absorption of short-wave radiation more near the water surface.
To complete this step, model F10_soil_heat_flux is loaded and run. Files for NDVI, surface temperature, surface albedo, and Rn are entered and the output files and model are saved. View the image for G/Rn and check the values for various land types (Table 2). Table 2. Estimates of G/Rn for various surfaces
Surface type Deep, clear water Snow Desert Agriculture Bare soil
G/Rn 0.5 0.5 0.2 – 0.4 0.05 – 0.15 0.2 – 0.4
Full cover alfalfa Rock
0.04 0.2 – 0.6
B. Sensible Heat Flux (H) Sensible heat flux is the rate of heat loss to the air by convection and conduction, due to a temperature difference. It is computed using the following equation for heat transport: H = (ρ × cp × dT) / rah (25)
where; ρ is air density (kg/m3), cp is air specific heat (1004 J/kg/K), dT (K) is the temperature difference (T1 – T2) between two heights (z1 and z2), and rah is the aerodynamic resistance to heat transport (s/m). The sensible heat flux (H) is a function of the temperature gradient, surface roughness, and wind speed. Equation (25) is difficult to solve because there are two unknowns,ahr and dT. To facilitate this computation, we utilize the two “anchor” pixels (where reliable values for H can be predicted and a dT estimated) and the wind speed at a given height. The aerodynamic resistance to heat transport (rah) is computed for neutral stability as:
25
rah
=
z ln 2 z1
(26)
u* × k
where; z1 and z2 are heights in meters above the zero plane displacement (d) of the vegetation, u* is the friction velocity (m/s) which quantifies the turbulent velocity fluctuations in the air, and k is von Karman’s constant (0.41). The friction velocity (u *) is computed using the logarithmic wind law for neutral atmospheric conditions:
u* =
ku x
z ln x z om
(27)
where; k is von Karman’s constant, ux is the wind speed (m/s) at height zx, and zom is the momentum roughness length (m). zom is a measure of the form drag and skin friction for the layer of air that interacts with the surface. Figure 4 shows how the following steps are used to compute sensible heat flux (H): 1. The friction velocity (u*) at the weather station is calculated for neutral atmospheric conditions using Equation (27). (An explanation of neutral conditions is given in Appendix 11) The assumption of neutral conditions for the weather station is reasonable if the weather station is surrounded by well-watered vegetation. The calculation of* u requires a wind speed measurement (ux) at a known height (zx) at the time of the
satellite image. This measurement is taken from the weather data that was set up earlier in the SEBAL process. The momentum roughness length (z om) is empirically estimated from the average vegetation height around the weather station using the following equation (Brutsaert, 1982): = zom 0.12h
(28)
where; h is the vegetation height (m). (At the Kimberly weather station in southern Idaho, h=0.3m was chosen to represent the weather station surrounded by a small 50 m expanse of 0.15 m grass with other taller agricultural crops upwind of that.) This calculation is done on a calculator or spreadsheet. 2. The wind speed at a height above the weather station where one can assume no effect from the surface roughness is calculated. This height is referred to as the “blending height” and 200 meters is used. u200 is calculated using a rearranged Equation (27):
200 z om
ln
u 200 = u *
k
26
(29)
where; u* is the friction velocity at the weather station. This calculation is also done on a calculator or spreadsheet. 3. The friction velocity (u*) for each pixel is computed using the ERDAS Model Maker tool. u200 is assumed to be constant for all pixels of the image since it is defined as occurring at a “blending height” unaffected by surface features. From Equation (27):
u* =
ku 200
(30)
200 ln z om
where; zom is the particular momentum roughness length for each pixel. Momentum Roughness. The momentum roughness length (z om) for each pixel can be computed by two methods: 1) Using a land-use map: When a land-use map is available, values of zom for non-agricultural surface features (such as water, forest, desert grass, snow, manmade structures, etc.) can be assigned using the values given in Appendix 2. For agricultural areas, zom is calculated as a function of Leaf Area Index (LAI): zom =LAI 0.018 ×
(31a)
To compute zom, model F11_surface_roughness, is loaded and run. The file names for LAI and the land-use map are specified. Only the section of the land-use map that corresponds to the area of interest is used, so one should define the proper extent of the land-use map (the map has been previously created outside of SEBAL). The extent is specified by transferring the four coordinates defining the extent of the LAI image into the land-use map using the following steps: View the LAI image and record the four coordinates from the image info icon. In the model, click on the land-use map and enter it from a CD or disk. Enter the four coordinates of the LAI image. Save the output file and the model. 2) Using NDVI and surface albedo data: When a land-use map is not available, the second method for computing the momentum roughness length using NDVI and surface albedo data by Bastiaanssen (2000) and modified by Allen (2001)4 can be used. In this method, zom is computed from the following equation: 4
Bastiaanssen, W.G.M. 2000. SEBAL-based sensible and latent heat fluxes in the irrigated Gediz Basin, Turkey. J. Hydrology 229:87-100. Allen, R.G. 2001. Prediction of regional ET for the Tampa Bay, Florida watershed using SEBAL. personal communication. 27
zom = exp[(a × NDVI/α) + b]
(31b)
where; a and b are correlation constants derived from a plot of ln(z om) vs NDVI/α for two or more sample pixels representing specific vegetation types. Equation (31b) is an empirical expression that must be fitted to the local vegetation and conditions. The use of albedo (α) to modify NDVI helps to distinguish between some tall and short vegetation types that may have similar NDVI. For example, trees generally have smaller albedo than agricultural crops. Other variations on Equation (31b) can be explored by the user to improve the predictive accuracy. To determine a and in Equation a series of sample representing vegetation types andbconditions of(31b), interest are selected frompixels the respective image and the associated values for NDVI and surface albedo are obtained. For each type or condition, a characteristic vegetation height (h) is assigned. A value forom z is computed using Equation (28), however the constant (0.12) may vary depending on the type of vegetation represented by the pixel (for agricultural crops, 0.12 is used; for an open forest, 0.2 may be more appropriate; for a tight canopy forest, a value less than 0.2 may be more appropriate). The ln(z om) is plotted against NDVI/α in a spreadsheet and the correlation coefficients, a and b, are derived by regression or visual fitting. Zom for each pixel is now computed using model F11_surface_roughness_using_NDVI that utilizes a fitted Equation (31b). The files for surface albedo and NDVI are input along with values for the correlation coefficients a and b. Save the output file and the model.
28
Weather station H for each pixel
u, zx, zom , u*
wind speed at 200 meters
friction velocity at each pixel
u 200
= u*
u* =
z 200 ln z 0m
k
ψm (blending) ψh (z1) ψh (z2)
z2 z1
ln
u*
u* × k
Cold Pixel
Hot Pixel
Hcold=Rn-G-λETcold
Hhot=Rn-G
dTcold=Hcold*rah/(ρCp)
dThot = Hhot*rah /(ρCp)
dT for eac h pixel
rah 3
z ln 200 z0m
=
dT
* s L = − ρ Cp u T kgH
ku200
rah for each pixel
rah
H = ρCp
=
k u 200
z ln 200 − Ψ m(200) z 0m
z2 − Ψ h (z ) + Ψ h (z )
ln
rah
z = 1
2
1
u* × k
dT=aTs+b
Figure 4. Flow Chart of the Iterative Process for the Calculation of Sensible Heat (H)
29
4. The aerodynamic resistance to heat transport (r ah) is computed in model F12 as is u * in step 3. A series of iterations is required to determine the value forahr for each period that considers the impacts of instability (i.e., buoyancy) on ah r and H. Assuming neutral atmospheric conditions, an initial ah r is computed using Equation (26). z1 is the height just above the zero plane displacement (d≅ 0.67 × height of vegetation) for the surface or crop canopy and z2 is some distance above the zero plane displacement, but below the height of the surface boundary layer. Based on experienced analysis, values of 0.1 meter for z1 and 2.0 meters for z2 are assigned. Steps 3 and 4 are completed by running model
F12 to compute initial values of u* and rah. The file for zom from model F11 is input along with u200, which was calculated above with equation (29). Save the output files and the model5. 5. To compute the sensible heat flux (H) from Equation (25), the near surface temperature difference (dT) for each pixel needs to be defined. This is given as dT = Tz1 – Tz2. The air temperature at each pixel is unknown, along with explicit values for Tz1 and Tz2. Therefore, only the difference dTis utilized. SEBAL computesdT for each pixel by assuming a linear relationship between dT and T s:
dT = b + aTs
(32)
where; b and a are the correlation coefficients. To define these coefficients, SEBAL uses the two “anchor” pixels where a value for H can be reliably estimated. Values for H and dT at these “anchor” pixels are computed in a spreadsheet format as explained in Appendix 8 and in the following steps. The linearity of the dT vs. Ts function is a major presumption of SEBAL. However, research by Bastiaanssen and others andby the University of Idaho at Kimberly, indicate that this assumption appears to fit a large range of conditions. a) At the “cold” pixel, we define the sensible heat flux as Hcold = Rn – G - λETcold from Equation (23). Experience in Idaho shows that the most “cold” (wet) agricultural fields have an ET rate about 5% larger than the reference ET (ET r). Hence, ETcold is assumed to be 1.05 × ETr (which was calculated earlier in the SEBAL process, Appendix 3). Hcold is now calculated in the spreadsheet as: Hcold = Rn – G – 1.05λETr. dTcold is computed in the spreadsheet using the inverse of Equation (25):
dTcold = Hcold × rah_cold / (ρcold × cp)
5
(33)
Note that because dT is defined as ∆T between z1 and z2, the estimate of Tair as Tair = Ts + dT is not a
real, representation of air temperature. This is because the T at z1 is generally somewhat less than Ts. In addition, Ts is a radiometrically derived surface temperature, whereas T at z1 is an aerodynamically representative temperature.
30
Note: An alternative method used by Bastiaanssen is to choose the “cold” pixel at a body of water where one assumes Hcold = 0; and therefore, λETcold = Rn – G and dTcold = 0. The use of Equation (33) and ETr helps to consider the impact of regional advection and wind speed on ET. It also helps to account for any error in image estimates for albedo and R n by SEBAL. b) At the “hot” pixel, Hhot = Rn – G - λEThot; where EThot is assumed to be zero for a “hot” (dry) agricultural field having no green vegetation and with dry soil surface layer. The weather data should be checked to see if this assumption is correct. If there was some precipitation 1-4 days before the image date, then ET hot should
be estimated using a water balance model and tracking the soil moisture of the “hot” pixel. Hhot is now calculated in the spreadsheet and dThot computed from Equation (33).
Plotting dTcold vs. Ts_cold and dThot vs. Ts_hot gives the following graph:
The correlation coefficients, b and a, are now computed for the linear relationship: dT = b + aTs. The temperature differential (dT) for each pixel can now be computed using the coefficients b and a and the surface temperature (T s). The determination of dT at the anchor pixels and b and a are done in the spreadsheet. Figure 4 shows a flow chart of the process. 6. An approximation for air temperature (Ta) for each pixel is computed as Ta = Ts – dT and an approximation for air density ρ ( ) is calculated in the spreadsheet for the anchor
pixels. 7. The sensible heat flux (H) for each pixel is computed in SEBAL using Equation (25). This is the first estimation of H assuming neutral atmospheric conditions. Model F13 is
31
used for this computation. The file for the initial value ofahr from model F12, the surface temperature file, and the correlation coefficients b and a are input. The output is saved as H1 and the model is saved. 8. In order to account for the buoyancy effects generated by surface heating, SEBAL applies the Monin-Obukhov theory in an iterative process (Figure 4). Atmospheric conditions of stability have a large effect on the aerodynamic resistance ah (r) and must be considered in the computation of sensible heat flux (H), especially for dry conditions. Appendix 11 gives details of different stability conditions that can occur in the atmosphere. SEBAL repeats the computation of H through a number of iterations, each
one correcting for buoyancy effects, until the value forahr stabilizes. The Monin-Obukhov length (L) is used to define the stability conditions of the atmosphere in the iterative process (note, this is not the same “L” as used in the SAVI computation). It is a function of the heat and momentum fluxes and is computed as follows: 3
L=−
ρc p u* Ts kgH
(34)
where; ρ is the density of air (kg/m2), cp is the air specific heat (1004 J/kg/K), u* is the friction velocity (m/s), Ts is the surface temperature (K), g is the gravitational constant (9.81 m/s2), and H is the sensible heat flux (W/m2). This computation is carried out in model F14. Values of L define the stability conditions of the atmosphere. If L<0, the atmosphere is considered unstable; if L>0, the atmosphere is considered stable; if L=0, the atmosphere is considered neutral (Appendix 11). Depending on the atmospheric conditions, the values of the stability corrections for momentum and heat transport (ψm and ψh) are computed by the model F14 as follows, 6 using the formulations by Paulson (1970) and Webb (1970) : If L<0; unstable conditions:
1 + x ( 200 m ) 2 1 + x ( 200 m ) − 2 ARCTAN (x ( 200 m ) ) + 0.5π + ln 2 2
Ψm ( 200 m ) = 2 ln
2 1 + x ( 2m ) 2
Ψh ( 2m ) = 2 ln
(35)
(36a)
6
Paulson, C.A. 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. . Appl. Meteorol. 9:857-861.
Webb, E.K. 1970. Profile relationships: the log-linear range, and extension to strong stability. Quart. J. Roy. Meteorol. Soc. 96:67-90. 32
2 1 + x ( 0.1m ) 2
Ψh (0.1m ) = 2 ln
(36b)
where;
x (200m )
= 1 − 16
x (2m ) x (0.1m )
= 1 − 16
200
2 = 1 − 16 L 0.1
0.25
L
(37a)
0.25
(37b)
0.25
L
(37c)
If L>0; stable conditions:
2
ψ m ( 200 m ) = −5 L
2 L 0 .1 ψ h (0.1m ) = −5 L ψ h (2 m ) = −5
(38)
(39a) (39b)
If L=0; neutral conditions:ψm and ψh = 0 9. A corrected value for the friction velocity (u *) is now computed for each successive iteration as:
u* =
u 200 k
200 − Ψm ( 200 m ) ln z 0m
(40)
where; u200 is the wind speed at 200 meters (m/s), k is von Karman’s constant (0.41), zom is the roughness length for each pixel (m), andψm(200m) is the stability correction for momentum transport at 200 meters (equation 35 or 38). This computation is carried out in model F14. 10. A corrected value for the aerodynamic resistance to heat transport ah (r) is now computed during each iteration as:
33
rah
=
z ln 2 − Ψh ( z ) 2 z1
+ Ψh ( z1 ) (41)
u* × k
where; z2 = 2.0 meters, z1 =0.1 meters, and ψh(z2) and ψh(z1) are the stability corrections for heat transport at 2 meters and 1 meters (equations 36 or 39). This computation is carried out in model F14. 11. Return to step 5 and the spreadsheet to compute new dT values for the “cold” and ah. New values for the correlation coefficients b and a are “hot” pixel usingand thethen corrected also computed dT forreach pixel is revised as dT = b + aT s.
12. Step 6 is repeated for a revised value for air temperature (T a) and air density (ρ). 13. Step 7 is repeated to compute a corrected value for H. 14. Step 8 is repeated to compute a new stability correction. 15. This iterative process is repeated until the successive values for dT hot and rah at the “hot” pixel have stabilized.
Refer to Figure 4 for a flow chart of the iteration process and to Appendix 8 for instruction on the use of the spreadsheet. 16. Next, compute the final corrected value for the sensible heat flux (H) at each pixel,
which will be used in the computation of the instantaneous ET at each pixel. C. Latent Heat Flux ( ET), Instantaneous ET (ETinst), and Reference ET Fraction (ETrF) Latent heat flux is the rate of latent heat loss from the surface due to evapotranspiration. It can be computed for each pixel using Equation (1):
λET = –RnG – H
(1)
2 where; λET is an instantaneous value for the time of the satellite overpass (W/m ).
1. An instantaneous value of ET in equivalent evaporation depth is computed as:
ETinst = 3600 λET λ
34
(42)
where; ETinst is the instantaneous ET (mm/hr), 3600 is the time conversion from seconds to hours, and λ is the latent heat of vaporization or the heat absorbed when a kilogram of water evaporates (J/kg). 2. The Reference ET Fraction (ETrF) is defined as the ratio of the computed instantaneous ET (ETinst) for each pixel to the reference ET (ETr) computed from weather data:
ETrF =
ETinst
(43)
ETr
where; ETinst is from Equation 41 (mm/hr) and ETr is the reference ET at the time of the image from the REF-ET software (mm/hr). ETrF is similar to the well-known crop coefficient, Kc. ETrF is used to extrapolate ET from the image time to24-hour or longer periods. One should generally expect ETrF values to range from 0 to 1. At a totally dry pixel, ET = 0 and ETrF = 0. A pixel in a well established field of alfalfa or corn can occasionally have an ET slightly greater than ETr and therefore ETrF > 1, perhaps up to 1.1. However, ETr generally represents an upper bound on ET for large expanses of wellwatered vegetation. Negative values for ETrF can occur in SEBAL due to systematic errors caused by various assumptions made earlier in the energy balance process. ETinst and ETrF are computed in model F25 along with 24-hour ET as shown in Section 5 below.
5. 24-Hour Evapotranspiration (ET24) Daily values of ET (ET24) are often more useful than instantaneous ET. SEBAL computes the ET24 by assuming that the instantaneous ETrF computed in model F25 is the same as the 24-hour average. Figure 5 shows a plot of ET and ET rF vs. time. This graph shows how ET varies throughout the day while ETrF is relatively constant. Finally, the ET24 (mm/day) can be computed as:
ET24 = ETrF × ETr _ 24
(44)
where; ETr-24 is the cumulative 24-hour ETr for the day of the image. This is calculated by adding the hourly ETr values over the day of the image. ET24 is computed in model F25. Enter the files for H, G, R n, Ts, and the values for ETr and ETr-24. Save the output files and the model.
35
4/18/1989, Sugar Beet Data by Dr. Jim Wr ight, Kimberly ID
1
1
)r 0.8 h / 0.6 m m ( 0.4 T E 0.2
0.8 0.6 0.4
F r T E
0.2
0
0 0 0 : 7
0 0 : 9
0 0 : 1 1
ETSB(mm) ETrF
0 0 : 3 1
0 0 : 5 1
0 0 : 7 1
ETr(mm) ETrF24
0 0 : 9 1
ETo(mm)
5/20/1989, Sugar Beet Data by Dr. Jim Wright, Kimberly ID
1
1
)r 0.8 h / m0.6 m ( 0.4 T E 0.2
0.8 0.6 0.4
F r T E
0.2
0
0 0 0 : 7
ETSB(mm) ETrF
0 0 : 9
0 0 : 1 1
0 0 : 3 1
0 0 : 5 1
ETr(mm) ETrF24
0 0 : 7 1
0 0 : 9 1
ETo(mm)
Figure 5. Plot of ET and ETrF vs. time (from observations of Dr. Jim Wright, Kimberly, Idaho
36
Figure 6 shows an example of a daily evapotranspiration map. The range of values for ET 24 can be colored using the ERDAS Imagine software so that ET intensity can be readily observed on the image.
Path 39: Am. Falls - 24-hour ET 8/07/00
Figure 6. A Map of 24-hour Evapotranspiration
37
6. Seasonal Evapotranspiration (ET seasonal) A seasonal evapotranspiration map that covers an entire growing season is often valuable. This can be derived from the 24-hour evapotranspiration data by extrapolating the ET 24 proportionally to the reference evapotranspiration (ET r). We assume that the ET for the entire area of interest changes in proportion to the change in the ET r at the weather station. ETr is computed for a specific location and therefore does not represent the actual condition at each pixel. This does not matter, however, since ET r is used only as an index of the relative change in weather,and therefore ET, for the imagearea. We also assume that the ETrF computed for the time of the image is constant for the entire period represented by the image. The following steps show the process for computing seasonal ET: 1. The first step is to decide the length of the season for which ET is desired (i.e., March 1 to October 31). 2. The second step is to determine the period represented byeach satellite image within the chosen season (i.e., if the first two images for the above season are on th March 15 and April 8, then the period represented by the March 15 image would be March 1 to March 27). 3. The third step is to compute the cumulative ETr for the period represented by the image. This is simply the sum of daily ETr values over the period. These 24-hour values can be computed using the University of Idaho REF-ET software described in Appendix 3. The same ETr method must be used through the SEBAL process. ETr should represent the alfalfa reference ET, which is a larger value than for the clipped grass reference (ET o). 4. The fourth step is to compute the cumulative ET for each period as follows:
ET period = ETr F period
n
∑ ET −
r 24
(45)
1
where; ETrFperiod is the representative ETrF for the period, ETr-24 is the daily ETr, and n is the number of days in the period. Units for ETperiod will be in mm when ETr-24 is in mm/day. 5. The fifth step is to compute the seasonal ET by summing all of the ETperiod values for the length of the season. The following difficulties encountered in the computation of seasonal ET must be understood:
38
1. If there is some cloud cover in one of the images used for the seasonal ET computation, then there can be no ET or ETrF values represented for this area during the period represented by the image. One can assign ET rF values to the cloud-covered areas by interpolating between the ET rF values for the images on either side of the cloud-covered image. 2. When the ETr at the weather station is not representative of the entire area (i.e., where the image spans many mountain valleys or includes both coastal and inland areas having different effects from clouds or wind that are not well correlated across the image), then the image should be broken down into subimages, each with ownsub-image. value of ETr. In this situation, simple SEBAL runs should be made foritseach 3. If an agricultural field is dry onthe day of the image and is then irrigated on a following day, the increase in ET for this field during the representative period of the image will not be considered in the computation of seasonal ET. This problem can be minimized if many satellite images are used for the computation of seasonal ET with each image representing a shorter time period. One can also apply a water balance model to each pixel that includes irrigation and precipitation as inputs. Figure 7 shows an example of a seasonal evapotranspiration map that can be created for the area of interest:
Seasonal ET for Bear River Basin
Figure 7. A Map of Seasonal Evapotranspiration
39
Appendix 1 Acquiring Header File Information
The header file of the satellite image contains important information for SEBAL processing. The following information must be obtained from the header file: • The satellite overpass date and time • The latitude and longitude of the image • The sun elevation angle,β, at the overpass time (used in the “Flat” Model, only) • Gain and Bias levels for each band (applicable for Landsat 7 ETM+, only) There are three methods for obtaining the header file information depending on the satellite type (Landsat 5TM or Landsat 7ETM+) and on the format used to store the image on the distribution disk. Method A: Method A is the most straightforward method for acquiring header file information. It is applicable for the following satellite types and data formats: • Landsat 5 TM, where the srcinal image format is NLAPS • Landsat 7 ETM+, where the srcinal image format is NLAPS or FAST
For the NLAPS format, explore the srcinal image disk. Find and open the small file titled “NLAPS Correction Processing Report.” This file contains all the information needed including the gains and biases for Landsat 7. For band 6 of Landsat 7, the top row of the gain and bias section is for the low gain image and the second row is for the high gain image. We recommend using the low gain image as explained in section 2 B of the main text. For Landsat 7 – FAST format, explore the srcinal image disk. Open the bottom two small FST files. These contain the gains and biases for bands 1-7 with biases in the first column and gains in the second column. The list with six entries is for bands 1-5 and 7. The list with two entries is for band 6 (low gain on top, high gain below). The satellite overpass time is not included in the FAST header files and Method B must be used to acquire exact image time from the internet. Print the header file information for future reference. Method B: Method B is used to obtain header file information via the internet. It is applicable for the following data formats and satellite types: • Original image format is NLAPS and satellite type is Landsat 5TM (if the header file is missing) • Original image format isFAST and satellite type is Landsat5TM (regardless the header file availability) • Original image format isFAST and satellite type is Landsat7ETM+ (only if header file is available)
40
• For all other formats (including GeoTIFF) where the satellite type is Landsat 5 TM (regardless the header file availability) (1) First, go to the Earth Explorer web page:http://edcsns17.cr.usgs.gov/EarthExplorer/ (2) Click “Enter as a guest”. (3) In the “Data Set Selection” box, select TM (Landsat 5) or ETM+ (Landsat 7) depending on the image, and click “Continue”.
(4) Define “Acquisition Date” and “Path Row Search”, and click the “search” button.
(5) When the search is finished, select the appropriate image and click “Result”.
41
(6) Click “Show” in the “Show All Fields” category, for the appropriate image.
(7) All necessary header file information including the satellite overpass time (use either start time or stop time) is presented. It is recommended that the header file information be printed for future reference. Method C: Method C is applicable for the following data formats and satellite types: • Original image format is NLAPS and satellite type is Landsat 7ETM+ (if the header file is missing) • Original image format isFAST and satellite type is Landsat7ETM+ (if the header file is missing)
all other formats (including GeoTIFF) when satellite type is Landsat 7ETM+ • For (regardless the header file availability) The following is the procedure for obtaining header file information (including the satellite overpass time) from the internet for the specified data types. This procedure was obtained from Dr. Satya Kalluri of Raytheon Company: A. The core metadata that contains much needed info like solar zenith and azimuth angles, gain settings etc. will not be in the header of the GeoTIFF files! You should get this as a separate ascii file when you get your L7 data. If you do not have this, you can get it online for each scene using the ECS data gateway. This is how you go about doing it: Go to http://edcimswww.cr.usgs.gov/pub/imswelcome/ Login as a guest or register Under Choose Search Key Words - Method 2 - click on DATA SET and click on SELECT On the next screen, select LANDSAT-7 LEVEL-0R WRS-SCENE V002 under Data Set list 2, and click on OK. This will take you back to the main page.
42
On the main page, select Type in Path/Row Range under Choose Search Area Enter your path and row Enter your start date and end date Click on Start Search! The search should return your granule (image) as a result. Click on "Granule attributes" button, and this will give you all the metadata need including the gain settings and solar angles for that particular you scene. When you follow these instructions, you obtain the equivalent of the header file information as shown below. Landsat7 ETM+ P25/R39
9/20/1999
All of the data in the following table are available in the "Gra nule attributes" (Header of the imag e has equiverent information too) DateandTime,inGMT Gaincombination Centeroftheimage SunElevation Band12345(66)7 Lat Lon Degree Start Date 20 Sep 1999, 16:43:32 HHHHH(HL)H 30.31° Lat -94.79° Lon 53.9359 Stop Date 20 Sep 1999, 16:43:32 Band 6 has both High and Low gains I recommend to use Low gain (but High gain is okay too) Comment (Lowgainhasbitlessresolutionbutsafer) We do not use band 8 then I do not includes band 8 (band 8 is always Low gain)
Note that Band 6 has both High and Low gain images (files). For SEBAL, we recommend using the Low gain image, even though it has slightly less resolution, since it is less likely to suffer from saturation. SEBAL does not use Band 8 (the Panchromatic Band), so header information for this band is not described here.
43
Appendix 2 Land-Use Map and Momentum Roughness Length A. Land-Use Map A land-use map is used in SEBAL for determining the momentum roughness length (zom) and for establishing G/Rn values for various land-use types. If a land-use map is not available for the area of interest, the user can develop one using satellite data. The landuse map helps the user develop a good working knowledge of the area being studied so he can better recognize the reliability of the SEBAL computations. The land-use map should × 30m). be in a raster format with the same pixel size as the satellite image (30m
The land-use map for the Snake River Plain of southern Idaho contains the following landuse types along with a given pixel value: Pixel value 0 1 2 3 4 5 6 7 8 9 10 11 12 13 21 22
Land-use type background water city and manmade structures vegetated field forest (relatively flat) grassland sage brush bare soil burned area salty soil basalt (dark gray) basalt (black) basalt (gray) basalt (light gray) forest (mountain) bare soil (mountain), dead grass and other small vegetation
B. Momentum Roughness Length The land-use classification is used to assign values of om z to specific land types. The following are examples of zom estimations given in a land-use map:
a) b) c) d)
Water; Cities; Forests; Grassland;
e) Desert with vegetation; f) Snow;
zom zom zom zom
= = = =
0.0005 m 0.2 m 0.5 m 0.02 m
zom = 0.1 m zom = 0.005 m
44
These values for zom can vary by region and their use should be carefully considered. For agricultural areas, zom is computed as a function of crop height (h) which in turn is defined as a function of Leaf Area Index (LAI). For crops such as alfalfa, potatoes, beans, beets, and peas, h = 0.15LAI. Other empirical relationships between h and LAI for specific crops are given below (based on the work of M. Tasumi, Univ. Idaho, using data from Dr. J.L. Wright of the USDA-ARS, Kimberly, Idaho):
2
Alfalfa: H = 0.009LAI + 0.076LAI, 3 2 Corn: H = 0.03LAI - 0.2194LAI + 0.7243LAI, 3 2 Potato and Tall beans: H = 0.0043LAI - 0.0694LAI + 0.3783LAI, 3 2 Beans, Beet and Peas: H = 0.0025LAI - 0.0417LAI + 0.2754LAI, 3 2 Spring Wheat: H = 0.0414LAI - 0.1848LAI + 0.3214LAI, 3 2 Winter Wheat: H = 0.0039LAI - 0.0411LAI + 0.1747LAI,
R2 = 0.9881 R2 = 0.9938 R2 = 0.9837 R2 = 0.9756 R2 = 0.9896 R2 = 0.9914
The data from the figures above are compared against the h = 0.15 LAI function in the following figure (from Univ. Idaho, 2002):
45
LAI vs Crop Height, Kimberly ID (Data by Dr. J.Wright USDA/ARS)
2 1.8
1971 Alfalfa (3seasons)
1.6
1972 Potato
1.4
) m ( 1.2 t h 1 g i e 0.8 H 0.6
1973 Beans 1974 Beans 1975 Beets 1976 Corn 1977 Pias
0.4 0.2
1978 Wheat
0
1979 Wheat
02468
1:0.15 Line
LAI Note that in an area dominated by corn production, a different height-LAI function should be developed.
46
Appendix 3 Reference Evapotranspiration, ETr, Calculation
The reference evaotranspiration (ETr) is the ET rate expected from a well-defined surface of full-cover alfalfa or clipped grass. ETr is used in SEBAL to estimate the ET at the “cold” pixel and to calculate the reference ET fraction (ET rF). ETr can be computed for a given weather station using the REF-ET software made available by the University of Idaho: Allen, R.G. (2000). REF-ET: Reference Evapotranspiration Calculation Software for FAO and ASCE Standardized Equations , University of Idaho www.kimberly.uidaho.edu/ref-et/ This software is included with this SEBAL manual and contains the ET r program and instructions. A. Requirements for running REF-ET: Weather data for hourly or shorter time periods surrounding the date and time of satellite image is required. The weather data should include air temperature, solar radiation, wind speed, and humidity (in the form of dewpoint temperature, relative humidity, or vapor pressure). Data for other days can be processed at the same time and will be useful for predicting ET for adjacent time periods in SEBAL.
The data file should be in “text” format (i.e., ASCII), and data columns can be separated by tabs, spaces, or commas. The data columns can have headings and can be in any order and can use a wide range of data units. The order of the weather parameters and their units are defined within REF-ET. The weather data file should be created before running REF-ET. The following figure is an example of hourly weather data with headings for use in REFET. Weather station information is defined in section B-5 of this appendix.
47
B. Steps for running REF-ET 1. Open REF-ET and click on “Proceed”
2. Use the file browser as shown in the above box, identify the name of the data file,
48
and click on “Open” 3. Provide the name for a new REF-ET definition file (i.e., ref_et.def) tobe created by REF-ET. The definition file defines the order of weather data in the data file and the associated measurement units. If a definition file has been previously created for the weather data set, then this definition file can be reused by REF-ET as shown in the following dialogue box.
49
4. Indicate the data file parameter order andunits by double clicking on the list of data parameters and measurement units in the box in the upper right hand of the following REF-ET screen. Selected parameters will appear in the list in the upper left hand side of the following screen (green area).
Click on “Continue”
50
5. Describe weather station and weather data file
Enter “1” (or other value) for “Initial Lines to be skipped” if the first line(s) of the weather data is/are headings. Click on “Continue”
51
6. Enter output modes and reference ET equations. The followingselections are recommended:
Select SI units Select file and screen output to create a disk file of results Select the “ETr - ASCE Penman-Monteith Standardized Form” equation. This method is recommended as it represents a tall alfalfa surface that best represents the ET from highly vegetated pixels. The “intermediate” files contain an echo of the data file and various parameters involved in calculating reference ET. Click on “Continue” 7. Save the definition file for future use
52
8. Save the output file using a name selected by the user. Results from REF-ET that are echoed on the screen should look similar to the following screen:
Click on “End”
53
Appendix 4 Preparing the Original Satellite Image for Use in SEBAL
The srcinal Landsat image needs to begeo-rectified for use in SEBAL. This can be done locally or by sending it to a commercial entity such as the Earth Satellite Corporation. It is recommended that the rectified image by saved in the GeoTIF format and that nearest neighbor resampling be used during rectification to preserve spectral information. SEBAL uses the following seven bands of the spectrum: Band 1 Band 2 Band 3 Band 4 Band 5 Band 6 Band 7
visible (blue) visible (green) visible (red) near infrared near infrared thermal near infrared
These bands must be layered inside ERDAS Imagine, in order, from 1 to 7 to create an image file for use in the SEBAL process. After that, a smaller subset image can be created for the area of interest. The steps for the layering a Landsat 7ETM+ image are as follows (Figure 4.1): 1. Insert the geo-rectified disk (GeoTIF format) 2. View the details The six large files with a “t” in the name are the bands 1 - 5 and 7, in order. The two large files with a “k” in the name are for band 6, low gain and high gain in that order (the low gain is recommended). The large file with a “p” in the name is for band 8 (PanChromatic) and is not used in SEBAL. 3. In ERDAS, go to the Interpreter tool – Utilities - Layer stack and follow these steps: a. The input file will come from thegeo-rectified CD in a geo-TIFFformat. b. Create an output file having a unique name. It is recommended that the name include the date, path, and row of the image. c. Add the seven layers in order beginning with band 1 (the first large “t” file). After adding band 5, add band 6 (the first large “k” file for low gain). Then add band 7 (the last large “t” file). 3. A layered image is now created.
For Landsat 5TM images, layering is similar to above except that spectral bands 1 – 7 are listed in order in the header file (Figure 4.2). To create a smaller subset image follow these steps: 1. In ERDAS, go to the Data Prep tool - subset image. 2. Input the image file created above. 3. Create an output file name.
54
4. Go to the Viewer tool and view the layered image. Go to Utility -Inquire Box and fit the box to the desired area. Click on ‘Apply’. 5. Return to Data Prep and click on ‘From Inquire Box’. 6. The subset image file is now created. The area of interest can now be viewed with the ERDAS Viewer tool. The image is shown in false color. ERDAS has many tools for working with the image which the user must be familiar with. Explaining all of these tools is outside the scope of this manual.
Figure 4.1 Layering the image for Landsat 7ETM+
Band 6 (low & high)
Bands 1-5,7
55
Figure 4.2 Layering the image for Landsat 5TM
Bands 1-7 in order
56
Appendix 5 Time Issues Around Weather Data and Reference Evapotranspiration
A. Converting Greenwich Mean Time (GMT) to Local Time
Satellite image times are given in GMT and must be converted to the local time of the area of interest. For SEBAL processing, use winter standard time (local time without applying any daylight savings time shift): timage (Local Time) = GMT + Correction
where timage (Local Time) is the image time using local time, and “Correction” is the time correction that depends on the time zone. This correction term is found on the internet. Table 5.1 gives samples of the correction term for some locations: Table 5.1. List of the correction term (From Date/Time Property of Microsoft Windows)
For example, if the weather station is in the Mountain Time zone, the correction term is –7 hours. Therefore, if the satellite image time is 17:25 GMT, the local time of the image is 17:25 – 7 = 10:25 AM. The correction can also be determined by dividing the longitude of the center of the time zone by 15o. For example, the Mountain Time zone center is -105o. Therefore, the Correction = -105/15 = -7 hours. The centers of Eastern, Central, Mountain and Pacific Time zones are –75O, -90O, -105O, and –120o respectfully.
57
B. Determining the Time for Weather Data and Reference Evapotranspiration (ET r)
This section describes how the instantaneous wind speed and the reference ET at the satellite overpass time are determined from hourly (or shorter) weather data. This procedure assumes that one has hourly (or shorter) meteorological data from a weather station. Refer to part 2D in “The Theoretical Basis of SEBAL” for the selection of weather stations and required weather data. B-1. Data Investigation
Any data from weather station must be studied determine what interval it represents. Thea table below shows weather datatoobtained from thetime following weather station: (6/20/2000, Aberdeen Idaho weather station in the Agrimet Public Weather System, http://mac1.pn.usbr.gov/agrimet/). ETr was computed using the REF-ET software (Appendix 3).
DoY HrMn Year
Tmax Tmin
172 172 172 172 172 172 172
0 100 200 300 400 500 600
0 0 0 0 0 0 0
C C 10.9 10.9 10.2 10.2 10.3 10.3 9.3 9.3 7.3 7.3 5.8 5.8 3.6 3.6
172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172
700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3.7 5.7 10.7 12.9 14.7 16.2 17.6 18.8 20.3 21.2 21.7 22.2 22.0 20.9 19.1 17.1 14.6
3.7 5.7 10.7 12.9 14.7 16.2 17.6 18.8 20.3 21.2 21.7 22.2 22.0 20.9 19.1 17.1 14.6
ASCE Wind DewP stPM ETr W/m2 m/s C mm/h 0 1.7 8.8 0.00 0 2.3 7.3 .01 0 4.1 3.6 .04 0 1.3 2.6 .01 0 .3 3.6 -.02 0 1.0 3.1 .02 1 1.2 2.0 .01 Rs
64 1.7 1.6 .05 156 1.9 2.6 .10 415 .4 2.8 .25 591 .7 1.2 .38 747 1.8 .3 .54 868 3.4 .9 .68 941 4.5 .5 .79 964 6.2 1.5 .87 939 4.5 .3 .87 869 4.8 -.3 .87 619 5.4 .2 .76 469 5.6 1.6 .69 402 5.7 .7 .66 226 3.9 1.4 .41 69 3.7 2.5 .12 3 3.8 .8 .10 1 3.3 2.0 .06
When using these data, one must ask two questions: 1. Was daylight savings time, DST, applied to the time of record?
58
2. What time interval does the data represent? (i.e., doesthe posted time representthe end of the measurement period, the center of the measurement period, or the beginning of the measurement period?) In the example above, DST was applied and the data represent averaged values over the hour ending at the reported time. For example, 11:00 AM data represents the average of values from 10:00 AM – 11:00 AM. B-2. Instantaneous Wind Speed and ETr Calculation This section describes how to determine the instantaneous wind speed at the image time using the weather data. Reference ET (ETr) can be determined by the same method.
(1) Calculate the satellite overpass time in local standard time (i.e., without applying the daylight saving time): timage (Local Time) = GMT + Correction
(2) Determine the weather time periods to use for interpolation. Quite likely, the weather parameter for the instantaneous image time will need to be interpolated from averages representing two time periods. A general formula to determine which time periods to use in the interpolation is the following:
t image (local time)
t1
= int
t2
= t 1 + ∆t
∆t
1
2
+ − Flag period ∆t + Flag DST
where; t1 = time identifier for the first weather period in hours t2 = time identifier for the second weather period in hours ∆t = length of weather period in units consistent with t1 and t2 (i.e., hours) timage (local time) = the satellite image time (local time, with no adjustment for daylight savings) Flagperiod = a flag for the way that the weather data period is represented by its“time label”. Flagperiod = 0 if the time label is the endpoint of the weather data period; = 1 if the time label is the beginning of the weather data period; = 0.5 if the time label represents the midpoint of the weather data period. FlagDST = a flag for use of DST in the weather data set. FlagDST = 0 if no DST is used and FlagDST = 1 hour if DST is used in the weather data set. (3) Estimate the instantaneous wind speed applying a simple linear interpolation assumption.
59
For interpolation of data we assume that the averages for the periods represent measurements occurring at the midpoints of the periods. A general interpolation formula based on the t1 and t2 period identifiers from step (2) is: Data image
= Data t 1 +
Data t
2
− Data t 1 ∆t t image ( Local time) − t1 − Flag DST + ∆t Flag period − 2 ∆t
where; Datat1
= the average value for the specific data parameter (i.e., wind or ETr) for the t1 time period
Datat2
= the average value for the specific data parameter (i.e., wind or ETr) for the t2 time period
Example: Using the weather data from section B-1. The Image Overpass Date&Time: 6/20/2000, 17:49 GMT Location of weather station: In Mountain Time zone (time correction is -7) Weather Data: Given in section B-1 (Aberdeen, Idaho; based on DST) timage (Local Time) = 17:49 – 7:00 = 10:49 am
For this weather data DST was applied and the time identifier (label) for the weather data represented the endpoint of the data period (Flag DST = 1, Flagperiod = 0). ∆t = 1 hour. Therefore, t1 and t2 are calculated as:
(10 + 49 / 60) 1 + − 0 (1) + 1 = int[11.3](1) + 1 = 12 hours 1 2
t1
= int
t2
= 12 + 1 = 13 hours
Therefore, the two weather data periods used to determine the wind speed and ET r at the image time are 1200 and 1300 hours, Mountain Standard Time. Calculate the instantaneous wind speed at 10:49 AM (Mountain Standard Time) by applying a linear interpolation using the reported time periods of 1200 and 1300 hours (Mountain Daylight Savings Time, MDT).
60
The wind speed is 3.4 m/s for the 1200 period and 4.5 m/s for the 1300 period (see weather data). The wind speed at 10:49 AM (i.e., image time) is calculated as:
Data image
= 3.4 +
(4.5 − 3.4) 1
1 (10 + 49 / 60) − 12 − 1 + 1(0) − 2 = 3.75 m / s
61
Appendix 6 Tables for Surface Albedo Computations
Table 6.1. LMIN and LMAX values for Landsat 5 TM (Left: for 1984 from Markham and Barker, 1986, Right: for 2000 calibrated by Univ. Idaho 2002)
After 15, Jan, 1984 -2
-1
Band W*m *ster Number LMIN 1 -1.500 2 -2.800 3 -1.200 4 -1.500 5 -0.370 6 1.238 7 -0.150
For 2000
-1
-2
* m LMAX 152.100 296.800 204.300 206.200 27.190 15.600 14.380
-1
Band W*m *ster Number LMIN 1 -1.765 2 -3.576 3 -1.502 4 -1.763 5 -0.411 6 1.238 7 -0.137
-1
* m LMAX 178.941 379.055 255.695 242.303 30.178 15.600 13.156
Table 6.2. LMAX and LMIN values for Landsat 7 ETM+ (Landsat 7 Science User Data Handbook Chap.11, 2002)
W*m-2*ster-1*µm-1 B e f o reJ u l y1 ,2 0 0 0 Band Number 1 2 3 4 5 6 7 8
L o wG a i n LMIN LMAX -6.20 297.50 -6.00 303.40 -4.50 235.50 -4.50 235.00 -1.00 47.70 0.00 17.04 -0.35 16.60 -5.00 244.00
H i g hG a i n LMIN LMAX -6.20 194.30 -6.00 202.40 -4.50 158.60 -4.50 157.50 -1.00 31.76 3.20 12.65 -0.35 10.932 -5.00 158.40
Af t e rJu l y1 ,2 0 0 0 L o wG a i n LMIN LMAX -6.20 293.70 -6.40 300.90 -5.00 234.40 -5.10 241.10 -1.00 47.57 0.00 17.04 -0.35 16.54 -4.70 243.10
H i g hG a i n LMIN LMAX -6.20 191.60 -6.40 196.50 -5.00 152.90 -5.10 157.40 -1.00 31.06 3.20 12.65 -0.35 10.80 -4.70 158.30
Band 8 (panchromatic band) is not used in SEBAL. The gain combination for the particular image must be clear when selecting the LMIN and LMAX values in Table 2. These calibration constants are best used for 1999to 2001 images. For the future use, it is recommended to access the “Landsat 7 Science User Data Handbook” to see if the latest calibration constants are available.
62
Table 6.3. ESUNλ for Landsat 5 TM (Markham and Barker, 1986), and for Landsat 7 ETM+ 2 (Landsat 7 Science User Data Handbook Chap.11, 2002), both are in W/m /µm
Landsat5 Landsat7
Band1 1957 1969
Band2 1829 1840
Band3 1557 1551
Band4 1047 1044
Band5 219.3 225.7
Band6 -
Band7 74.52 82.07
Band6 -
Band7 0.011 0.012
Note: a dummy value of 1 is entered for band 6.
Table 6.4. Weighting coefficients, Landsat5 Landsat7
Band1 0.293 0.293
Band2 0.274 0.274
Band3 0.233 0.231
Band4 0.157 0.156
Band5 0.033 0.034
Note: a dummy value of 0 is entered for band 6. Table 6.5. Number of the day of the year. Number of the day in the year (J) Day
January
1 2 3 4 5 6 7
1 2 3 4 5 6 7
32 33 34 35 36 37 38
60 61 62 63 64 65 66
91 92 93 94 95 96 97
121 122 123 124 125 126 127
152 153 154 155 156 157 158
8 9 10
8 9 10
39 40 41
67 68 69
98 99 100
128 129 130
159 160 161
11 12 13 14 15 16 17 18 19 20
11 12 13 14 15 16 17 18 19 20
42 43 44 45 46 47 48 49 50 51
70 71 72 73 74 75 76 77 78 79
101 102 103 104 105 106 107 108 109 110
131 132 133 134 135 136 137 138 139 140
162 163 164 165 166 167 168 169 170 171
52 53 54 55 56 57 58 59 (60) -
80 81 82 83 84 85 86 87 88 89 90
111 112 113 114 115 116 117 118 119 120 -
141 142 143 144 145 146 147 148 149 150 151
172 173 174 175 176 177 178 179 180 181 -
21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 *add 1 if leap year
February
March*
63
April*
May*
June*
TABLE 6.5 (continued) Number of the day in the year (J) Day
July*
August*
1 2 3 4 5 6 7 8 9 10
182 183 184 185 186 187 188 189 190 191
213 214 215 216 217 218 219 220 221 222
11 12 13 14 15 16 17 18 19 20
192 193 194 195 196 197 198 199 200 201
21 22 23 24 25 26 27
202 203 204 205 206 207 208
28 209 29 210 30 211 31 212 *add 1 if leap year
September*
October*
November*
December*
244 245 246 247 248 249 250 251 252 253
274 275 276 277 278 279 280 281 282 283
305 306 307 308 309 310 311 312 313 314
335 336 337 338 339 340 341 342 343 344
223 224 225 226 227 228 229 230 231 232
254 255 256 257 258 259 260 261 262 263
284 285 286 287 288 289 290 291 292 293
315 316 317 318 319 320 321 322 323 324
345 346 347 348 349 350 351 352 353 354
233 234 235 236 237 238 239
264 265 266 267 268 269 270
294 295 296 297 298 299 300
325 326 327 328 329 330 331
355 356 357 358 359 360 361
240 241 242 243
271 272 273 -
301 302 303 304
332 333 334 -
362 363 364 365
The number of the day in the year, J, between 1 (January 1) and 365 or 366 (December 31) can be calculated as follows:
J = DM − 32 + INT 275
M 3 M Mod (Y ,4 ) − + 0.975 + 2 INT + INT 9 4 M + 1 100
where; DM = the day of the month (1-31) M = the number of the month (1-12) Y = the number of the year (1996 or 96) The INT function finds the integer number of the argument in parenthesis by rounding downward. The Mod(Y,4) function finds the modulus (remainder) of the quotient Y/4.
64
Table 6.6. Typical albedo values.
Fresh snow Old snow and ice Black soil Clay White-yellow sand Gray-white sand Grass or pasture
0.80 – 0.85 0.30 – 0.70 0.08 – 0.14 0.16 – 0.23 0.34 – 0.40 0.18 – 0.23 0.15 – 0.25
Corn field 0.14 0.22 Rice field 0.17 – – 0.22 Coniferous forest 0.10 – 0.15 Deciduous forest 0.15 – 0.20 Water 0.025 – 0.348 (depending on solar elevation angle) (Hsrcuchi, Ikuo (Ed.) 1992. Agricultural Meteorology, Buneidou, Tokyo, Japan )
65
Appendix 7 Methods for Selecting “Cold” and “Hot” Pixels
The selection of the two “anchor” pixels requires knowledge and skill. Below is a description of the steps for selecting these pixels, followed by an illustrative example.
A. “Cold” Pixel Selection Philosophy of the Cold Pixel. The cold pixel is used in SEBAL to define the amount of ET (or more directly, the amount of sensible heat, H) occurring from the most well-watered and fully vegetated areas of the image. It is presumed that these areas represent instances where the maximum amount of available energy is being consumed by evaporation. The impact of the near maximum evaporation rates is to cool the surface below that of areas having less ET.
In the traditional application of SEBAL by Bastiaanssen, the cold pixel is generally selected from a water body and it is assumed that ET = Rn – G (i.e., that all availableenergy is used to evaporate water, so that H = 0 for the cold pixel. In the Idaho extension to SEBAL, we have elected to presume that ET at the cold pixel is closely predicted by the ET rate from a large expanse of alfalfa vegetation. We therefore assume that ET = 1.05 ET r at the cold pixel, where ETr is the rate of ET from the alfalfa reference, defined using the ASCE Standardized Penman-Monteith equation for alfalfa. The ET r is calculated using weather data measured at a weather station within or very close to the image. The use of 1.05 ET r to represent the ET from the coldareas pixeland helps incorporate effects advection heat energy from surrounding desert thetoimpact of wind speedofon ET. Theof 1.05 factor is applied because it is likely, in a large image, that some fields may have a wet soil surface beneath the vegetation canopy that might increase the total ET rate about 5% above that of the ETr standard. H for the cold pixel is then calculated as H = R n – G – 1.05 ETr. The cold pixel should therefore be selected to represent an agricultural area that has vegetation at “full cover” and that is well irrigated. Full cover means that the LAI is greater than about 3. Instructions The following steps necessary for the selection of the “cold” pixel: 1. Open the image or subset image and view it in true color (layer combination of3,2,1) and in false color (layer combination of 4,3,2). 2. Open the Ts image in a new viewer with the Pseudo Color option. Link the two images. 3. On the Ts image use “Raster – Attributes” to observe the range of T s values.
66
4.
Make a first guess for Tcold by looking at some fully covered agricultural fields (dark green). Insert a color gradation around the selected temperature.
5. View the colored Ts image and select a point that represents a “cold” (wet) agricultural field. Do not select an extreme cold point but one that is representative of a well-watered full cover crop. For the best results, the cold pixel should have a surface albedo in the range of 0.22 to 0.24 (corresponding to a full covered reference alfalfa field). It should have a Leaf Area Index in the range of 4 to 6 (corresponding to full covered agricultural field). Finally, the cold pixel should be located nearby the weather station so that the weather conditions are similar. 6. Observe other areas of the image for similar temperatures. 7. Select Tcold from a pixel in the center of the chosen field that represents a very cold, but not extremely cold, point in the image. 8. Record the coordinates and temperature of the “cold” pixel.
Illustrative Example (Landsat 7 ETM+, 4/8/2000)
The following example illustrates some of the thought processes employed during selection of a cold pixel. In the following example, many well-watered agricultural fields, during an initial viewing, are found to have surface temperatures of around 290 K. This can serve as an initial guess for the cold pixel temperature. Next, we create a graduated color map for Ts over a range of 284-295 K (Figure 1) that includes the 290 K temperature. In this picture, the blue colored pixels are near 284 K and the red pixels are near 295 K. The image in Fig. 1 is an early April image, and not many fields are well-watered for this period in this region. One can see that most of the agricultural fields are still bare fields and cultivation has not been started. In the coloreds T image, bluish fields are more likely to be well-watered fields having some vegetation cover.
67
Fig.1. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image(left, 284K:Blue to 295K:Red) and the true color image (right).The coordinates of the inquire cursor are x/y = 436892/194372)
First, let’s focus on the northwest region. There are many cold pixel candidate fields in this region. Some of the candidate fields are circled in Figure 2.
Fig.2. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (left, 284K:Blue to 295K:Red) and the false color image (right).The coordinates are x/y = 435805/188719)
One of the candidate fields is displayed (marked by the cross-hairs) at the center of the next image (Figure 3). The temperature is quite homogeneous in the field, around 289.5 K. This is a good candidate field because the false color image indicates a high intensity of vegetation. In addition, there are some fields in other areas having similar temperatures. Based on the local area (south-central Idaho) and time of this image, this field is probably full cover, well-watered, winter wheat.
68
Fig.3. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (left, 284K:Blue to 295K:Red) and the false color image (right). The coordinates of the inquire cursor are x/y = 439452/194132)
There is another center pivot irrigated field in Fig.3 to the left of the candidate field. The false color image indicates that the intensity of the vegetation is not very strong. It is interesting that the Ts image shows some variation of temperature in the field, and it looks as if the image has captured the impact of wet vegetation near the center pivot lateral. It is expected that the temperature be relatively homogeneous within a “dense” wheat field irrigated by center-pivot, and that the impact of the wet soil/vegetation be more apparent in a field having less than full-cover. The surface temperature of the second field near the center pivot lateral is 290 K, which is slightly higher than for the dense wheat field. This is understandable, because thedue dense, well-watered wheat fieldThe should have higher ET than a sparse, well-watered field to higher transpiration rates. higher ETa has decreased the surface temperature by 0.5 K. Next, we will look at a different region within the area of interest. In the center region, it appears that there are not as many well-watered agricultural fields (Figure 4).
Fig.4. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (left, 284K:Blue to 295K:Red) and the false color image (right). The coordinates are x/y = 488062/151266 )
69
The circled field in Fig.4 (486446/147416) appears to be the only good candidate field in this area. The temperature is 289.5 K. We are evaluating an April image. If this were a July image, we would not select this field as a candidate because there would be many well-watered fields at this time and this field, with a lower temperature than surrounding fields, would probably represent an “extreme” condition. Since this is an early April image, we do not expect many of the fields to be in full cover and this field may be a good candidate for the “cold” pixel. The temperatures of the two “cold” pixel candidates from two different parts of the image are the same. This is good confirmation that we have located an expected temperature condition representing typical full cover vegetation that is well-watered. It is common for similar temperatures to occur in well-watered, well-developed agricultural fields. However, it is not necessary that the “cold” pixel candidates in different areas have similar temperatures. The temperatures of “cold” pixel candidates can vary by region due to regional climate differences within an image, vegetation development, and soil type. However, if the regional weather is known to vary significantly within an image, the SEBAL user should carefully investigate how to best select the “cold” pixel, because one basic assumption in SEBAL is that climate conditions are uniform within the area of interest. Where weather varies, one may consider splitting the image into two or more parts and apply SEBAL to each part separately, with unique cold and hot pixels for each part. In the northeast region5). of The our example image, there arethat some “cold” bare fieldsirrigated (the blue circled fields in Figure false color image shows many center-pivot bare fields in this area were just starting to be irrigated. One could select such a wet bare field; however, it is best if the “cold” pixel is within a well-watered field having a full cover crop. This is because we relate ET for the cold pixel to ET r representing ET from a full cover alfalfa crop. There is only one field in this area that meets this criterion in Figure 5 (circled green).
70
Fig.5. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (left, 284K:Blue to 295K:Red) and the false color image (right). The coordinates of the inquire cursor are x/y = 538010/186916) o higher than the two other The surface temperature in this field is 291.1 K, which is 1.5 candidate fields selected earlier. If this image was for July or early August rather than early April, more fields would be well-watered. In this case, we would not select this field because this may represent an extremely cold condition (the coldest in the region). In our April situation, however, the false color image tells us that this field is one of the few fully vegetated fields in the area and that it may be a good candidate for a cold pixel.
Finally, we chose the candidate field from Figure 4 with T cold = 289.5 K for our “cold” pixel. This field is close to the weather station (20km) and should experience similar weather conditions. If the “cold” pixel is far from the weather station, it is possible that the field is “cold” because of varying weather conditions. After selecting the “cold” pixel, we record Tcold and its coordinates for future reference.
B. “Hot” Pixel Selection Philosophy of the Hot Pixel. The “magic” of SEBAL is in its prediction of sensible heat (H) for each pixel of an image. H is needed to solve forλET at each pixel using the energy balance equation λET = Rn – G – H. Values for H are “distributed” across an image in SEBAL according to the surface temperature (Ts). This is done using a “dT vs. Ts” function where dT is the difference between the air temperature very near the surface (at 0.1 m above the zero plane displacement height, d) and the air temperature at 2 m above the zero plane displacement height. SEBAL presumes a linear change in dT with sT. The linear equation for dT vs. T is developed by using the dT values for the cold and hot pixels. The s dT at the “hot” pixel is determined from H, as described in the main text, by assuming that H = Rn – G. Therefore, the hot pixel should beone where there is no ET.
71
Instructions
The following steps are necessary for the selection of the “hot” pixel: 1. The selection of the “hot” pixel follows the same procedure as for the “cold” pixel. It is more difficult, however, because there is a broader range of temperatures for “hot” pixel candidates to select from. 2. assume The “hot”there pixelis should be located in a dry and bare agricultural fieldwherethat oneone cannot no evapotranspiration taking place. It is recommended use a hot desert area, an asphalt parking lot, a roof, a south facing mountain slope, or other such extremely hot areas, since the relationship between dT and sTmay not follow the same linear relation as for agricultural and bare agricultural soils. In addition, the prediction of soil heat flux, G, is less certain for man-made structures and desert. Desert is difficult because the thin canopy creates a two-source system for heat transfer (canopy and underlying soil). 3. The “hot” pixel should have a surface albedo similar toother dry and bare fields in the area of interest. It should have a LAI in the range of 0 to 0.4 (corresponding to little or no vegetation). In the Mountain Model, the “hot” pixel should have a slope of less than one percent. As with the “cold” pixel, it should be located near the weather station. 4. If there has been precipitation within the pastthree or four days prior to the image, then it is possible that the hot pixel may exhibit some ET. In this case, it cannot be assumed that the H = Rn – G at the hot pixel. Instead, H must be predicted as R n– G – ETbare soil, where ETbare soil is predicted by operating a daily soil water balance model of the surface soil using ground-based weather measurements. An example model is included in Annex 8 of FAO-56 available at: http://www.fao.org/docrep/X0490E/X0490E00.htm A downloadable spreadsheet from Annex 8 is available from: http://www.kimberly.uidaho.edu/water/fao56/index.html Generally, the ET from a bare soil can be roughly estimated to be equal to 0.8 ET r, 0.5 ETr, 0.3 ETr, 0.2 ETr, and 0.1 ETr, 1, 2, 3, 4, and 5 days following a substantial rain event (15 mm or more). Illustrative Example (Landsat 7 ETM+, 4/8/2000)
In the following image, a simple daily water balance model predicted that bare soil in this area is totally dry so that we can assume zero ET for the “hot” pixel.
72
We open the image or subset image and select true and false color. We open the sTimage using Raster Options-Pseudo Color and link the two viewers. We go to Raster – Attributes in the Ts viewer to view the attribute table shown below:
The hottest pixel in the image is around 316 K. However, the purpose is not to select an extreme “hot” pixel but to select a “dry” pixel where ET is zero. The problem is that the surface temperatures of “totally dry surfaces” are not necessarily uniform. For example, Ts of a dry black asphalt road might be extremely high while T s of a dry desert might be lower but still quite high. T s of a very dry agricultural field might be moderately high. ET from these three different surfaces could all be zero or very close to zero. The selection of an extremely high temperature pixel for the “hot” pixel could cause SEBAL to predict positive ET for other “dry “ pixels. For example, if the “cold” pixel is in a field having Ts = 18o C and our “hot” pixel is selected in the middle of a desert with sT= 43o C, SEBAL would regard any totally dry and fallow agricultural field having sT= 30o C as a field with some moisture. This is because the dT for this dry agricultural field would be predicted to be less than that for the hot pixel based on our dT vs. T s function, and therefore H would be under predicted. Therefore, it is best to limit the land surface condition for candidate hot pixels to agricultural fields. Our first step in selecting a “hot” pixel is to find a representative, totally dry agricultural field in the area of interest. It is important to check precipitation records in order to get an idea of how wet an area is. One can apply a water balance model to predict the residual moisture in non-irrigated bare soil. For our example, many fallow fields are evident in the image and
73
weather records indicate that there was no rainfall prior to the image time. We can then assume that there should be many dry and “hot” agricultural fields in the area. We do not want to look at the extremely “hot” pixels so we ignore anything higher than 311 K (from the attribute table). These values represent mostly desert areas. Next, we will colorize the Ts image using a blue-red gradation (blue for 302 K and red for 310 K):
Fig.6. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (302K:Blue to 310K:Red).
It appears that our first guess for a temperature range that includes the hot pixel (302-310 K) is too high, because the red color appears only in the desert areas. There is some bluish color in the agricultural areas, so wewill change our range to300-306 K. However, we are not yet certain that the hot pixel is below 306 K so we will assign a light blue color to 307 K, yellow to7).308.5 K and pink to 310 K in addition to the blue-red gradation for 300-306 K (Figure
Fig.7. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (300K:Blue to 306K:Red, 307K:light blue, 308.5K:yellow, and 310K: pink).
74
Now, we look at the northwest area (Figure 8).
Fig.8. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (300K:Blue to 306K:Red, 307K:light blue, 308.5K:yellow, and 310K: pink). The coordinates of the cursor are x/y = 440393/175789.
In the area shown in Fig.8, Ts of the yellow-circled field is around 302.7 K. Actually, the values for Ts in this specific field vary from 301.7-303.2 K, but we select 302.7 K as the representative value. `We assume that some fields in this area are very dry. The green-circled field in Figure 8 appears to be a fallow field with a lower temperature of 301.5 K. We must decide if this field is very dry or if it has some moisture and ET. The temperature of the pink-circled field in Fig.8 is 301.6 K, which is very similar to the green-circled field (301.5 K). However, this field contains a little vegetation since it is slightly reddish in the false color image so we assume that the pink colored field is relatively dry but probably has some moisture and ET. There are two possible explanations for why the pink-circled field has a lower temperature than the yellow-circled field: (1) The green-circled field may have some surface moisture since it has a lower temperature than the yellow-circled field (which may be very dry). (2) The green-circled field may be totally dry and bare, but has a lower surface temperature than the yellow-circled field because it is surrounded by more humid fields than the other field.
75
The pink-circled field is darker than other fields in the area (see Figure 9). It is possible that this field is not a bare field but is partially covered with dark desert-type vegetation.
Fig.9. P40/R30, Landsat 7 ETM+ 4/8/2000, true color image, x/y = 442873/179798.
We decide that the best candidate for the “hot” pixel in the northwest region is the yellowcircled field (302.7 K). Next, we will look at the southeast region (Figure 10).
Fig.10. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (300K:Blue to 306K:Red, 307K:light blue, 308.5K:yellow, and 310K: pink). The coordinates of the cursor are x/y = 518565/159301. (Note that the cross-hairs point to the eastern part of the city of Burley. Cities should be avoided during anchor pixel selection).
It appears that this region is hotter than the northwest region. The light blue, yellow and pink areas are where the temperature is 307-310 K. which is much hotter than our previous “hot” pixel candidate (302.7 K) selected from the other area. This is explained as follows: 1) The regional weather or soil conditions are different between thenorthwest and southeast regions of this image. If so, we might think about separating the image into two areas, since SEBAL assumes that the regional weather is relatively uniform within the area of interest. 76
(2) The regional weather or soil conditions are the same but there is another reason why the surface temperature conditions are so different between two areas of this image. We must now decide how to proceed. First, we doubt that the northwest region is colder because of rainfall but we are not completely sure of this since the precipitation data is for the center of the area where the weather station is located. We therefore investigate T s in the center of the image where the weather station is located and where we know that there has been no observed precipitation (Figure 11).
Fig.11. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (300K:Blue to 306K:Red, 307K:light blue, 308.5K:yellow, and 310K: pink). Center of the plain (the coordinates are x/y = 466958/160477)
The temperature range for the center of the image is very similar to northwest region which indicates that our “hot” pixel candidate from the northwest region (302.7 K) might be a good selection and representative of bare soils having ET = 0. Now, we go back to the southeast region (Figure 12). If we look at the region from the center to the east, it is clear that surface temperatures increase toward the east. Since the east tends to contain more sandy soil, this temperature increase might be due to soil type.
77
Fig.12. P40/R30, Landsat 7 ETM+ 4/8/2000, Colored Ts image (300K:Blue to 306K:Red, 307K:light blue, 308.5K:yellow, and 310K: pink). Center to the east of the plain.
When we scan the entire area of interest, we find that the surface temperature in the southeast region is higher. A type of “heat island” can be recognized in the followings T image (Figure 13).
Fig.13. P40/R30, Landsat 7 ETM+ 4/8/2000, Ts image (Black and white, entire snake river plain)
78
We conclude that this “heat island” is a local phenomenon, occurring during the early spring image, and that our “hot” pixel of 302.7 K from a more typical part of the image is valid. This value for Thot is probably not a good value to apply for the southeast region. It would be best to separate this region into a separate SEBAL application. In our application, we will include this region, but we will ignore the ET estimates for this region since it is outside our area of interest. After selecting a “hot” pixel, Thot should be recorded along with its coordinates.
79
Appendix 8 Iteration Spreadsheet Used for the Atmospheric Stability Correction
The SEBAL process described in this manual includes a spreadsheet named “Calculations_for_05_H_Automatic_GMD_model.xls” that is used to facilitate the iteration process used in SEBAL for the correction of the aerodynamic resistance to heat transport (rah) due to buoyancy effects in the atmosphere. In SEBAL, first, initial values for the friction velocity (u*) and rah are computed assuming neutral atmospheric conditions. The sensible heat flux (H) is computed next. Then, using the Monin-Obukhov theory, *uand rah are computed again using the stability correction. The correctedahr is then used to compute a new dT function and a new value for H. This iteration is repeated until dT andahr at the “hot” pixel stabilize. In the SEBAL process and ERDAS Imagine Model Maker, the iterative process can be done manually, i.e., an updated Imagine model (for example, “F13_h_1.gmd”) is operated during each iteration. Alternatively, the iterative process for Hcan be done automatically for all iterations using the “F13toF24_h_automatic.gmd” model. In both cases, the spreadsheet is useful for calculating the correlation coefficients b and a of the dT function that is used for the calculation of H inside SEBAL. The spreadsheet calculates values for Hcold, dTcold, Hhot and dThot. The correlation coefficients b and a for the dT function are computed and a corrected value for aTand ρ are then calculated. As the iteration continues with stability corrections, new values for dT cold and dThot are calculated along with new coefficients, b and a. On the following page is an example of the iterative spreadsheet followed by instructions for the manual operation and the automatic operation of the process.
80
Part 1: using the manual iteration process
1. Go to the spreadsheet. 2. Enter values in the colored, marked boxes for the following parameters: a. The “blending height” forhigh level wind speed (generally 200 m) b. z1 in m (generally 0.1 m) c. z2 in m (generally 2 m) d. ETr at the time of the image in mm/hour 3. For the “cold” and “hot” pixels, enterthe following information in thecolored, marked boxes: a. The x and y coordinates of both pixels b. Elevation in meters of both pixels (elevation ofthe weather station for theFlat Model) c. The expected “ETrF” (i.e,. fraction of ETr ) at both pixels (this is generally 1.05 for the “cold” pixel and 0 for the “hot” pixel) d. Surface temperatures for both pixels (Ts) in Kelvin e. Elevation corrected surface temperatures (Ts_dem) for both pixels. (For the Flat
81
Model, Tsis used) f. Net Radiation (Rn) for both pixels in W/m2 at the time of the image (computed by SEBAL) g. Soil Heat Flux (G) for both pixels in W/m2 at the time of the image (computed by SEBAL) h. Momentum roughness length (zom) for both pixels i. Wind speed at the blending layer height (u200m) in m/s. This value is generally the same for both pixels. 4. Following the entry of the data in step 3, the spreadsheet will automatically calculate the b and a coefficients that are used in the dT = a T s + b equation. These coefficients are entered into an updated “F13_h_1.gmd” model, one set at a time (i.e., per iteration). a) Load model “F13_h_1.gmd”. b) Specify the file names for Ts and initial rah, files along with the elevation and enter the first set of coefficients b and a. c) Save the output files for h_1, airdensity_1, and the model. 5. Corrected values of u* and rah are computed using the Model Maker: a) Load model ” F14__u_star_2_rah_2.gmd.” b) Specify file names for h_1, airdensity_1, initial u*, Ts, zom, and the value for u200. c) Save output files for u_star_2 and rah_2. 6. Modify model “F13_h_1.gmd” to use new filenames and enter the second(or third, or fourth, so on) set of values for b and a.: a) Enter rah_2, and new coefficients b and a. b) Save output as h_2 and airdensity_2 c) Save model as; model “F15_h_2.gmd” 7. Corrected values of u* and rah are computed (similar to step 5). a) Use model “F14” from step 5 b) Enter file names for h_2, airdensity_2, and u_star_2 c) Save output as u_star_3 and rah_3 d) Save model as; model “F16__u_star_3_rah_3.gmd.” 8. Repeat steps 6 and 7, etc. using corrected values. 9. Iteration continues until dThot and rah for the hot pixel stabilize (usually 3 to 7 iterations). 10. Finish with a final run of the model for H. Save the output files for h and airdensity along with the model. .
82
Part 2: Using the automatic iteration process
The model “F13 to 24_h_automatic.gmd” is used with the support spreadsheet, “Calculations_for_05_H_Automatic_GMD_model.xls”. This model should be operated after F12_initial_u_star_rah.gmd. Since this one large model gives the final H value in W/m2, one can skip all the iteration models. Therefore, after running this model, only the final model “F25_et24.gmd” need be run. The following are the steps necessary to setup and run this model: 1. thesupport “cold” and “hot” pixels 2. Select Open the spreadsheet and input the necessary values Enter a) b) c) d)
values in the colored boxes for the following parameters: The “blending height” forhigh level wind speed (generally 200 m) z1 in m (generally 0.1 m) z2 in m (generally 2 m) ETr at the time of the image in mm/hour
For the “cold” and “hot” pixels, enter the following information in the colored boxes: a) The x and y coordinates of both pixels b) Elevation in meters of both pixels (elevation of the weather station for the Flat Model) c) The expected “ETrF” (i.e,. fraction of ETr ) at both pixels (this is generally 1.05 for the “cold” pixel and 0 for the “hot” pixel) d) Surface temperatures for both pixels in Kelvin (read from the sTimage) e) Elevation corrected surfacetemperatures for both pixels. (read from Ts_dem image; for the Flat Model, use Ts values) f) Net Radiation (Rn) for both pixels in W/m2 at the time of the image (read from Rn image) g. Soil Heat Flux (G) for both pixels in W/m2 at the time of the image (read from Rn image ) h. Momentum roughness length (zom) for both pixels (read from zom image) i. Wind speed at the blending layer height (u200m) in m/s. (calculated in the support spreadsheet “SEBAL Model constants for flat model.xls”)
Once the above values are input, the correlation coefficients b and a for 7 iteration steps are computed. One must make sure that the rah and dT values of the “cold” and “hot” pixels are numerically stable (i.e., they do not diverge with the iteration time). This can be checked by looking at the red-circled area in the following snapshot:
83
The blue-circled cell tells you how many iterations are recommended. The “F13to24_h_automatic.gmd” model automatically runs 7 iterations with the 7 sets of coefficients b and a. As long as the “rah” and “dT” values for “cold” and “hot” pixels converge, the number of iterations is not important. If these parameters diverge, we recommend increasing the wind speed at 200m to 4 m/s, regardless what was computed using the actual weather data 3. Setup and Run the “F13to24_h_automatic.gmd” model Complete the following steps to run the model: a) Load model “F13to24_h_automatic.gmd”. b) Define the following input values and files: • Surface temperature image • Representative elevation value • Representative Elevation value • rah_1 image • zom image • calculated wind speed (m/s) at 200m height c) Enter all seven sets of coefficients b and a into the appropriate tables. (Note that the first entry in the tables is numbered as “0”). d) Save the output files for h and airdensity along with the model.
84
Appendix 9 SAVI and LAI for Southern Idaho
Note: Equations (14) and (15) forSAVI and the LAI are not well calibrated for universal use with all vegetation types, and some uncertainties remain on their use. They are, however, utilized in SEBAL for two reasons: (1) LAI might be a better index than NDVI, since it has a physical meaning, and (2) LAI might be a better index for estimating the soil heat flux (G) and the momentum roughness length (z om).
Some of the development or validation of these equations in Idaho is shown below: Change of SAVI number w ith various L values for 5 soil / dead g rass conditions in Southern Idaho, 1989
0.15
0.01
0.14
0.009
0.13 0.008
0.12
I 0.11 V 0.1 A S0.09
0.007
0.006
0.08
0.005
0.07 0.004
0.06
n io t a i v e D d r a d n a t S
#1 #2 #3 #4 #5
0.05
0.003 0
1 . 0
2 . 0
3 . 0
4 . 0
5 . 0
6 . 0
7 . 0
STDEV
L SAVI vs. L, where L is the coefficient in the SAVI equation for five soil or dead grass conditions in southern Idaho; #1: Dry bare field with bright color (4/18/89), #2: Wet bare field with dark color (5/4/89), #3 to #5: Totally dead grassland from desert area with different soil color. The standard deviation (STDEV) between the five soil locations or conditions becomes minimum when L = 0.1. This indicates that L = 0.1 tends to reduce the impact of soil wetness on the vegetation index. The value of L for SAVIID was chosen as 0.1 due to this analysis.
85
A value of L for a particular area can be derived from an analysis of multiple images for the area. Using three or more images, look at some pixels of vegetation that is not likely to change much (for example; sage brush during mid to late summer, conifer trees, etc.). Compute SAVI using a range of values for L and plot SAVI vs. image date. Select the L that provides for a consistent value of SAVI over time. Compareison of calculated and measured LAIs 1989 Sugar Beets 6 5 4
I A3 L 2 1
Estimated by SEBAL Measured LAI (Beet)
0 9 8 / 2 / 4
9 8 /
9 8 /
9 8 /
9 8 /
5
7
8
0 1
2 2 /
1 1 /
0 3 /
9 1 /
9 8 / 8 / 2 1
Appendix 10
G/Rn Computations for Water Surfaces and Snow
SEBAL uses an empirical function (equation 24) to compute G/R n based on surface temperature, surface albedo, and the NDVI. For surfaces such as water, snow, and wetlands other methods for defining G/R n are needed. A. G/Rn for Water Surfaces The heat balance at water surfaces works very differently than at land surfaces. Generally, albedo values for clear water are very small (less than .05) when the sun angle is more than 10o above the horizon. Shortwave radiation penetrates into a water body according to the transparency of the water and is absorbed at a range of depths below the surface where it is converted into heat (G). For oceans and rivers, a large amount of energy is transferred horizontally with the currents, and the vertical heat balance equation is no longer consistent. Hence, the equation for G/R n for a water surface is difficult to define. The function for G/Rn that is most suited for a specific water surface must be evaluated.
For the inland states such as Idaho, the primary water surfaces are lakes where the depth of solar radiation penetration will vary with the sun angle and the turbidity of the water body. In the 2000 Idaho project, a large, deep lake (Bear Lake) was the primary water surface in the study area. Therefore, a deep lake condition was focused on during development of the
86
equation for G. The following are examples that may provide ideas on how to compute G for lakes: 1. G equation for deep and clear lakes Yamamoto and Kondo reported the seasonal trends in the heat balance for a deep lake from a set of actual measurements (Yamamoto and Kondo, 1968). Figure 10.1 shows a graph of the monthly averages for Rn and G values that they measured in a deep lake (mean depth was 21m) in Japan. This graph indicates how the ratio G/R n varies with the time of year. In the spring and early summer it is large since the water is cold and the air is warm. In the fall it is small since the water is warm and the air is cool. One can note also that for a deep, clear body, the net daily value for G can be strongly positive or negative due to change in heat storage, whereas, it is generally assumed for vegetation that the net G for a 24-hour period is nearly zero (i.e., no net change in storage for the day).
Figure 10.1 The plot of the monthly averages of R n and G values measured in a deep lake (mean depth of 21m) in Japan (Yamamoto and Kondo, 1968).
In our 2000 project (Tasumi et al., 2000), we separated the year into two periods; JanuaryJune and July-December. We derived a Rn-G function for each period using Yamamoto’s measurements. Figure 10.2 is for January-June. One extreme data point in January was rejected. Figure 10.3 is for July-December. One extreme data point in July was rejected.
87
y = 0.8694x - 35.286 2
150
R = 0.8174
) 100 2 m / 50 W ( G 0 0
50
-50
100
150
200
Figure 10.2. A Rn-G curve and measurements for estimating G for water during the January to June period (based on Yamamoto and Kondo, 1968).
Rn (W/m2)
y = 1.0656x - 114.28 50
) 0 2 /m-50 W ( G-100 -150
2
R = 0.9472 0
50
100
150
Figure 10.3. A Rn-G curve and measurements for estimating G for water during the July to December period (based on Yamamoto and Kondo, 1968)
Rn (W/m2)
This data represents the monthly averages and G may be larger if the day time is focused on. In addition to Yamamoto’s work, Amayreh (1995) measured G for Bear Lake, Idaho/Utah using eddy covariance and heat storage change and reported that daily values of G during the summer and fall could be predicted using: G = 0.984Rn – 62 . This equation does not take into account the seasonal difference in G since Amayreh’s work was limited to summer and fall. However, Amayreh’s measurements are similar to those of Kondo’s in both magnitude and trends and this simple prediction approach based on n R is therefore assumed to be valid for deep clear water bodies in general. The following equations were developed for estimating G for water surfaces based on Yamamoto’s and Amayreh’s work: From July to December, the instantaneous G during midday (and at the satellite image time): Gwater = R –n90
(10.1)
Gwater = R–n100
(10.2)
and for 24-hours:
From January to June, the instantaneous G during midday (and at the satellite image time): Gwater = 0.9 –R40 n
88
(10.3)
and for 24-hours: Gwater = 0.9 –R50 n
(10.4)
The penetration of solar radiation into water decreases as the depth of the water body decreases and/or as the turbidity of the water increases. Highly turbid bodies of water and very shallow bodies of water will have a value of G that is close to the value of nRand the offsets of the above equations will be reduced or may completely disappear. 2. G equation for shallow (less than 2m) or turbid lakes The value of G for shallow lakes is estimated to be between the value for deep lakes and for wetlands. This results in the following equation developed from Table 10.1: for instantaneous at midday: Gwater = 0.50Rn and for 24-hours: Gwater = 0
89
(10.5)
Table 10.1. Calculated G/Rn ratio for deep lake and wetland for various Rn ranges. For deep lakes, G/Rn equation for Bear Lake was used. For Wetlands, Burba et al (1999: G = -51 + 0.41Rn) was used.
Rn 200 250 300 350 400 450 500 550 600
G/Rn Deep Lake Wetland 0.63 0.16 0.69 0.21 0.73 0.24 0.76 0.26 0.79 0.81 0.82 0.83 0.84
0.28 0.30 0.31 0.32 0.33
B. G/Rn for Snow Surfaces In the application of SEBAL to Idaho, the value of G for snow surface was predicted for daytime periods as:
Gsnow = 0.5 Rn
(10.6)
This is a very crude estimate for G and assumes that one half of the net radiation incident on the snow surface penetrates in the form of light and is absorbed into the snow mass as G. For 24 hours, the average G was assumed to be zero. Additional research and literature review is needed in this area.
Estimated LAI values from eight TM images using equation (15) compared with actual LAI values measured by Dr. J.L.Wright, USDA/ARS, in a sugar beet field at Kimberly, Idaho 1989.
Equation (15), when used with SAVIL=0.1, closely estimates the measured LAI of sugar beets at Kimberly, especially during the crop-growing season. The Idaho application has shown that the LAI equation (15) is quite sensitive to SAVI after LAI > 3. This confirms the known phenomenon that the NDVI (which is similar to SAVI) “saturates”, i.e., it does not increase, once the LAI reaches a certain amount (Kondoh, 1999). The sensitivity problem of LAI values greater than 3 does not affect the G or H estimation. The accuracy of equation (15) has been found to be extremely sensitive to the value for L used to compute SAVI. Field verification is needed fornew areas or for new vegetation.
90
Appendix 11 The Stability of the Atmosphere Air has three stability conditions; unstable, neutral, and stable. Stability conditions must be considered during the computation of sensible heat flux (H) because they affect the aerodynamic resistance to heat transport (rah). In SEBAL, a neutral atmospheric condition is initially assumed and a stability correction is later applied using the Monin-Obukhov length (L) as the indicator. The three stability conditions for air are illustrated in Figure 11.1. If a dry air mass at 10o C is suddenly moved upward 100 meters, its temperature will decrease about 1o C due to the lapse effect. If the surrounding air temperature is lower than this air mass, an upward force on the air mass will occur. This condition is called “unstable”. If the surrounding air temperature is the same as the air mass, no force will occur and the condition is called “neutral”. If the surrounding air temperature is higher than the air mass, a downward force on the air mass will occur. This condition is called “stable”. Generally, air temperature decreases by about 0.65 oC when elevation increases by 100 m under neutral stability conditions. In places where positive sensible flux because (H) occurs, temperature with elevation becomesheat smaller thethe air change mass is of heated by the positive H. In this case, vertical air movement is easier and therefore aerodynamic resistance becomes smaller as H becomes higher. In general, air is in a neutral condition at a well-watered agricultural field at noon and an unstable condition at a dry bare field at noon. The stable condition will most likely occur at night and sometimes in the afternoon in irrigated areas surrounded by desert.
91
Tendency of vertical movement of air at near land surface, o
assuming 6.5 C/km air lapse rate
o
o 19.35 C 19.7 C
o
19.35 C
o
19.35 C
o
19.35 C
o
19.1 C
100m o
20 C
o
20 C
o
20 C
o
20 C
o
20 C
S tabl e
N eutra
U nst abl e
o
20 C
: Direction of Force : Direction of Sensible Heat Flux
Relative shapes of eddy at near land surface
U nst abl e
N eutra
Figure 11.1. Three Conditions of Atmospheric Stability
92
S tabl e
Appendix 12 The SEBAL Mountain Model SEBAL provides accurate estimations of ET for relatively flat agricultural areas. For use in mountainous areas where there is significant relief and a wide range of slopes and aspects, the SEBAL Mountain Model was developed. The SEBAL Mountain Model contains modifications that correct for slope, aspect, and elevation in the computations. The 24-hour solar radiation is typically much higher on a south slope than on a north slope in the Northern Hemisphere. The result of not making the corrections for slope, aspect, and elevation is that a higher than actual ET will be computed for north slopes and high elevations because SEBAL incorrectly interprets the cooler surface temperatures as indicating “wet” (therefore cool) areas. The cooler surface temperatures on north facing slopes and high elevations is not due to “wetness” however, but is from the temperature lapse caused by higher elevation and from less incident solar radiation caused by the aspect. Predictions of air flow and wind speed over mountainous areas are much more difficult. The effects from the orthographic drainage of air caused by cooling, the acceleration of air streams passing over mountains due to the venturi effect, and the impacts of drag due to undulating topography all combine to make wind speed prediction difficult and complicated. The following describes the modifications that are made to the SEBAL Flat Model in order for it to be used in mountainous areas (SEBAL Mountain Model): A. Data Requirements A Digital Elevation Model (DEM) is required to provide elevation data. The units of the DEM should be in meters with a pixel size the same as that of the satellite data (30m× 30m). As in the flat model, a layered image must first be created followed by a smaller subset image if desired (Appendix 4). A DEM subset image must then be created which is the same size as the subset image. To do this, one must go to the ERDAS Data Prep-subset image tool and input the DEM from a disk and enter the coordinates of the subset image. The output will be the DEM subset image. B. Slope and Aspect Images The Mountain Model requires slope and aspect information for its computations. The slope and aspect for each pixel are derived using the DEM subset image and the ERDAS Interpreter–topographic analysis–slope/aspect tool. The slope and aspect images are then saved as the output. C. Calculation of the Cosine of the Solar Incidence Angle (cos ) The solar incidence angle is the angle between the solar beam and a vertical line perpendicular to the land surface. In the flat model, we assume that the land surface is horizontal and the calculation of cosθ is very simple and is a constant over the area of interest. In the Mountain Model, cosθ is different for each pixel depending on the slope and
93
aspect of the land surface. SEBAL uses the following equation for computing cos θ (Duffie and Beckman, 1991): cos θ
= sin(δ ) sin(φ ) cos( s) − sin(δ ) cos(φ ) sin(s) cos(γ ) + cos(δ ) cos(φ ) cos( s) cos(ω ) + cos(δ ) sin(φ ) sin( s) cos(γ ) cos(ω ) + cos(δ ) sin(φ ) sin( s) sin(ω )
(12-1)
where;
δ= φ=
declination of the earth (radians; positive in summer in northern hemisphere) latitude of the pixel (radians; positive for northern hemisphere) s= slope (radians) where; s=0 is horizontal and s=π/2 is vertical downward (s is always positive and represents a downward slope in any direction) γ = surface aspect angle (radians) where; γ = 0 for due south, γ = -π/2 for east, γ = +π/2 for west and γ = ± π for north (Figure 12.1). ω= hour angle. ω = 0 at solar noon, ω is negative in morning andω is positive in afternoon
Figure 12.1. The definitions of slope and aspect in SEBAL
In order to solve equation (12-1), the following steps are taken: 1. The sine and cosine of the slope (s) and aspect (γ) are calculated. This is done with model_01_sin_cos_slope_aspect. The slope and aspect image files are input and the sine and cosine of the slope/aspect are output and saved along with the model. 2. The declination (δ), latitude (φ), and hour angle (ω) are computed. The declination (in radians) is computed as:
δ = .409 sin{(2π/365 × DOY) – 1.39}
(12-2)
The latitude (in radians) is computed as:
φ = latitude in degrees × π/180
(12-3
The hour angle (in radians) is computed as:
ω = π/12{(t + 0.06667(Lz - Lm) + Sc) - 12}
94
(12-4)
where; t is the standard clock time for the satellite overpass time (hours; i.e. for 14:30, t = 14.5), Lz is the longitude of the center of the local time zone (degrees west of Greenwich; i.e., Lz = 75, 90, 105, and 1200 for the Eastern, Central, Rocky Mountain, and Pacific time zones respectively), L m is the longitude of the center of the satellite image (degrees west of Greenwich), and S c is a seasonal correction for solar time (hours). The seasonal correction for solar time is computed as: Sc= 0.1645 sin(2b) – 0.1255 cos(b) – 0.025 sin(b) b = 2π(DOY – 81)/364
(12-5)
(12-6)
3. Import and run model_02_cos_theta. Inputs arethe sine and cosine of the slope/aspect, the declination, the latitude, and the hour angle. The output file for the cosine of the incidence angle is saved along with the model. In model_02, an initial value for cos θ is computed and is then divided by the cosine of the slope to correct it to a horizontal equivalent. This assures that all of the energy fluxes computed by SEBAL will be analyzed on a horizontal equivalent basis. D. Reflectivity In the model for reflectivity, the input of cosθ is a file rather than a constant. E. Transmissivity The transmissivity is no longer a constant but varies for each pixel as a function of the elevation. It is computed with the DEM subset image as the input. F. Surface Albedo The surface albedo is computed using the transmissivity file as input. G. Incoming Shortwave Radiation The incoming shortwave radiation (RS↓) is computed for each pixel since cosθ is now variable for each pixel. H. Surface Temperature Generally, air temperature decreases 6.5oC for each 1 km elevation increase under neutral stability conditions. Since surface temperatures are in strong equilibrium with air temperature, one can usually observe similar decreases in surface temperature.
During the prediction of the near surface temperature difference (dT), SEBAL assumes a linear relationship between dT and the surface temperature (T s). The surface temperature, however, needs to be adjusted to a common reference elevation for accurate prediction of dT. Otherwise, high elevations that appear to be “cool” may be misinterpreted as having high ET. Therefore, in the Mountain Model, a “lapsed” (artificial) surface temperature map is created for purposes of computing dT by assuming that the rate of decrease in surface temperature due to elevation increase is the same as that for a typical air profile. The
95
elevation data is provided by the DEM. This“artificial” lapse-adjusted surface temperature is referred to as a DEM corrected surface temperature (T s-dem). The DEM corrected surface temperature is calculated by the following equation:
Ts _ dem = Ts + 0.0065∆z
(12-7)
where; ∆z is the elevation of each pixel minus the elevation of adatum (meters). The term ∆z is positive if the elevation of a pixel is higher than the datum. The representative elevation for the image, usually at the weather station, is used for the datum elevation. I. Outgoing Longwave Radiation The outgoing longwave radiation (RL↑) is computed as in the Flat Model using the uncorrected surface temperature (Ts). J. “Cold” and “Hot” Pixel Selection The corrected surface temperature (Ts-dem) is used for the selection of the two “anchor” pixels. The Ts image shows high and low pixels rather than wet or dry pixels because of the elevation effects. For higher than datum pixels, the lapse correction is added ands-dem T is higher than Ts. Using the Ts-dem data allows us to find truly wet and dry pixels. The methods described in Appendix 7 are followed except that the T s-dem image is used rather than the Ts image. K. Incoming Longwave Radiation The incoming longwave radiation (RL↓) is computed since the transmissivity of air is no longer a constant for all pixels. The reference temperature (T cold) is adjusted for the lapse
effect in a similar way as in equation (12.7): Tcold(each pixel) = Tcold + 0.0065∆z
(12-8)
where; ∆z is the elevation of the “cold” pixel minus the elevation of each pixel. The elevation of the “cold” pixel is taken from the DEM subset image by entering its coordinates. L. Net Surface Radiation The net surface radiation (Rn) is computed in model_14 similar to the Flat Model except that incoming longwave radiation (RL↓) and incoming shortwave radiation (RS↓) are now files rather than constants. M. Soil Heat Flux The soil heat flux (G) and G/Rn are computed as in the Flat Model using the uncorrected surface temperature (Ts).
96
N. Momentum Roughness Length The momentum roughness length (zom) is initially computed as in the Flat Model using a land use map and LAI as input. It is then adjusted for the slope using the following equation:
zom(mnt) = zom × [1 + (slope - 5)/20]
(12-9)
O. Initial Friction Velocity and Aerodynamic Resistance to Heat Transport The wind velocity at 200 meters (u200) is adjusted for elevation effects before computing the friction velocity (u*) and the resistance to heat transport (rah). The elevation effects on wind
speed are complicated and difficult to quantify. A mountain wind speed weighting coefficient (ϖ) is used to adjust u200 as follows:
ϖ = 1 + 0.1[ ({Elevation – Elevationstation}/1000)]
(12-10)
where; Elevation is the elevation of each pixel and Elevation station is the elevation where the wind speed is measured. u200 for each pixel is then adjusted by multiplying it byϖ. P. Sensible Heat Flux The sensible heat flux (H) is computed as in the Flat Model. The uncorrected surface temperature (Ts) is used to compute a new air density. The corrected surface temperature (Ts-dem) is used for the computation of dT. Q. Atmospheric Stability Correction The iterative process for the atmospheric stability correction is similar to the Flat Model and uses the same spreadsheet. R. Applications and Precautions For The SEBAL Mountain Model
The Mountain Model is intended to improve the computation of components of the energy balance for sloping areas. However, some uncertainties remain, so that while the Mountain Model does provide improved estimates of ET in mountainous areas, the ET estimates in mountains have higher uncertainty thanthose for flatter areas. When integrated over a sub-basin scale, the ET from the SEBAL mountain model is probably reasonably accurate. This is because some of the uncertainties (due to unknown wind speed and direction for example) tend to cancel. At the present time, the “SEBAL mountain model” is offered by the University of Idaho for a purchase fee. The fee is intended to recoup development costs for the model and to support user inquiries and needs, which may be extensive, due to the higher level of sophistication in the mountain model.
97
98