ABSTRACT Title of Dissertation: MAPPING PHOTOSYNTHETICALLY ACTIVE RADIATION (PAR) USING MULTIPLE REMOTE SENSING DATA Tao Zheng, Doctor of Philosophy, 2007 Dissertation directed by: Dr. Shunlin Liang, Professor Department of Geography Incident Photosynthetically Active Radiation (PAR) is an important parameter for terrestrial ecosystem models. Presently, deriving PAR using remotely sensed data is the only practical approach to meet the needs for large scale ecosystem modeling. The usefulness of the currently available PAR products is constricted by their limited spatial and temporal resolution. In addition, the applicability of the existing algorithms for deriving PAR using remotely sensed data are limited by their requirements for external atmospheric information. This study develops new algorithms to estimate incident PAR using remotely sensed data from MODIS (Moderate Resolution Imaging Spectroradiometer), GOES (Geostationary Operational Environmental Satellite), and AVHRR (Advanced Very High Resolution Radiometer). The new PAR algorithms differ from existing algorithms in that the new algorithms derive surface properties and atmospheric optical properties using time-series of at-sensor radiance without external atmospheric information. First, a new PAR algorithm is developed for MODIS visible band data. The validity of the algorithm?s underpinning theoretical basis is examined and associated errors are analyzed in light of their impact on PAR estimation accuracy. Second, the MODIS PAR algorithm is adapted to AVHRR in order to take advantage of the long data acquisition record of AVHRR. In addition, the scaling of remote sensing derived instantaneous PAR to daily PAR is addressed. Last, the new algorithm is extended to GOES visible band data. Two major improvements of GOES PAR algorithm over that of MODIS and AVHRR are the inclusion of the bi-directional reflectance distribution function for deriving surface reflectance, and the procedure for excluding cloud-shadowed pixels in searching for observations made under clear skies. Furthermore, the topographic impact on PAR is accessed and corrected. To assess the effectiveness of the newly developed PAR algorithms, validation efforts have been made using ground measurements made at FLUXNET sites. The validations indicate that the new PAR algorithms for MODIS, GOES, and AVHRR are capable of reaching reasonably high accuracy with no need for external atmospheric information. This work is the first attempt to develop a unified PAR estimation algorithm for both polar-orbiting and geostationary satellite data. The new algorithms developed in this study have been used to produce incident PAR over North America routinely to support the North America Carbon Program. MAPPING PHOTOSYNTHETICALLY ACTIVE RADIATION (PAR) USING MULTIPLE REMOTE SENSING DATA By Tao Zheng Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2007 Advisory Committee: Dr. Shunlin Liang, Chair Dr. John Townshend Dr. Rachel Pinker Dr. Ralph Dubayah Dr. Si-Chee Tsay ? Copyright by Tao Zheng 2007 ii TABLE OF CONTENTS LIST OF TABLES............................................................................................................. iv LIST OF FIGURES ............................................................................................................ v Chapter 1 Introduction ........................................................................................................ 1 1.1 Importance of PAR ............................................................................................. 1 1.2 Current methods for PAR estimation.................................................................. 3 1.2.1 Statistical Models........................................................................................ 3 1.2.2 Physical models .......................................................................................... 4 1.2.3 Satellite based PAR estimation methods..................................................... 6 1.3 Issues with the existing remote sensing based PAR estimation methods........... 9 1.3.1 The conversion factor between SW irradiance and PAR............................ 9 1.3.2 Diffuse and direct PAR ............................................................................ 10 1.3.3 Topographic effect on PAR...................................................................... 10 1.3.4 The impact of snow and ice on PAR estimation........................................ 12 1.5 Objectives of this study..................................................................................... 13 1.6 Structure of the dissertation .............................................................................. 15 Chapter 2 Estimating PAR using MODIS data................................................................. 16 2.1 Introduction............................................................................................................. 16 2.2 The theoretical basis of the new PAR algorithm .................................................... 17 2.2.1 Top of Atmosphere (TOA) radiance................................................................. 18 2.2.2 Spectral downward flux at surface level.......................................................... 19 2.2.3 Exploring the relation between TOA radiance and PAR................................. 20 2.2.4 MODTRAN4 simulation and PAR/TOA radiance relation....................... 22 2.3 MODIS PAR algorithm description........................................................................ 27 2.3.1 Generating Look-Up Tables ............................................................................ 29 2.3.2 Searching look-up table to calculate incident PAR ......................................... 31 2.4 Error budget analysis .............................................................................................. 32 2.5 The uncertainty associated with atmosphere model, ozone, and water vapor........ 35 2.5.1 Uncertainty caused by variation in atmospheric ozone concentration ........... 38 2.5.2 Uncertainty caused by variation in atmospheric water vapor......................... 39 2.5.3 Uncertainty caused by variation in atmospheric gases composition............... 41 2.6 Validation................................................................................................................ 42 2.7 Calculating Daily PAR from MODIS derived instantaneous PAR ........................ 51 2.8 Mapping incident instantaneous PAR over the Greater Washington D.C. area ..... 52 2.9 Summary................................................................................................................. 53 Chapter 3 Estimating PAR from AVHRR data ................................................................ 55 3.1 Introduction............................................................................................................. 55 3.2 The new AVHRR PAR algorithm .......................................................................... 56 3.3 Error budget analysis .............................................................................................. 58 3.4 AVHRR visible band Data...................................................................................... 59 3.5 Scaling instantaneous PAR to daily average PAR.................................................. 60 3.6 Validation................................................................................................................ 62 3.7. Summary and discussion........................................................................................ 67 Chapter 4 Estimating PAR from GOES data.................................................................... 68 4.1 Introduction............................................................................................................. 68 iii 4.2 Algorithm description ............................................................................................. 70 4.2.1 GOES-12 visible band at-sensor radiance ...................................................... 70 4.2.2 Examine the relation between PAR and GOES visible band TOA radiance ... 72 4.2.3 Excluding cloud-shadowed observations......................................................... 74 4.2.4 Fitting BRDF model to calculate surface reflectance ..................................... 75 4.2.5 Error budget analysis ...................................................................................... 78 4.3 Validation................................................................................................................ 78 4.4 Correcting topographic impacts on PAR ................................................................ 83 4.4.1 Digital Elevation Dataset ................................................................................ 84 4.4.2 Angular effects on direct PAR.......................................................................... 86 4.4.3 Sky view factor for diffuse PAR ....................................................................... 88 4.4.4 Reflected direct and diffuse flux from neighboring terrains............................ 89 4.4.5 The results of the topographically corrected PAR........................................... 90 4.5 Discussion............................................................................................................... 91 4.6 Summary................................................................................................................. 94 Chapter 5 Conclusions ...................................................................................................... 96 5.1 Estimating PAR using MODIS data ....................................................................... 96 5.2 Estimating PAR using AVHRR data ...................................................................... 98 5.3 Estimating PAR using GOES data.......................................................................... 99 5.4 Combining MODIS, AVHRR, and GOES data for daily PAR estimation........... 101 5.5 Issues for future research ...................................................................................... 103 Reference ........................................................................................................................ 105 iv LIST OF TABLES Table 2-1 Solar and satellite viewing angles used to form illuminating/viewing geometry for MODTRAN4 parameterization........................................................................... 23 Table 2-2 Cloud extinction coefficients (km -1 ) at 550 nm, thickness and based heights . 25 Table 2-3 Aerosol Types and Visibility (km)................................................................... 25 Table 2-4 Look-up table that links TOA radiance to atmospheric profile index. 0 ? : solar zenith angle, ? : sensor zenith angle, ? : relative azimuth angle. ? : atmospheric albedo........................................................................................................................ 30 Table 2-5 Look-up table that links PAR with atmospheric profile index 0 ? : solar zenith angle, ? : atmospheric spherical albedo. .................................................................. 31 Table 2-6 Water vapor concentration and the integrated transmittance over the PAR region ........................................................................................................................ 40 Table 2-7 Summaries of regression analysis at six validation sites.................................. 50 Table 3-1 Validation results of instantaneous PAR at three FLUXNET tower sites........ 66 Table 3-2 Validation results of daily integrated PAR at three FLUXNET tower sites. ... 67 Table 4-1 Summary of site locations and validation statistics at the four validation sites. ................................................................................................................................... 83 Table 5-1 Statistical summary of validation results of instantaneous PAR estimated using MODIS, AVHRR, and GOES. N is number of samples. ....................................... 101 v LIST OF FIGURES Figure 2-1 Relationship between MODIS blue TOA radiance and PAR from MODTRAN4 simulations under various atmospheric conditions and solar zenith angles. ....................................................................................................................... 27 Figure 2-2 Procedure of the new MODIS PAR algorithm ............................................... 33 Figure 2-3 Relationship between MODIS band 3 (blue band) TOA radiance and PAR from MODTRAN4 simulations. The figure in the right panel represents the relation from results of all atmospheric conditions. The figure in the left panel represents the relation from only one cloud type (altostratus) and one aerosol model (rural). ....... 34 Figure 2-4 Error caused by the linear interpolation scheme of look-up tables entries generated with selective atmospheric conditions...................................................... 37 Figure 2-5 Change of PAR with the change of atmospheric volume ozone. Atmosphere condition is represented by mid-latitude summer model with rural aerosol model with visibility of 50 km, altostratus cloud with cloud extinction coefficient of 5 km -1 , water vapor concentration of 2.0g/cm2. Visible region broad band visibility is set at 0.05. Solar zenith angle is 0 o ..................................................................................... 39 Figure 2-6 Impacts of atmospheric water vapor on PAR. Atmosphere condition is represented by mid-latitude summer model with rural aerosol model with visibility of 50 km, altostratus cloud with cloud extinction coefficient of 5 km -1 , volume ozone of 350 Dobson Units. Visible region broad band visibility is set at 0.05. Solar zenith angle is 0 o . ...................................................................................................... 40 Figure 2-7 Impact of atmospheric gases composition on PAR. The five error bars presents the five standard atmosphere models included in MODTRAN4: 1: tropical, 2: mid- latitude summer, 3: mid-latitude winter, 4: sub-arctic summer, 5: sub-arctic winter. Under each atmosphere model, water vapor column varies from 0 to 8 g/m 2 , and ozone column varies from 0 to 800 Dobson Units. .................................................. 42 Figure 2-8 Combined impacts of variations in atmospheric gases composition, water vapor, and ozone concentration on PAR. Water vapor column ranges from 0 to 8 g/m 2 , ozone volume ranges from 0 to 800 Dobson units. Cloud is represented by altostratus cloud with extinction coefficients of 1, 5, 20, and 50 km -1 in the figures of the four rows respectively......................................................................................... 43 Figure 2-9 Validation of the estimated instantaneous PAR at Black Hills, ND. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. ................................................................................................................. 45 vi Figure 2-10 True color composite MODIS imagery of greater Washington, D. C. area for 8 days (May 22, May 25, May 29, May 31, June 5, June 7, June 8, and June 10, 2003). ........................................................................................................................ 52 Figure 2-11 Retrieved incident PAR maps from the corresponding MODIS images in Figure 2-17................................................................................................................ 53 Figure 3-1 Relationship between AVHRR-14 band-1 TOA radiance and PAR for a range of various atmospheric profiles and solar zenith angles. .......................................... 57 Figure 3-2 Linear interpolation of PAR based on the relation between PAR and TOA radiance established from simulations using only one cloud type and one aerosol type. The left panel is the TOA radiance-PAR relation based on one cloud/aerosol type simulation, the right panel is based on multiple cloud/aerosol type. ................ 59 Figure 3-3 Number of satellite overpasses by NOAA-14 AVHRR over the North America continent for day 143 in 1996. .................................................................................. 62 Figure 3-4 AVHRR band-1 data and estimated instantaneous PAR value covering the continental U.S.. The image in the upper panel is TOA radiance at 20:05 GTM, day 225 in 2004. The image in the lower panel is the estimated instantaneous PAR for the same time. The white pixels on the TOA radiance image are clouds. The unit of TOA radiance is Wm -2 sr -1 s -1 . The unit of instantaneous PAR is umol/m 2 /s............. 64 Figure 3-5 Validations of PAR derived using AVHRR band 1 data from NOAA-14 at Walker Branch, Tennessee. The left graph is a comparison between retrieved (X- axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis). ............................... 65 Figure 3-6 Validations of PAR derived using AVHRR band 1 data from NOAA-14 at Metolius, Oregon. The left graph is a comparison between retrieved (X-axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis)............................................... 65 Figure 3-7 Validations of PAR derived using AVHRR band-1 data from NOAA-14 at Park Falls, Wisconsin. The left graph is a comparison between retrieved (X-axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis)............................................... 66 Figure 4-1 Relationship between GOES12 band-1 TOA radiance and PAR at a range of different atmospheric profiles and solar zenith angles. Surface reflectance is 0.05, sensor zenith angle and relative azimuth angle are all set to 0 o ................................ 73 Figure 4-2 An examples of the TOA radiance converted surface reflectance (triangles) and the BRDF simulated surface reflectance (crosses) during day 191-198, 2004 at Lost Creek, Wisconsin.............................................................................................. 77 vii Figure 4-3 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at Canaan Valley, West Virginia. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. ............................... 80 Figure 4-4 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at the Lost Creek site in Wisconsin. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. ............................... 81 Figure 4-5 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at the Willow Creek site, Wisconsin. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. ......................... 82 Figure 4-6 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at Metolius site, Oregon. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. ............................................. 83 Figure 4-7 USGS GTOPO30 DEM projected to Lambert Azimuthal Equal Area projection, with center of projection at 90W/30N (Longitude/Latitude). The study area covers the continental U.S. and part of southern Canada.................................. 85 Figure 4-8 Effect of topography on direct PAR based on the angle between solar illumination direction and the normal to the terrain slope. Dark colors represent pixels where low percentage of incident direct PAR is received. Bright colors represent pixels where high percentage of incident PAR is received. Note that the topographic effect on direct PAR is dependent on time of day and Julian day of year. This image is of 19:05 GMT, May 04, 2004. ........................................................... 87 Figure 4-9 Histogram of values of sky view factor covering the same area as Figures 4-8. The sky view factor ranges from 0.3 to 1.0 with mean value of 0.98 (X axis) and standard deviation of 0.019. The Y axis is the percentage value.............................. 88 Figure 4-10 GOES visible band data and topographically corrected PAR covering the study area. The 5 images on the left panel are the TOA radiance at 20:00 GMT, day 191, 192, 193, 194, 195, year 2004. The five images on the right hand panel are the topographically corrected PAR at corresponding times. The impact of topographic impact is more pronounced at absence of cloud than at presence of cloud. Mountainous areas experience more topographic impact than relatively flat areas. The white pixel on the TOA radiance images are clouds. The unit of TOA radiance is W/m2/rad/nm. The unit of instantaneous PAR is umol/m 2 /s................................ 93 Figure 5-1 Validation of the instantaneous PAR derived using MODIS data at the all six FLUXNET sites used in chapter 2. The solid line the is 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression..................................... 97 viii Figure 5-2 Validation of the instantaneous PAR derived using AVHRR data at the three FLUXNET sites used in chapter 3. The solid line is the 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression..................................... 99 Figure 5-3 Validation of the instantaneous PAR derived using GOES data at the four FLUXNET sites used in chapter 4. The solid line is the 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression................................... 100 1 Chapter 1 Introduction Incident Photosynthetically Active Radiation (PAR), defined as illuminated solar radiation at the Earth?s surface between 400 and 700 nm, is an important variable controlling terrestrial green vegetation?s biological productivity. To meet the needs of vegetation involved ecological modeling, advanced methods have been developed to estimate incident PAR using remotely sensed data. The availability of remotely sensed data acquired by Moderate Resolution Imaging Spectroradiometer (MODIS) and other sensors offers the possibility for new approaches in estimating PAR with regional to continental applicability. The primary objective of this study is to develop new algorithms for estimating PAR, using data from MODIS and other sensors, such as Advanced Very High Resolution Radiometer (AVHRR) and Geostationary Operational Environmental Satellite (GOES). 1.1 Importance of PAR Change in terrestrial biological productivity is the most fundamental measure of ?global change?. Biological productivity is the primary source for food, fiber and fuel (Running et al., 2004), and plays an important role in studying global climate and the biogeochemical circle (Myneni et al., 1995; Sellers and Schimel, 1993). The enormous spatial variability of Net Primary Productivity (NPP) at global scale is likely to change with increased atmospheric carbon dioxide and the related global climate change (Myneni et al., 1997). 2 APAR (Absorbed Photosynthetically Active Radiation) regulates the intensity of vegetation intensity and therefore NPP (Goward et al., 1985; Landsberg et al., 1996), along with other physical and biological factors, such as canopy structure (Beringer, 1994; Reich et al., 1995), respiration cost (Lavigne and Ryan, 1997; Maier et al., 1998), canopy temperature (Schwarz et al., 1997), and water availability (Will and Teskey, 1997). APAR can be calculated by multiplying FPAR (Fraction of PAR absorbed by vegetation) with incident PAR (Prince, 1991; Prince and Goward, 1995; Running et al., 2004). Algorithms to estimate FPAR using remotely sensed data have been developed, such as algorithms utilizing the NDVI-FPAR correlation (Goward and Huemmrich, 1992; Sellers, 1985), or those based on radiative transfer models (Myneni et al., 1997). With FPAR known, PAR is the most critical parameter for estimating APAR and terrestrial modeling (Dai et al., 2004). Incident PAR is the solar radiation from 400 to 700 nm, and is commonly defined as ? = ? 7.0 4.0 2 )()( ?? dIWmPAR 1-1 where )(?I is spectral irradiance, and ? is the wavelength (Asrar, 1989). Incident PAR can be measured on ground using spectroradiometers, pyranometers, or sensors based on silicon photodiodes (Ross and Sulev, 2000). Although incident PAR is measured at observation stations across the U.S. and other countries (Charlock and 3 Alberta, 1996; Hicks et al., 1996; Ohmura, 1998), a worldwide observation network has not been established to meet the needs for terrestrial modeling at a regional and global scale. Without the sufficient ground measurement at suitable temporal and spatial scale, many methods have been developed to estimate surface incident PAR through other avenues. These PAR estimation approaches can be roughly grouped into three categories: statistical models, physical models, and satellite based methods. 1.2 Current methods for PAR estimation 1.2.1 Statistical Models Statistical models estimate surface incident PAR by establishing empirical relations between PAR and some atmospheric/meteorological factors, without explicitly exploring the interaction between solar radiation and various atmospheric components. Designed to explore the statistical relation between satellite observations and ground irradiance measurements, the Heliosat model was developed to estimate surface solar shortwave irradiance using Meteosat visible band data based on cloud and sky clearness indices (Cano et al., 1986). Subsequent studies improved the original Heliosat model by introducing an atmospheric backscatter term, normalization with clear sky irradiance, and cloud geometrical correction (Beyer et al., 1996). The improved Heliosat model was applied to derive shortwave irradiance for Africa based on Meteosat images (Beyer et al., 1997). Direct and diffuse PAR are estimated by applying the empirical model developed 4 by (Alados-Arboleda et al., 2000) in conjunction with the direct and diffuse shortwave irradiance resulting from the model developed in (Iqbal, 1983). The primary advantage of statistical models lies in their simplicity. These statistical models enable the estimation of PAR without knowledge of the underlying physical mechanism involved in atmosphere radiative transfer processes. On the other hand, statistical models? main drawback is that they are typically site-specific. The statistical relationship established between PAR and other geophysical variables are difficult to be applied to areas other than from where the relationship is derived from, and therefore greatly limits the applicability of statistical models. 1.2.2 Physical models Physical models for PAR estimation explicitly simulate the interaction between solar radiation and the Earth?s atmosphere. The solar radiation/atmosphere interaction processes that are included in physical models generally include Rayleigh scattering, water vapor absorption, ozone absorption, aerosol attenuation, and cloud scattering and absorption (Gu and Smith, 1997). One physical model, CPCR2, was developed by (Gueymard, 1989) to estimate the diffuse and direct component of incident PAR based on extraterrestrial solar irradiance and the effect of various absorption and scattering processes. Because the irradiance estimation output from CPCR2 encompasses the wavelengths between 290 nm and 700 nm that are 5 broader than the spectrum of PAR, an empirical formula has been developed to compute the incident PAR from the CPCR2 output. Another PAR model, developed (Gueymard, 1993) based on similar physical principles, is able to compute diffuse and direct PAR directly because the bandwidth used in the model ranges from 400 to 700 nm. Both models have been applied to estimate PAR under cloudless sky conditions (Alados- Arboleda et al., 2000). In estimating the impacts of smoke and cloud cover on incident PAR, PAR is treated as the top of atmosphere solar flux attenuated by the interactions with atmosphere (Kobayashi and Matsnaga, 2004). Gas absorption, Rayleigh scattering, aerosol scattering and absorption, as well as cloud scattering are included the model. A salient advantage of physical models is that their accuracy improves as better knowledge of radiative atmospheric transfer and other physical processes is gained. Through explicit accounting for physical mechanism, physical models typically are not site-specific and can be applied to extended areas and atmospheric conditions. However, the accuracy of physical models is heavily influenced by the accuracy of input parameters and how well the models represent the real physical processes. The availability of high quality input data is critical for achieving expected modeling accuracy. It is not uncommon that accuracy of a well developed physical model is seriously compromised due to the lack of accurate and reliable inputs. 6 1.2.3 Satellite based PAR estimation methods Algorithms that estimate incident PAR based on satellite remotely sensed data have become increasingly important due to their ability to provide PAR information at sufficient scale and temporal/spatial resolution required by ecosystem modeling. Sensors aboard on geostationary and solar-orbiting satellites have been widely used for PAR estimation. Sensors aboard on geostationary satellites, such as Meteosat and GOES, acquiring data at high temporal resolution, are able to better capture the diurnal variation of PAR. Research has been conducted to estimate PAR (or total shortwave irradiance) based on imagery acquired by these sensors. A model proposed by (Iqbal, 1983) was adopted and modified by Rubio to estimate the hourly direct and diffuse PAR using Meteosat observations (Rubio et al., 2004). The method developed by (Pinker and Laszlo, 1992) to estimate PAR is based on a relationship between atmospheric transmissivity and top of atmosphere reflectivity at the visible region of the spectrum. Atmospheric transmissivity is determined by matching TOA (Top of Atmosphere) reflectivity resulting from model simulations to that observed by the sensor. This method was applied to ISCCP C1 (International Satellite Cloud and Climatology Project) data and generated the first global map of monthly PAR. 7 (Eck and Dye, 1991) developed a method to estimate PAR using the ultraviolet radiance of TOMS (Total Ozone Mapping Spectrometer). Based on the fact that cloud reflectivity is fairly constant and cloud absorption is not significant across the ultraviolet and PAR wavelength, the method parameterizes cloud at PAR region as a simple linear function of TOMS ultraviolet reflectance. Gautier et al. (Gautier et al., 1980) developed a physical model to estimate the total shortwave isolation at surface level using GOES observations. The model calculates surface reflectance based on clear sky TOA radiance and by parameterizing effects of other physical processes, including Rayleigh scattering, water vapor and ozone absorption, cloud absorption, and aerosol attenuation of incoming solar irradiance. This model was later modified by (Gu and Smith, 1997) to estimate PAR surface flux over the BOREAS study area based on GOES observations. To accommodate the climatic and biological difference between the high latitude area and tropical rainforest, this model was further modified to calculate the total shortwave insulation and PAR using GOES data from the Amazon basin (Gu et al., 2004). The major modification involved is focused on the increased cloud absorption factor to better fit the tropical rainforest environment. Visible band data acquired by GOES8 also have been applied to estimate total shortwave flux using a simplified physical model based on the shortwave radiative transfer equations (Ceballos et al., 2004). The impact of cloud cover on incident PAR is approximated by setting minimal and maximal values of TOA spectral reflectance, with 8 Rmin =0.09 representing the clearest sky condition, and Rmax=0.465 representing the most overcast condition. Remotely sensed data acquired by solar orbiting sensors have also been used to estimate PAR. A simplified radiative transfer model was developed by (Van Laake and Sanchez- Azofeifa, 2004) to estimate incident PAR using MODIS atmospheric products. The major atmospheric attenuation processes, including Rayleigh scattering, water vapor absorption, and aerosol scattering, are treated individually in the model. All the input data and parameters are from standard MODIS atmospheric products, including Angstrom turbidity coefficient, atmospheric water content, cloud optical thickness, cloud top pressure, and total atmospheric ozone. This approach was later applied to estimate PAR in Costa Rica (Van Laake and Sanchez-Azofeifa, 2005). Besides estimation of PAR over land, incident PAR over ocean has been estimated based on simplified models using MODIS data (Carder et al., 1999) and SeaWiFS data (Frouin et al., 2000). Because there are many more ground stations measuring shortwave irradiance than PAR (Ross and Sulev, 2000), and more algorithms have been developed to estimate shortwave irradiance using satellite remotely sensed data (Pinker et al., 1995), another commonly used PAR estimation method is to derive PAR from measured or estimated total solar shortwave insolation (Mottus and Sulev, 2001). For PAR estimation, the major advantages of satellite based methods are: first, because these methods draw on remotely sensed data?s advantage in spatial and temporal 9 coverage, they are typically formulated to be applied to the extended area and time span; second, satellite based methods are capable of combining the advantages of both statistical models and physical models. Satellite based methods, especially those hybrid models, typically utilize statistical approaches to supplement physical approaches in cases where physical processes are not well understood or the necessary input data are not available. 1.3 Issues with the existing remote sensing based PAR estimation methods 1.3.1 The conversion factor between SW irradiance and PAR The conversion ratio between total shortwave insolation and PAR varies with time and location as it is affected by many factors including atmospheric pressure, solar elevation, turbidity and precipitable water (Alados et al., 1996). A conversion factor of 0.45 is reported suitable for the Eastern United States during summer season (Pinker et al., 1995). But study also suggests that variability in the ratio of PAR to SW (total shortwave) insolation is not negligible on a global scale (Pinker et al., 1992). Empirical models have been developed to explore the variability of the conversion factor. For instance, a model was developed to estimate the conversion ratio based on Angstrom coefficient for aerosol transmittance and optical air mass (Al-Shooshan, 1997; Alados-Arboleda et al., 2000). (Alados et al., 1996) developed two empirical models for the SW/PAR conversion ratio. The first model estimates the ratio based on four parameters: clearness of sky, brightness 10 of sky, solar elevation angle, and dew point temperature. The second model includes only the first three parameters, excluding dew point temperature. The comparison with ground measurements indicate that the first model performed better (with R 2 of 0.81) than the second model (R 2 = 0.739). In another study, Al-Shooshan utilized neural networks to estimate the conversion ratio, with global irradiance, solar zenith angle and sunshine duration as input data (lopez et al., 2001). 1.3.2 Diffuse and direct PAR The ratio of the two components of the incident PAR, direct and diffuse PAR, has impact on the intensity of vegetation synthesis. Through research on different forest types, (Gu et al., 2003) eported that increased diffuse radiation results in higher light use efficiency and is less likely to cause canopy photosynthesis saturation (Gu et al., 2002). The increased diffuse radiation caused by the Mount Pinatubo eruption also led to enhanced photosynthesis (Gu et al., 2003). Remote sensing based PAR estimation algorithms that sufficiently address the seperation of diffuse and direct PAR are not common. In this study, algorithms centered on extensive MODTRAN4 simulations are developed to estimate PAR using remotely sensed data. Because MODTRAN4 is capable of calculating both direct and diffuse downward spectral flux at surface level separately, it enables the new PAR algorithms to estimate direct and diffuse PAR separately. 1.3.3 Topographic effect on PAR 11 Most algorithms used to generate the satellite based PAR products have not accommodated topographic impact on incident PAR, which could be considerable especially at higher spatial resolution (Pinker and Laszlo, 1992). For example, in validating the VP-RAD (Vapor Radiation) model, Winslow reported that ISCCP-PL PAR product was underestimated in mountainous areas by comparing it with long-term radiation climatologic data (Winslow et al., 2001). Topography has several major effects on surface incident PAR. First, surface elevation affects atmospheric transmittance by changing the actual atmospheric path length (Van Laake and Sanchez-Azofeifa, 2005). Second, the surface slope affects the angle between the surface normal and sun zenith angle, which is proportional to the fraction of light intercepted by the inclined surface (Corripio, 2003). Third, incident PAR value at one location could be affected by the shadow and reflecting radiation of neighboring terrain. Efforts to account for topography have been made (Duguay, 1995; Kumar et al., 1997; Varley et al., 1996). In the procedure used to estimate PAR based on MODIS atmospheric products, a digital elevation model (DEM) was included to account for the effect on the shortened atmospheric path and surface exposure factor caused by elevation and slope (Van Laake and Sanchez-Azofeifa, 2005). To calculate incoming solar radiation at mountainous area using Iqbal?s parametric model, (Corripio, 2003) used DEM to determine whether each pixel was sunlit or in the shade of surrounding pixels. The angle of incidence of the sun on the inclined surfaces is also considered in this study. 12 1.3.4 The impact of snow and ice on PAR estimation In the visible region, snow and ice have spectral characteristics similar to reflective cloud: both snow/ice and cloud appear to be bright in the visible region. This similarity imposes a challenge for estimating PAR using remotely sensed data because of the difficulty in differentiating cloud from snow/ice. For cases in which snow/ice is mis- identified as cloud, PAR would be underestimated. Conversely, when cloud is treated as snow/ice, PAR would be overestimated. Therefore it is vital to differentiate cloud from snow/ice for estimating PAR using remotely sensed data. In this study, in order to account the effect of snow/ice on estimating incident PAR, snow and non-snow pixels will be treated separately. Snow/ice products that are used to this end include MODIS snow cover product (Hall et al., 2001), and snow and ice mapping systems (Ramsay, 1998). More details related with this issue are given in Chapter 4. 1.4 The need for high spatial resolution PAR products Because high-resolution incident PAR over land is not a standard EOS product, MODIS NPP group (MOD17) uses the DAO (Data Assimilation Office) incident PAR products multiplied with MOD15 FPAR (Fractional PAR) to calculate the needed APAR for NPP modeling. While the MOD15 FPAR product is of 1 km resolution, DAO PAR products are of much coarser spatial resolution (2 by 2.5 degree). The resolution disparity and the ensuing disaggregation of DAO PAR impact final NPP product?s accuracy. The availability of PAR product of compatible spatial resolution will help in this regard. 13 To estimate NPP using the GLO-PEM model (GLObal Productivity Estimation Model) (Prince and Goward, 1995), APAR was calculated by multiplying the FPAR product generated from AVHRR NDVI to the TOMS (Total Ozone Monitoring Spectrometer) incident PAR product (Eck and Dye, 1991). The TOMS incident PAR product consists of monthly average estimation, at spatial resolution of 1 by 1 degree from 66 North to 66 South latitude. In addition to the concern for the coarse spatial resolution, the monthly average PAR compromises the accuracy of the modeling because significant variation of PAR may occur over the span of a day. To mitigate this problem, the monthly average PAR was interpolated to generate the 10-days average. 1.5 Objectives of this study The objective of this study is to develop new algorithms to estimate PAR using remotely sensed data acquired by MODIS, AVHRR, and GOES. An important feature of the new PAR algorithms is that the atmospheric information is derived solely from the satellite TOA radiance. Most of the existing satellite based algorithms compute PAR in two major steps: deriving the atmospheric parameters, and computing PAR based on the derived atmospheric parameters (Pinker et al., 1995; Pinker and Laszlo, 1992; Pinker et al., 1992). Of the multiple atmospheric parameters that affect surface flux, the existing algorithms require some, such as water vapor and ozone concentration, to be known, in order to 14 derive the others, such as cloud extinction coefficient and aerosol optical depth. In the new algorithms, the atmospheric parameters are derived solely based on the satellite observations without requiring any external atmospheric knowledge. In addition, for those existing algorithm that derive atmospheric parameters based on the linkage between atmospheric condition and satellite observations at the top of atmosphere, land surface information are also needed due to the fact such linkage is dependent on land surface property. In the new algorithms, land surface reflectance are derived through utilizing the time series of satellite observations, and therefore. Because the new algorithm derived both atmospheric and land surface information solely from satellite observation, they have wider applicability than the existing algorithms in situation where external information for atmospheric or land surface are lacking. This feature makes the new algorithms different from the existing algorithm in that: the existing algorithms all requires some information of land surface or the atmosphere in order to derive the needed atmospheric parameters as the basis for computing The combined use of MODIS, AVHRR, and GOES takes advantage of the temporal resolution and spatial coverage offered by the three sensors to generate multiple years of PAR product covering the whole North America continent as support for ecosystem modeling needs. The core of the new PAR algorithms developed in this study is the derivation of atmospheric and land surface conditions using satellite TOA radiance, which are based on lookup tables constructed using MODTRAN4 simulation results. MOTRAN4 15 simulations are used to establish the relation between satellite TOA radiance and PAR under a variety of atmospheric and land surface conditions. Given the critical role played by MODTRAN4 simulations in the new algorithms, it is inevitable that the uncertainty of MODTRAN4 affects the effectiveness and accuracy of the new PAR algorithms. As it is beyond the scope of this study to analyze the uncertainty of MODTRAN4, the impacts caused by such uncertainty are not analyzed. 1.6 Structure of the dissertation The rest of this dissertation is organized as follows. The second chapter addresses the development of a new PAR algorithm for MODIS data. The third chapter describes using AVHRR visible band data to derive both instantaneous and daily PAR values. In the fourth chapter, a new algorithm to estimate PAR using GOES visible band data is developed, evaluated, and validated. At last, a summary of this work is given in the concluding chapter. 16 Chapter 2 Estimating PAR using MODIS data 2.1 Introduction The amount of solar radiation that reaches the Earth?s surface is determined by several factors, including extraterrestrial solar radiation, interactions between solar radiation and various atmosphere constituents, and interactions between solar radiation and the Earth?s surface. A clear understanding of these interactions and accurate information of the relevant atmospheric and surface properties enables radiative transfer models to compute PAR. The first factor that affects PAR is the variation of extraterrestrial solar radiation. Extraterrestrial solar radiation refers to the solar radiation emanated from the Sun, before reaching the Earth atmosphere. Because extraterrestrial solar radiation is the source of PAR, its variation directly impacts PAR. At present, reliable and accurate databases for extraterrestrial solar spectral flux are readily accessible and have been incorporated in widely used atmospheric radiative transfer models (Berk et al., 1998; Vermote et al., 1997). In terms of estimating PAR, extraterrestrial solar radiation variation is well constrained. The second and most important factor that affects PAR is the Earth?s atmosphere. Here, the atmosphere is used as a broader term that encompasses clouds, aerosols, atmospheric gases (including ozone), water vapor, and all other atmospheric constituents. After emanating from the Sun, solar radiation interacts with various atmospheric constituents as 17 it travels through the Earth atmosphere. These interactions, including absorption and scattering by atmospheric gases, water vapor, aerosol, and cloud, are the primary PAR modulators. Last, the Earth surface also affects PAR, although to a lesser degree than the atmosphere. The Earth surface?s impacts on the magnitude of PAR are primarily through multiple scattering. Multiple scattering refers to the phenomenon in which part of the solar radiation that reaches the Earth?s surface through the atmosphere is reflected back into the atmosphere, interacting with the atmosphere, and eventually reaching the Earth?s surface again. Because of the re-entry into the atmosphere of the reflected solar radiation, the magnitude of multiple scattering?s contribution of PAR is dependent on not only Earth surface conditions but also on atmospheric optical properties. 2.2 The theoretical basis of the new PAR algorithm As mentioned in the last section, atmospheric radiative transfer models are capable of calculating PAR if they are provided with required information for both atmospheric and Earth surface optical properties. A method has been developed to estimated PAR using the standard atmospheric and surface products from MODIS (Van Laake and Sanchez- Azofeifa, 2004). A primary drawback of this approach is its dependence on external input for atmospheric and surface conditions, that limit the applicability of the methods where required external inputs are not available. 18 In this section, the possibility of deriving PAR from at-sensor radiance without external atmospheric information is explored. If proved valid, such an approach will substantially extend the applicability of PAR estimation using remote sensing data. 2.2.1 Top of Atmosphere (TOA) radiance The at-sensor radiance measured by space-born remote sensor is also referred to the top of atmosphere (TOA) radiance. In this study, the at-sensor radiance and TOA radiance are used interchangeably. TOA radiance consists of two components: path radiance, and radiation that is reflected from the Earth surface. The first component, path radiance, refers to the radiance detected by the sensor as the result of backscattering into space by particles and molecules in the atmosphere (Kaufman, 1993). The magnitude of path radiance is solely controlled by atmospheric optical properties, and independent of the Earth surface properties. The second component of TOA radiance is formed by the part of solar radiation that reaches the Earth?s surface and is then reflected into atmosphere and eventually captured by the sensor positioned at the top of the atmosphere. Assuming a flat Lambertian surface, spectral TOA radiance can be expressed by the following equation (Chandrasekhar, 1960; Liang, 2004): 000 00 (,,) (,,) ( )() 1 s s r II E r ? ?? ? ?? ? ? ? ?? ? =+ ? ? 2-1 19 Where 0 (,,)I ? ?? is at-sensor spectral radiance for a given illuminating/viewing geometry: solar zenith angle is 0 ? , sensor zenith angle is ?, ( cos( )? ?= ), and relative azimuth angle is? . The first item on the right side of the equation, 00 (,,)I ? ?? , is path radiance. The second item on the right side of the equation is surface reflected radiance. As described in the last paragraph, the second component of TOA radiance is affected by atmospheric optical properties ( 0 ()? ? : transmittance from the sun to ground, ()? ? : total transmittance from the surface to the sensor, ? : atmospheric spherical albedo), and earth surface properties (surface reflectance s ? ). 2.2.2 Spectral downward flux at surface level Over a Lambertian surface, spectral downward flux at solar zenith angle 0 ? can be calculated by formula 2-2 (Chandrasekhar, 1960; Liang, 2004). 000 00 () () () 1 s s r FF E r 0 ? ? ??? ? =+ ? ? 2-2 Where 0 ? is solar zenith angle, 0 ? = cos( 0 ? ), )( 00 ?F is the downward spectral flux without contributions from the surface, s r is surface reflectance, ? is the spherical albedo of the atmosphere, is the extraterrestrial solar irradiance, and 0 E )( 0 ?? is total transmittance (direct plus diffuse). 20 2.2.3 Exploring the relation between TOA radiance and PAR By definition, PAR equals the integration of spectral downward flux at surface level over the spectrum of 400-700nm. From section 2.2.1 and 2.2.2, it is evident that both PAR and TOA radiance are attenuated extraterrestrial solar radiation as modulated by atmosphere and earth surface. The fact that both PAR and TOA radiance are tied to atmospheric optical properties and surface properties offers the possible utilization of such relationship in order to derive PAR from TOA radiance. As discussed in Chapter 1, the needs for atmospheric optical property related information limit the applicability of existing remote sensing based PAR algorithms. To overcome this problem, this study set out to explore the possibility of deriving PAR from TOA radiance directly without the need for external atmospheric optical properties information. The fundamental basis of the new method is as follows. Because atmospheric optical properties modulate both TOA radiance and spectral downward flux (and therefore PAR), a distinct atmospheric profile (as an combination of all relevant atmospheric properties) will result in a PAR value and a TOA radiance value. In theory, N number of distinct atmospheric profiles will result in N pairs of PAR/TOA radiance values. When the value of N is large enough to exhaust all the possible atmospheric profiles, it will enable us to discern the relation between TOA radiance and PAR. 21 Then the question whether PAR can be derived using only TOA radiance without external atmospheric information hinges on whether the relation between TOA radiance and PAR meets the following criterion: Criterion: Under a given illuminating/viewing geometry and surface property, there are no two distinct atmospheric profiles that yield the same TOA radiance value, but different PAR value. Meeting this criterion means that a TOA radiance value unambiguously points to a PAR value. This feature enables the derivation of PAR solely based on TOA radiance without knowledge of the actual underlying atmospheric profile. There have been studies that derived PAR from TOA radiance based on the empirical/statistical relation established using PAR ground measurements (Alados et al., 2002; Alados et al., 2003). Two major drawbacks of these approaches are: first, the empirical relationship between TOA radiance and PAR is site-specific and therefore not applicable to other areas, second, the empirical relation is based on limited amounts of in situ observations, and there is no vigorous examination whether TOA radiance and PAR meet the criterion described in the preceding paragraph. In this study, a different approach is employed to thoroughly and rigorously examine the relation between TOA radiance and PAR under various atmospheric profiles. One of the major difficulties in establishing TOA radiance/PAR from in situ observations is that 22 only a very limited amount of atmospheric profiles can be included. To avoid this problem, this study chooses to examine the TOA radiance/PAR relation using radiative transfer model simulations, which theoretically can be carried out under any possible atmospheric profile. 2.2.4 MODTRAN4 simulation and PAR/TOA radiance relation In this study, the atmospheric radiative transfer model MODTRAN4 is chosen to conduct the simulations. MODTRAN4 is considered the most comprehensive atmospheric radiative transfer model (Berk et al., 1998). The latest version of MODTRAN4 is capable of simulating a wide range of atmospheric radiative transfer processes. Fully parameterized, MODTRAN4 is capable of calculating the top of atmospheric upwelling radiance for any specific sensor and surface spectral downward flux, that can be integrated to calculate PAR. The downward spectral flux at surface level resulting from MODTRAN4 simulation can be integrated between 400 and 700 nanometers to calculate for PAR value. And the upwelling spectral radiance at sensor altitude (at the top of atmosphere) can be integrated with the sensor?s channel response function to calculate the channel radiance, which is TOA radiance measured by a given channel of a given sensor. 2.2.4.1 Illuminating/viewing geometry 23 To simulate for TOA radiance and surface downward spectral flux, MODTRAN4 needs to be parameterized for atmospheric gases (including ozone), water vapor, aerosols, clouds, and surface reflectance along with line-of-sight geometry. The line-of-sight geometry, often referred to as illuminating/viewing geometry, includes solar zenith angle, sensor zenith angle, and relative azimuth angle. Because both TOA radiance and surface downward spectral flux change with illuminating/viewing angles, the same set of atmospheric profile has to be simulated for different illuminating/viewing geometry. Because it is impossible to exhaust the indefinite number of possible combinations of solar and satellite angles, only representative values of the angles are used to form the illuminating/view geometry for MODTRAN4 simulations. These representative angles are listed in table 2-1. Table 2-1 Solar and satellite viewing angles used to form illuminating/viewing geometry for MODTRAN4 parameterization. Solar zenith angle 0 o , 20 o , 30 o , 40 o , 50 o , 60 o , 70 o , 80 o Satellite viewing angle 0 o , 10 o , 20 o , 30 o , 40 o , 50 o , 60 o , 70 o , 80 o Relative azimuth angle 0 o , 30 o , 60 o , 90 o , 120 o , 150 o , 180 o 2.2.4.2 Atmospheric parameters In MODTRAN4, the parameters that specify the atmospheric optical properties are organized into three categories. The first category is for an atmospheric model, which quantifies the atmospheric gases compositions, water vapor, and ozone concentrations. 24 The second category describes aerosol optical properties. The last category specifies cloud optical properties. With regard to radiative transfer within the visible spectrum, these atmospheric profile parameters exert impacts of various magnitudes. Within the visible spectrum, clouds and aerosols are the dominating modulators, outweighing impacts of water vapor, ozone and other atmospheric gases. Because each atmospheric profile is a combination of the possible values of the whole set of parameters, an exhaustive treatment of all possible atmospheric profiles is too voluminous to be practical. In order to render the numbers of atmospheric profiles manageable while minimizing the error, this study chooses to fully parameterize only cloud and aerosol optical properties using default values for atmospheric gases, water vapor, and ozone. Sensitivity studies through atmospheric radiative transfer model simulation indicates that atmospheric transmittance over the PAR region changes from 1.0 to 0.9603 with ozone concentration varying from 0 to 17.6862 (g/m2), and transmittance changes from 1.0 to 0.9875 when water vapor concentration varies from 0.0 to 6.7258 (g/m2) (Liang et al., 2006). The major effects of cloud on atmospheric transmittance are through the scattering process, while absorption within the PAR region is negligible (Frouin and Pinker, 1995). Thus, the major hurdle in estimating PAR using remotely sensed data is to account for the absorption and scattering caused by aerosol and cloud, which are highly variable in both spatial and temporal domains. 25 In MODTRAN4, cloud is specified by cloud type, and cloud extinction coefficient. Aerosol is specified by aerosol type and atmospheric visibility, a quantity that is convertible to aerosol optical depth. Table 2-2 lists the values used for cloud parameterization, and table 2-3 lists the values used for aerosol parameterization. Table 2-2 Cloud extinction coefficients (km -1 ) at 550 nm, thickness and based heights Altostratus cloud Stratus cloud Stratus/stratocumulus cloud Nimbostratus Cloud 1 1 1 1 5 5 3 5 20 20 10 10 50 50 15 30 Cloud extinction coefficient (km -1 ) 128 56.9 38.7 45 Thickness, km 0.6 0.67 1.34 0.5 Base height, km 2.4 0.33 0.66 0.16 Table 2-3 Aerosol Types and Visibility (km) Rural Aerosol Tropospheric Aerosol Urban Aerosol 5 5 5 10 10 10 20 20 20 50 50 50 Visibility (km) 100 100 100 2.2.4.3 Relation between PAR and MODIS visible band TOA radiance Using the atmospheric profiles described in the previous sections, simulations were conducted using MODTRAN4 and the resulting top of atmospheric spectral radiances were integrated to calculate MODIS band 3 and band 4 radiance, and the resulted surface 26 spectral flux were integrated to compute PAR. Under each illuminating/viewing geometry, each simulation resulted in MODIS TOA radiance and PAR, that are depicted as scatter plots. Figure 2-1 shows the scattering plots between PAR and MODIS band 3 (blue band, 459-479 nm). Solar zenith angle are 0 o , 20 o , 30 o , 40 o , 50 o , 60 o , 70 o , and 80 o , respectively for the 8 figures. Sensor viewing angles are all 0 o , and relative azimuth angles are all 0 o also. Surface reflectance is 0.05. As shown in Figure 2-1, PAR value decreases with the increase of MODIS blue band TOA radiance. Although the decrease of PAR with the increase of MODIS is not strictly monotonic, it is close to monotonic. Examination of the results under all other illuminating/viewing geometry yielded the near-monotonic decrease of PAR with MODIS blue TOA radiance. The near-monotonic decrease of PAR with MODIS visible band TOA radiance meet the criteria stipulated in section 2.2.3, providing the basis for utilizing TOA radiance to derive PAR without explicit knowledge of the underlying atmospheric profile. The development of the new algorithm based on the above described idea is described in detail below. 27 Figure 2-1 Relationship between MODIS blue TOA radiance and PAR from MODTRAN4 simulations under various atmospheric conditions and solar zenith angles. 2.3 MODIS PAR algorithm description 28 In theory, for a given illuminating/viewing geometry, results from a series of MODTRAN4 simulations with surface reflectance increasing at infinite small steps will enable the at-sensor radiance to surface reflectance conversion. But in practice, this method is extremely time-consuming given the large number of possible illuminating/viewing geometry combinations. To overcome this difficulty, an alternative solution is developed. At a flat Lambertian surface, spectral upwelling radiation, as measured by a sensor at top of the atmosphere, can be expressed by equation 2-1, which links at-sensor radiance ),( ,0 ???I with spectral surface reflectance through three variables , s r d F ),( ,00 ???I , and ? . Thus if the values of three variables are known, the relationship between at-sensor radiance and surface reflectance can be determined using equation 2-1 for a given atmosphere condition and illuminating/viewing geometry. The three pairs of surface reflectance/at-sensor radiance values produced by the above described simulations are sufficient to solve for the three variables: , d F ),( ,00 ???I , and ? . In this study, three random values for surface reflectance: 0.0, 0.5, and 0.8, are chosen to simulate the MODIS at-sensor radiance. The resulting values of at-sensor radiance, in conjunction with respective surface reflectance, are used to solve for the three variables. Thus, under a clear atmosphere, MODIS visible band at-sensor radiance can be converted to surface reflectance at any given geometry using equation 2-1 and the values of the three variables. 29 Because MODIS visible band at-sensor radiance increases with the increase of atmospheric turbidity over most natural surfaces except for snow and ice, observations taken under clear atmospheres will have low values in converted surface reflectance. Under most conditions, it is valid to assume that within a reasonably long period of time, there are observations taken under clear atmosphere. Therefore, for a series of observations at a given pixel, sorting for lowest values in converted surface reflectance will lead to those observations taken under clear atmosphere. At a pixel, a time-series of at-sensor radiances are first converted to surface reflectance under assumption of clear atmospheric conditions. The resulting surface reflectance is referred to as nominal surface reflectance thereafter. High value of nominal surface reflectance represents two possible conditions: (1) actual high surface reflectance caused by snow and ice, or (2) heavy cloud structure. Low value of the converted surface reflectance represents observations taken under clear atmospheric conditions. 2.3.1 Generating Look-Up Tables Although the above described method is capable of retrieving both surface reflectance and atmospheric parameters, the extremely time-consuming on-line MODTRAN4 simulations make it impractical for large scale applications. To overcome this obstacle, the look-up tables (LUT) approach is used. 30 To create the look-up tables, MODTRAN4 simulations are conducted for a series of representative configuration of illuminating/viewing geometry and atmospheric conditions. In this study, values of illuminating/viewing angles shown in Table 2-1 are used. Values used to specify cloud and aerosol optical properties are listed in Table 2-2 and 2-3 respectively. MODTRAN4 simulations and subsequent derivation of the three variables of equation 2- 1 created entries for the first look-up table, which links at-sensor radiance to atmospheric surface reflectance through the values of three atmospheric parameters: , d F ),( ,00 ???I , and ? . In this look-up table, each record (row) comprises values of the 7 variables (columns): cloud extinction coefficient/atmospheric visibility, solar zenith angle, sensor viewing angle, relative azimuth angle, , d F ),( ,00 ???I , and ? . A typical look-up table of this kind is shown in Table 2-4. Table 2-4 Look-up table that links TOA radiance to atmospheric profile index. 0 ? : solar zenith angle, ? : sensor zenith angle, ? : relative azimuth angle. ? : atmospheric albedo Atmospheric profile 0 ? ? ? d F ),( ,00 ???I ? ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. In a similar fashion, a second look-up table is created for calculating PAR with known surface reflectance and atmospheric profile. Surface downward spectral flux at a Lambertain surface can be calculated by equation 2-2. (See Section 2.2.2) 31 Treating the product of )( 000 ??? E as one single variable, at a given solar zenith angle and atmospheric condition, Equation 2-2 links surface radiation with surface reflectance through three variables: )( 00 ?F ,? , and )( 000 ??? E . A typical look-up table of this kind is given in Table 2-5. Table 2-5 Look-up table that links PAR with atmospheric profile index 0 ? : solar zenith angle, ? : atmospheric spherical albedo. Atmospheric profile 0 ? )( 00 ?F )( 000 ??? E ? ?. ?. ?. ?. ?. ?. ?. ?. ?. ?. Using the first look-up table in conjunction with the second one, MODIS visible band TOA radiance is linked to corresponding PAR: The first look-up table links a TOA radiance value to atmospheric profile index, and the second look-up table link the atmospheric profile index to a PAR value. This linkage between MODIS TOA radiance and PAR serves as the basis for look-up table searching and PAR derivation, as described in detail in the next section. 2.3.2 Searching look-up table to calculate incident PAR The procedure for calculating incident PAR based on the look-up tables comprises: (1) For a time-series of observations, each at-sensor visible band radiance is first converted 32 to nominal surface reflectance under clear atmospheric condition using entries from the first look-up table and equation 2-1. (2) Cloud pixels are identified based on nominal surface reflectance and reference to the MODIS snow/ice map. (3) Sort the time-series of nominal surface reflectance, and treat the lowest positive 10 percent values of nominal surface reflectance are treated as actual surface reflectance. (4) The values of actual surface reflectance along with their illuminating/viewing angles information are used to compute surface reflectance value at observations taken under non-clear atmosphere. (5) With surface reflectance calculated at each observation, at-sensor radiance at each atmospheric condition, from the clearest to the most heavily cloudy is calculated based on entries of the first look-up table. The resulting series of simulated values of at-sensor radiance are compared with actual at-sensor radiance value to retrieve the atmospheric profile index. (6) With values of both surface reflectance and atmospheric profile index retrieved, the algorithm proceeds to calculate incident PAR by searching the second look- up table. A schematic diagram of the algorithm is presented in Figure 2-2. 2.4 Error budget analysis In the MODIS PAR algorithm, look-up table entries for every representative illuminating/viewing geometry consist of results from MODTRAN4 simulations under a set of atmospheric conditions. As indicated in section 2.2, using more possible atmospheric profiles in generating look-up tables enables a better capture of radiance/PAR relations. But using all possible atmospheric profiles in look-up tables is not practicable in actual algorithm implementations for two major reasons: first, using a 33 large amount of atmospheric profiles will result in large size look-up tables, second, searching of look-up table, as described in section 2.3 requires a monotonic relation between at-sensor radiance and PAR. As Figure 2-1 shows, the relation between at-sensor radiance and PAR is not strictly monotonic. To have the required monotonic relationship between TOA radiance and PAR, this algorithm uses the reduced amount of representative atmospheric profiles which introduces errors through approximation and linear interpolation. Lookup table I (connecting atmospheric condition to TOA radiance) Actual TOA Radiance Derived Surface Reflectance Lookup table II (connecting atmospheric condition to surface radiation Atmospheric Parameter Surface Spectral Flux PAR Figure 2-2 Procedure of the new MODIS PAR algorithm 34 Because using the reduced amount of atmospheric and linear interpolation will increase errors, it is important to evaluate the magnitude of such errors. Figure 2-2 shows the differences in at-sensor radiance/PAR relation between using different amounts of atmospheric conditions. The left panel of Figure 2-3 is generated by using atmospheric conditions presented by only one cloud type (altostratus cloud) and one aerosol type (rural aerosol model), while the right panel is generated by using atmospheric conditions represented by four different cloud types and three aerosol types as listed in Tables 2-2 and 2-3. Figure 2-3 Relationship between MODIS band 3 (blue band) TOA radiance and PAR from MODTRAN4 simulations. The figure in the right panel represents the relation from results of all atmospheric conditions. The figure in the left panel represents the relation from only one cloud type (altostratus) and one aerosol model (rural). In this algorithm, PAR values corresponding to those atmospheric conditions that are not present in look-up tables are linearly interpolated. Such linear interpolation causes error 35 equivalent to the difference between the actual PAR value and the linearly interpolated PAR value. While it is not practical to examine the error caused by linear interpolation under all possible atmospheric conditions, examining the magnitude of errors under representative atmospheric conditions gives a meaningful measurement of the possible error range. Figure 2-4 shows the error caused by the linear interpolation method through comparing the PAR values resulted from MODTRAN4 simulations and PAR values resulted from linear interpolations. The scatter plots in the leaf column of Figure 2-4 show the comparison between PAR values from MODTRAN4 simulations and PAR values from linear interpolations under varying solar zenith angle. The histograms in the middle column show the relative error in percentage introduced by linear interpolations. The right error bars in the right column indicate the relative error in percentage grouped by magnitude of PAR value. A very comprehensive analysis of error caused by linear interpolation was conducted, although only results from some representative conditions are shown here. Based on the error budget analysis, the relative errors caused by linear interpolation under most conditions are under 6%. 2.5 The uncertainty associated with atmosphere model, ozone, and water vapor Precise description of an atmosphere profile requires specifying many parameters, due to the complexity of atmospheric composition. The most commonly used parameters include those for atmospheric gases, water vapor, ozone, cloud and aerosol. 36 37 Figure 2-4 Error caused by the linear interpolation scheme of look-up tables entries generated with selective atmospheric conditions. 38 To best examine the TOA radiance/PAR relation, all the possible combinations of all atmospheric parameters should be used in the MODTRAN4 simulation. But the number of possible combinations is too large to be conducted. In addition, variations of some atmospheric parameters have only insignificant effect on atmospheric transmittance at the visible spectrum. In the visible region, aerosol and cloud have the most significant impacts, and therefore their impacts are thoroughly examined by varying value for the respective parameters. To render the method practical, the PAR algorithm generates look-up table entries from an atmospheric profile with fixed value of water vapor, atmospheric gases composition (including ozone). Therefore the uncertainty caused by such simplification should be examined. 2.5.1 Uncertainty caused by variation in atmospheric ozone concentration The uncertainty in PAR estimation caused by using default value for atmospheric ozone is assessed by MODTRAN4 simulations for atmospheric profiles with varying amount of atmospheric volume ozone. Figure 2-5 shows the range of PAR with atmospheric volume ozone changing from 0 to 700 Dobson Units for a given atmospheric profile. 39 Figure 2-5 Change of PAR with the change of atmospheric volume ozone. Atmosphere condition is represented by mid-latitude summer model with rural aerosol model with visibility of 50 km, altostratus cloud with cloud extinction coefficient of 5 km -1 , water vapor concentration of 2.0g/cm2. Visible region broad band visibility is set at 0.05. Solar zenith angle is 0 o . 2.5.2 Uncertainty caused by variation in atmospheric water vapor While water vapor has significant impact in the infrared region, its impact in the visible region is much smaller. Like that of atmospheric ozone, the effects of atmospheric water vapor on PAR estimation are assessed by MODTRRAN simulations with varying water vapor values. Figure 2-6 shows the change of PAR with water vapor concentration 40 changes from 0 to 8 g/m 2 . Table 2-6 shows the changes of integrated atmospheric transmittance over PAR region with change of water vapor concentration. Figure 2-6 Impacts of atmospheric water vapor on PAR. Atmosphere condition is represented by mid-latitude summer model with rural aerosol model with visibility of 50 km, altostratus cloud with cloud extinction coefficient of 5 km -1 , volume ozone of 350 Dobson Units. Visible region broad band visibility is set at 0.05. Solar zenith angle is 0 o . Table 2-6 Water vapor concentration and the integrated transmittance over the PAR region Concentration (g/m2) Transmittance 0.0 1.0 0.3703 0.9993 2.0162 0.9957 4.0735 0.9919 6.1309 0.9885 6.7258 0.9875 41 2.5.3 Uncertainty caused by variation in atmospheric gases composition. Uncertainty in PAR estimation caused by using default atmospheric gases profiles was examined by MODTRAN4 simulations. In these simulations, the five atmospheric models included in MODTRAN4 were used: tropical, mid-latitude summer, mid-latitude winter, sub-arctic summer, and sub-arctic winter atmospheres. Figure 2-7 shows the range of PAR variation under the five atmosphere models with varying water vapor and ozone concentrations. Figure 2-8 shows the combined impacts of atmospheric gases composition, water vapor, and ozone concentration under varying solar zenith angles. The four figures in the left side column of Figure 2-8 show the ranges of PAR with varying atmospheric profiles represented by volume ozone ranging from 0 to 800 Dobson units, water vapor ranging from 0 to 8 g/m 2 , five atmosphere models used in Figure 2-7. The four figures in the right column show the uncertainty in percentage caused by the combined effect of varying atmospheric gases composition, water vapor, and ozone concentration. Cloud (altostratus) extinction coefficients are 1, 5, 20, and 50 km-1, respectively for figures in the four rows in Figure 2-8. Each figure in Figure 2-8 is grouped by solar zenith angle, which varies from 0 o to 80 o with a 10 o increment. As indicated by Figure 2-8, the uncertainty caused by the combined impacts of atmospheric gases composition, water vapor, and ozone concentration are below 2% under most situations. 42 Figure 2-7 Impact of atmospheric gases composition on PAR. The five error bars presents the five standard atmosphere models included in MODTRAN4: 1: tropical, 2: mid- latitude summer, 3: mid-latitude winter, 4: sub-arctic summer, 5: sub-arctic winter. Under each atmosphere model, water vapor column varies from 0 to 8 g/m 2 , and ozone column varies from 0 to 800 Dobson Units. 2.6 Validation The process of validation was conducted using ground measurements of incident PAR at six FLUXNET sites: Fort Peck (Montana), Lost Creek (Wisconsin), Oak Ridge (Tennessee), Walker Branch (Tennessee), Santarem (Brazil), and Black Hills (South Dakota). For each site, a 3 by 3 window of the MODIS TOA radiance (MOD02) and angular values were extracted from the MODIS data sets ordered from NASA?s Earth Observation System (EOS) gateway. The ground measurements collected at 30-minute or 60-minute intervals were compared with the retrieved values. The measurement values closest to the MODIS data acquisition were used without any interpolation. 43 Figure 2-8 Combined impacts of variations in atmospheric gases composition, water vapor, and ozone concentration on PAR. Water vapor column ranges from 0 to 8 g/m 2 , ozone volume ranges from 0 to 800 Dobson units. Cloud is represented by altostratus cloud with extinction coefficients of 1, 5, 20, and 50 km -1 in the figures of the four rows respectively. 44 The validation results have been reported in a published paper (Liang et al., 2006). In the published paper, a robust regression method called the least trimmed squares regression is used to characterize the fitness between the measured PAR and estimated PAR. In this dissertation, the ordinary least squares regression is used to replace the least trimmed square regression for the following reasons: first, the ordinary least square regression provides a better basis for the cross-comparison with results from other studies; second, to keep the consistency with statistics used in AVHRR (Chapter 3) and GOES (Chapter 4), both of which use the ordinary least squares regression. The Black Hills site in North Dakota is located at 44? 9' 28.8" N, 103? 39'W. The flux tower is situated in an evergreen needle-leaf forest and temperate landscape. The half- hour measurement data in 2002 from days 182-304 were used for the validation, which is shown in Figure 2-9. The solid line in the figure is the 1:1 line, and the dashed line is the fitted line resulting from the ordinary least square regression. The relation statistics are reported in Table 2-7. The Fort Peck, Montana, station is located on the Fort Peck Tribes Reservation, approximately fifteen miles north of Poplar, Montana. The geographic coordinates are 48? 18.473' N and 105? 6.032'W. The three-meter tower was installed in November of 1999 over grassland with a surface elevation of 634m. The comparison of the retrieved and measured half-hour incident PAR in 2002 (the whole year) is shown in Figure 2-10. The overall fit is very good, although both overestimation and underestimation of PAR values were documented. 45 The Lost Creek site (46? 4.9' N, 89? 58.7' W) is located in northern Wisconsin. The elevation is about 480 m above sea level. The biome is mixed forest with canopy height of 1-2 meters. This site was established on Sept. 4, 2000. The measurement data from both 2000 and 2001 were used here and the validation results are shown in Figure 2-11. Figure 2-9 Validation of the estimated instantaneous PAR at Black Hills, ND. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. 46 Figure 2-10 Validation result at Fort Peck, South Dakota. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. Figure 2-11 Validation of the estimated instantaneous PAR at Lost Creek, Wisconsin.. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. 47 The Howland forest (Main Tower) site is located at Howland, Maine (45.20407? N, 68.74020? W). Topographically, the region varies from flat to gently rolling with a maximum elevation change of less than 68 m. The landscape is composed of deciduous evergreen needle forest, boreal/northern hardwood, old coniferous, hemlock, Douglas fir, and evergreen coniferous. Measurement data in 2002 from dates 74-365 were downloaded. The validation results are shown in Figure 2-12. It appears the worst one in all cases where the retrievals are biased. The Walker Branch Watershed site (35? 57' 31.56" N, 84? 17' 14.76" W ), located at Oak Ridge, Tennessee, is characterized by temperate climate, mixed-species, broad-leaved deciduous forest of oak and hickory. The validation of total PAR shown in Figure 2-13 indicates a good agreement between the estimated and measured values. However the validation of both the diffuse and direct PAR result in much more pronounced disagreement between the estimated and measured values (Figure 2-14). One possible reason for the larger disagreement could be the use of the fixed conversion factor between the energy unit (W/m2) and quantum unit (umol/m2/s). MODTRAN4 simulations output both the diffuse and direct flux in the energy unit, to match the quantum unit that are used in the ground measurements, a conversion factor (W/m2=4.6 umol/m2/s) is used. Study reports that the conversion factor varies with atmospheric condition (Dye, 2004). The use of the fixed value of 4.6 for the conversion factor could contribute to the poor validation results for the diffuse and direct PAR. 48 The Walker Branch Watershed site at Oak Ridge is the only site where the direct and diffuse components of the total PAR are validated separately, because this is the only site where ground measurements are available for the diffuse and direct PAR. At the other five sites, only total PAR measurements are available. Figure 2-12 Validation result at Howland Forest, Maine. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. 49 Figure 2-13 Validation result of total PAR at Oak Ridge, Tennessee. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. Figure 2-14 Validation results of direct and diffuse PAR at Oak Ridge, Tennessee. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. 50 Figure 2-15 Validation results at Santarem, Brazil. The solid line is the 1:1 line, and the dashed line is the fitted line from ordinary least square regression. The Tapajos National Forest site (Santa Rem ? Km83) is located in logged evergreen tropical broadleaf forest in Brazil. The geographic coordinates of the tower is (3? 1' 4.905372" S, 54? 58' 17.166324" W). The measured PAR data dates from 1999, but only data in 2001 were used in this study. The comparison is shown in Figure 2-15. The overall fit is also good, although overestimation occurs at large PAR value. Table 2-7 Summaries of regression analysis at six validation sites sites interce pt slope Multiple R2 Relative error (%) N bias RMSE Black Hills 128.60 0.8645 0.9013 6.61 110 76.96 150.91 Fort Peck 193.13 0.8385 0.8804 15.1 353 -38.77 215.41 Lost Creek 66.81 0.9227 0.9247 10.2 441 16.02 171.32 Howland Forest 99.60 0.7394 0.8949 21.9 208 198.94 230.38 Oak Ridge -3.513 1.0202 0.9881 4.1 217 -1.954 15.39 Santarem 182.59 0.8195 0.8661 7.6 167 88.667 156.01 51 2.7 Calculating Daily PAR from MODIS derived instantaneous PAR To meet the NPP and other models? need for daily PAR product, Liang et al (2006) has developed method to compute daily average PAR value using the instantaneous PAR value estimated from MODIS data. The method is based on deriving an empirical relationship between the daily average PAR and MODIS derived instantaneous PAR values through regression analysis. The instantaneous PAR value estimated using both the Aqua and Terra are included in the regression to take advantage of the multiple observations per daytime. Figure 2-16 are reproduced from the published paper (Liang et al., 2006) with permission from the authors. The daily average PAR values are reported in the energy unit (W/m2) and there are no statistics included in the published paper. Figure 2-16 Comparison of the measured daily average PAR and estimated daily average PAR. The solid line is the 1:1 line. 52 2.8 Mapping incident instantaneous PAR over the Greater Washington D.C. area To test the method for regional mapping, PAR over the greater Washington D.C. area was mapped. Figure 2-17 shows the images of the true color composite images of MODIS visible bands acquired on 8 day of 2003 (142, 145, 149, 151, 158, 159, 161), and Figure 2-18 shows the retrieved PAR values. Images are of 800 by 800 pixels. Each pixel?s nominal size is 1 km by 1 km. Visual inspection indicates that the cloud patterns shown in TOA radiance images are well captured on the corresponding instantaneous PAR images. Figure 2-10 True color composite MODIS imagery of greater Washington, D. C. area for 8 days (May 22, May 25, May 29, May 31, June 5, June 7, June 8, and June 10, 2003). 53 Figure 2-11 Retrieved incident PAR maps from the corresponding MODIS images in Figure 2-17 2.9 Summary A new algorithm for deriving PAR using MODIS TOA radiance values without external atmospheric information was described in this chapter. The validity of the new method?s theoretical basis was examined, and associated errors are analyzed. The application of the new algorithm and validation against ground measurement indicate that the new algorithm performs well. Compared with existing satellite based algorithms, the new MODIS PAR algorithm derives atmospheric and land surface information from time-series of TOA radiances 54 with no need for external information for atmosphere. This work has been published as a journal article (Liang et al., 2006), of which I am a co-author. My contribution to this work and publication include: implementation of the MODIS PAR algorithm and validation. Based on the instantaneous PAR estimates using the new MODIS PAR algorithm, Wang (Wang et al., 2007) developed an adjusted sinusoidal model to integrate daily PAR. I am one of the co-authors of the article manuscript submitted to IEEE Transactions on Geophysics and Remote Sensing. 55 Chapter 3 Estimating PAR from AVHRR data 3.1 Introduction In the last chapter, a new algorithm was developed to derive instantaneous PAR values from TOA radiance measured by MODIS visible bands. The two MODIS sensors aboard Terra and Aqua were launched into space on 1999 and 2002, respectively. Therefore, the longest PAR historical record that can be possibly generated from MODIS data covers the most recent 7 years at most. Because one of the primary uses of PAR is for terrestrial ecosystem modeling, which typically requires long-term PAR data for model initiation, parameterization, and validation, this study explores the feasibility of taping into the 20-plus years of continuous observations taken by Advanced Very High Resolution Radiometer (AVHRR). Around this topic, this chapter is organized as follows: first, extending the MODIS-based PAR algorithm to the visible band data of AVHRR; second, examine the error budget in AVHRR PAR algorithm; third, use an adjusted sinusoidal model to scale PAR from instantaneous values to daily average values; And last, validations of both instantaneous and daily average PAR values against ground measurements are presented. 56 3.2 The new AVHRR PAR algorithm Similar to the MODIS PAR algorithm, the prerequisite of the AVHRR PAR algorithm is the existence of the monotonic change of PAR value with the change of the AVHRR visible band at-sensor radiance value. Therefore it is necessary to examine the existence and validity of such relations between PAR and the AVHRR visible band TOA radiance values. Using the same set of atmospheric profiles described in chapter 2 (Table 2-1, 2-2, and 2- 3), the relation between NOAA-14 AVHRR band 1 at-sensor radiance and PAR was examined. Figure 3-1 shows the relation between AVHRR band-1 radiance and PAR for various illuminating/viewing geometries. The figures indicate that PAR decreases with increasing AVHRR band 1 at-sensor radiance in a near-monotonic fashion. Results from all other illuminating/viewing geometries yield similar patterns. Thus, there is a valid basis to develop a new PAR algorithm based on AVHRR visible band data, similar to the MODIS PAR algorithm presented in Chapter 2. 57 Figure 3-1 Relationship between AVHRR-14 band-1 TOA radiance and PAR for a range of various atmospheric profiles and solar zenith angles. 58 3.3 Error budget analysis To derive PAR from AVHRR TOA radiance, a look-up table approach, similar to that used for MODIS, is used. Because of the complexity of including all possible atmospheric profiles in the look-up table, and the near linear relation between TOA radiance and PAR, a simplified method that uses linear interpolation is employed. The simplified linear interpolation method consists of two steps. First, look-up table entries for various illuminating/viewing geometry are generated. Second, linear interpolation is used to derive PAR based on AVHRR band 1 TOA radiance values. In constructing the look-up table for each illuminating/viewing geometry, ideally all atmospheric profiles should be included in order to best capture the varying relation between PAR and TOA radiance. But for the implementation, it is not realistic to include all atmospheric profiles that have been used for the simulations in the exploratory stage. The reasons for not using all atmospheric profiles are: first, the large of number of atmospheric profiles will render the look-up table too large which will make the subsequent look-up table search too slow; second, because the relation between TOA radiance and PAR is not strictly monotonic, a simplified one-to-one relation is needed. The simplification scheme uses results from MODTRAN4 simulations under only one cloud type (altostratus cloud) and one aerosol type (rural aerosol model). Because such simplification introduces error, it is important to quantify the possible errors caused by the linear interpolation under varying circumstances. The left panel of 59 Figure 3-2 shows the relation between TOA radiance and PAR from MODTRAN4 simulations using single cloud/aerosol type. The right panel of Figure 3-2 shows the TOA radiance/PAR relation from simulations using multiple cloud/aerosol types. Figure 3-2 Linear interpolation of PAR based on the relation between PAR and TOA radiance established from simulations using only one cloud type and one aerosol type. The left panel is the TOA radiance-PAR relation based on one cloud/aerosol type simulation, the right panel is based on multiple cloud/aerosol type. 3.4 AVHRR visible band Data The algorithm described in section 3.2 was applied to data acquired by AVHRR aboard the NOAA-14 satellite. There are two forms of AVHRR data: the GAC of 4km spatial 60 resolution, and LAC of 1km spatial resolution. For the purpose of achieving better spatial match between satellite image pixel and point measurements taken on the ground, it would be more advantageous to use 1km LAC data for PAR estimation. But this study chooses to use GAC data in the new AVHRR algorithm, for two reasons: first, because LAC data only cover a set of localized areas, it is difficult to find collocated validation sites; second, the less consistent coverage of LAC data does not meet the requirement of frequent repeat observations by the new AVHRR PAR algorithm. 3.5 Scaling instantaneous PAR to daily average PAR As discussed in the beginning of this chapter, needs for daily PAR values (versus instantaneous PAR values) necessitate the instantaneous-to-daily scaling. Scaling methods have been developed for scaling, but none of them have universal applicability and absolute superior accuracy. In this study, the adjusted sinusoidal model developed for MODIS (Wang et al., 2007) is used. The adjusted sinusoidal model infers PAR distribution function over the course of a daytime by assuming the atmospheric condition changes between individual instantaneous observations are smooth. The instantaneous PAR distribution functions over time are expressed by equation 3-1 for situations when there is only one observation over a daytime, or by equation 3-2 for situations when there are two or more observations. 61 sunrisesunset sunriseoverpass sunrisesunset sunrise overpass TT Tt TT Tt TInstPARtInstPAR ? ? ? ? = ? ? )( sin )( sin )()( 3-1 where is the instantaneous PAR at time t , is the sunrise time, is the sunset time, is the satellite overpass time (time of instantaneous PAR value derived from satellite observation). )(tInstPAR sunrise T sunset T overpass T )()()( 21 12 1 12 2 tInstPAR TT Tt tInstPAR tT tT tInstPAR TT ? ? + ? ? = 3-2 where and are the time of the two satellite overpass, respectively. and are the instantaneous PAR value at time t as calculated using the sinusoidal model based on the first satellite overpass time and second overpass time, respectively. 1 T 2 T )( 1 tInstPAR T )( 2 tInstPAR T With instantaneous PAR distribution modeled, the daily accumulative PAR can be calculated by integrating instantaneous PAR value over a daytime (Equation 3-3). ? ? ? ? = ++ + + + ? ? + ? ? += 1 1 11 1 1 11 1 )()( ) ) ()( T T N i T T T Ii i T ii i T sunrise I I iI dttInstPAR TT Tt tInstPAR TT tT dttInstPARDailyPAR 3-3 ? + sunset n N T T T dttInstPAR )( Where is the daily accumulative PAR. DailyPAR 62 Applying the adjusted sinusoidal model to AVHRR data consists of two steps. First instantaneous PAR values are estimated from individual AVHRR band-1 radiance. Second, all instantaneous PAR values for a given day are used to compute the daily PAR value using the adjusted sinusoidal model. The number of AVHRR satellite daytime overpasses varies with location, although typically there are more frequent overpasses at higher latitudes than mid and lower latitudes. Figure 3-4 illustrates the number of daytime overpasses in North America by AVHRR onboard NOAA-14, during day 143 of 1996. Figure 3-3 Number of satellite overpasses by NOAA-14 AVHRR over the North America continent for day 143 in 1996. 3.6 Validation The newly developed method is implemented and resulting PAR estimations panel were examined for accuracy. Figure 3-5 shows the estimated instantaneous PAR covering the 63 continental U.S., along with corresponding at-sensor radiance of AVHRR-14 band-1 image. It can be observed from the imageries that the PAR estimation method successfully captured the cloud?s impact on PAR. To quantitatively assess the accuracy of the algorithm, the derived PAR values are validated against ground measurements at three FLUXNET sites: Metolius in Oregon, Park Falls in Tennessee, and Walker Branch in Wisconsin. The estimated and measured PAR values for all three sites are from January to March, 1996. The flux tower at Metolius in Oregon is located at a site composed of mixed-aged trees with open spaces. The elevation of the site is 916 m. The scatter plot, depicting the comparison between the measured and estimated instantaneous PAR values along with linear regression line, is in Figure 3-6. The estimated and retrieved PAR values agree well although both underestimations and overestimations are documented. The second site, Park Falls site in Wisconsin, is managed for timber extraction and biodiversity conservation. Figure 3-7 compares the estimated and measured instantaneous and daily PAR values at this site. At the Walker Branch in Tennessee, the dominant trees are oak, maple, and hickory. Figure 3-8 shows the comparisons between estimated and measured PAR values. 64 Table 3-1 and Table 3-2 summarize the validation results at the three sites. The statistics at all three sites indicates that the daily PAR value is more accurate than their instantaneous counterparts. Figure 3-4 AVHRR band-1 data and estimated instantaneous PAR value covering the continental U.S.. The image in the upper panel is TOA radiance at 20:05 GTM, day 225 in 2004. The image in the lower panel is the estimated instantaneous PAR for the same time. The white pixels on the TOA radiance image are clouds. The unit of TOA radiance is Wm -2 sr -1 s -1 . The unit of instantaneous PAR is umol/m 2 /s. 65 Figure 3-5 Validations of PAR derived using AVHRR band 1 data from NOAA-14 at Walker Branch, Tennessee. The left graph is a comparison between retrieved (X-axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis). Figure 3-6 Validations of PAR derived using AVHRR band 1 data from NOAA-14 at Metolius, Oregon. The left graph is a comparison between retrieved (X-axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis). 66 Figure 3-7 Validations of PAR derived using AVHRR band-1 data from NOAA-14 at Park Falls, Wisconsin. The left graph is a comparison between retrieved (X-axis) and measured (Y-axis) instantaneous PAR. The right graph is a comparison between retrieved (X-axis) and measured daily PAR (Y-axis). Table 3-1 Validation results of instantaneous PAR at three FLUXNET tower sites. Metolius Walker Branch Park Falls R 2 0.8851 0.9083 0.8511 Slope 0.9542 0.9370 0.8663 Intercept 65.21 -35.588 46.04 Bias -16.70 105.82 65.5 RMSE 197.40 163.52 199.3 Relative Error 22.97 32.57 40.80 RMSE/Average 18.35 16.205 25.91 N 558 455 516 67 Table 3-2 Validation results of daily integrated PAR at three FLUXNET tower sites. Metolius Walker Branch Park Falls R 2 0.9675 0.9275 0.9299 Slope 0.8383 0.8212 0.8290 Intercept 559.86 616.51 660.08 Bias 896.48 1030.19 567.02 RMSE 1310.28 1468.06 1419.87 Relative Error 18.19 22.41 32.46 RMSE/Average 15.67 17.94 21.47 N 341 333 322 3.7. Summary and discussion Deriving PAR from historical AVHRR data sets has the potential to provide terrestrial ecosystem modeling community with the much needed long time series of input parameter. To realize this potential, a new method for mapping PAR from AVHRR data was developed in this chapter. This new method was applied to AVHRR data of 1999 to generate both instantaneous and daily PAR estimates. The results of the initial validation are promising, though more thorough and systematic validations are needed. Based on the result of the initial validation, further refinements of the AVHRR PAR algorithm has been envisioned: (1) to improve the lookup tables, which should host more entries with finer intervals for both solar and sensor zenith angels; (2) to fuse the AVHRR derived instantaneous PAR value with those from other sensors, including MODIS, GOES, and SeaWiFS, to improve the integration of daily PAR. 68 Chapter 4 Estimating PAR from GOES data 4.1 Introduction In Chapter 3, scaling of instantaneous PAR to daily average value was discussed. One of the key issues of improving the daily PAR scaling accuracy is the number of instantaneous PAR per day. With everything else equal, the more instantaneous PAR values used for scaling, the more accurate the resulting daily PAR value will be. This improved accuracy is primarily because the more instantaneous PAR values distributed over a day enable a better capturing of the diurnal variation of PAR. Both MODIS and AVHRR are only capable of acquiring low numbers of day time observations per day because of they are aboard polar-orbiting satellites. By contrast, sensors aboard geostationary satellites such as GOES and METEOSAT, have much higher temporal resolution and therefore are better equipped to tackle the instantaneous- to-daily PAR scaling. Geostationary satellites have long been used to estimate surface solar insolation and PAR. Gautier et al. (Gautier et al., 1980) developed a simple physical model to estimate incident solar radiation using the first generation GOES visible data, which was later improved by including ozone absorption and introducing an empirical correction for sub- pixel clouds (Diak and Gautier, 1983). GOES observations were used to retrieve total solar radiation and PAR (Gu and Smith, 1997), and surface net radiation at the Boreal Ecosystem-Atmosphere Study (BOREAS) (Gu et al., 1999) and Amazon basin (Gu et al., 69 2004). GOES-8, 10, and 12 were used to calculate surface insolation and the results were validated using insolation data taken at the U.S. Climate Reference Network (USCRN) (Oktin et al., 2005). GOES visible band data were used to estimate incoming solar radiation over northwest Mexico (Stewart et al., 1999) and Brazil (Ceballos et al., 2004). Diner et al. (Diner et al., 1998) estimated daily solar radiation from GOES-8 visible imagery in an effort to apply satellite data in agricultural management. Incoming solar radiation estimated from GOES-8 imagery was used as part of input for modeling evapotranspiration from wetlands at north central Florida (Jacobs et al., 2002). To retrieve PAR, all existing algorithms require some external knowledge of either surface or atmospheric conditions. This is due to the fact that at-sensor radiance includes a mixed contribution from both surface and atmosphere, and the PAR retrieval requires separation of the two. The dependence of these PAR algorithms on external information limites their applicability, especially when the required information is lacking or erroneous. To address this problem, the new PAR algorithm developed for MODIS and AVHRR in the previous chapters was modified to be applied to data acquired by GOES imager?s visible channel. Similar to the PAR algorithm for MODIS and AVHRR, the GOES PAR algorithm is based on the simultaneous derivation of land surface and atmospheric parameters without external atmospheric information. Compared with MODIS PAR algorithm, there are two major improvements for the GOES PAR algorithm. The first improvement is the utilization of spatial relation for excluding the adverse effects of 70 cloud shadow on surface reflectance retrieval. The second improvement is the incorporation of a Bi-directional Reflectance Distribution Function (BRDF) for characterizing surface reflectance. In addition, the impact of topography on surface incident PAR also were assessed in order to provide the end-user of PAR product with topographically corrected PAR data. The rest of this chapter is organized as follows. Section 4.2 gives a detailed description of GOES PAR algorithm. Section 4.3 evaluates the performance of the new algorithm through validation against ground measurements. Section 4.4 examines the impacts of topographic impacts on PAR. Section 4.5 discusses some problems and possible further improvements. 4.2 Algorithm description The central theme of the GOES PAR algorithm is deriving PAR using GOES visible band data through the look-up table approach. The details of the algorithm are addressed in the following sections. 4.2.1 GOES-12 visible band at-sensor radiance GOES visible band data used in this algorithm were acquired by GOES-12. GOES-12 which was launched into orbit (W 75, over the equator) in 2001 to replace GOES-8 as the operational GOES East Satellite. An imager and a sounder are aboard GOES-12. GOES- 71 12 imager?s visible channel (0.55-0.72 um) has a nominal spatial resolution of 1km by 1km at nadir. It is referred to as nominal resolution in the sense that the east-west spatial resolution is actually 0.57 km at nadir due to over-sampling. The visible band data used in this study are taken by GOES-12 imager at 30-minute interval. Like its predecessors, GOES-12 imager?s visible band does not have an on-board calibration device. Sensor responsivity degradation and post-launch vicarious calibration have been a constant issue for GOES visible band observations (Weinreb et al., 1997; Knapp and Haar, 2000). In this study, the 16 bits integer count of GOES-12 visible band data are converted to nominal albedo values. Nominal albedo values were calculated by calibrated radiance and solar radiance. It is called nominal albedo because it is neither corrected for the solar zenith angle, nor the variation of Sun-Earth distance. Because visible band radiance is a required input for the GOES PAR algorithm, the visible band radiance is retrieved from converted nominal albedo values by applying the look-up table provided by NOAA. Because the visible band data used in this work were acquired in mid 2004, more than two years after the launching of the GOES-12 satellite, the on-orbit sensor degradation needs to be taken into account. Vicarious calibration conducted by NOAA using star- viewing data put GOES-12 imager visible channel response degradation rate at 5.68% per year during the period from January 22, 2003 to January 23, 2005. Based on this calibration result, the visible band radiance values were adjusted by using the degradation rate and number of days since the launching. 72 4.2.2 Examine the relation between PAR and GOES visible band TOA radiance Because the GOES PAR algorithm is built on retrieving PAR from TOA radiance with no external atmospheric information, its validity is based on the existence of PAR value?s monotonic change with the change of GOES visible band TOA radiance, as discussed in Chapter 2. Therefore, the relation between GOES visible band TOA radiance and PAR under various atmospheric conditions has to be examined first. Using the same set of atmospheric profiles detailed in Chapter 2, MODTRAN4 simulations were conducted. After the simulations are completed, the resulting top-of- atmospheric spectral radiances are integrated to calculate GOES-12 imager?s band-1 radiance, and the resulting surface spectral flux was integrated to PAR. Then, simulated resulting GOES-12 band-1 TOA radiance and PAR were examined. Figure 4-1 shows the scatter plots comparing PAR and GOES band-1 at-sensor radiance. Solar zenith angle are 0 o , 20 o , 30 o , 40 o , 50 o , 60 o , 70 o , and 80 o , respectively for the 8 figures. Both sensor viewing angles and relative azimuth angles are all 0 o . 73 Figure 4-1 Relationship between GOES12 band-1 TOA radiance and PAR at a range of different atmospheric profiles and solar zenith angles. Surface reflectance is 0.05, sensor zenith angle and relative azimuth angle are all set to 0 o . 74 As shown in Figure 4-1, PAR values decrease with the increasing GOES-12 band-1 TOA radiance. Although the decrease of PAR with increase of GOES TOA radiance is not strictly monotonic, it is close to monotonic. Examination of the results under all other illuminating/viewing geometries yield the near-monotonic decrease of PAR in respect to the increase in GOES-12 band-1 TOA radiance. The near-monotonic decrease of PAR with change of TOA radiance meets the criteria stipulated in section 2.2.3, providing the basis for utilizing TOA radiance to derive PAR without explicit knowledge of underlying atmospheric profile. 4.2.3 Excluding cloud-shadowed observations Cloud shadow detection algorithms have been developed (Simpson et al., 2000) for cloud screening and other purposes. In this study, cloud shadowed observations were excluded based on the fact that cloud shadowed pixels are adjacent in space to cloud pixels. The distance from a cloud pixel to its shadow on the ground is a function of cloud height, solar zenith angle, and solar azimuth angle (Simpson and Stitt, 1998). The knowledge of these angles allows the calculation of the exact position of cloud pixel?s shadow on the ground surface. Due to the complexity of deriving cloud height, a simpler method was used to exclude cloud-shadowed observations in this study. In this study, the highest solar zenith angle of observations used for deriving surface reflectance is 75 o . Visual interpretation of cloud edge and its shadow edge indicated that most shadows are with 10 pixels (10 km at nadir) away from their corresponding cloud structures. 75 Thus, a two-step cloud shadow exclusion approach was used: First, at a given pixel?s time series, high values of nominal surface reflectance were identified and the NOAA snow/ice cover map is used to exclude snow/ice pixels. Within the time series, observations with high nominal surface reflectance but not identified as snow/ice were treated as observation under cloud. Second, for each identified cloud pixel, every pixel within a 10 pixels of radius was labeled as possible cloud-shadowed pixel, and was excluded from subsequent sorting for clear observations. 4.2.4 Fitting BRDF model to calculate surface reflectance The sorting of nominal surface reflectance will result in a set of lowest positive value from the time series of a given pixel. These lowest 10 percent of positive values of nominal surface reflectance, excluding those that fall within a 10 pixel radius of identified cloud pixels, are regarded as actual surface reflectance. The actual surface reflectance are used to fit a modified three-parameter linear Ross-Li model to characterize the surface bi-direction feature, which in turn will be used to calculate surface reflectance value for the observations taken under non-clear atmosphere. The modified three-parameter linear Ross?Li model can be written as: ),,(),,(),,( ????????? vsgeogeovsvolvolisovss KfKffr ?+?+= 4-1 where: s ? is solar zenith angle, v ? is sensor zenith angle, ? relative azimuth angle, is isotropic coefficient, is volumetric coefficient, is geometric coefficient, iso f vol f geo f 76 ),,( ??? vsvol K is volumetric kernel function, and ),,( ??? vsgeo K is geometric Kernel function. After the fitting of the BRDF model, the procedure then compares the values of the standard deviation ( fit ? ) of the fit against a predefined threshold value max ? ( 025.0 max =? in absolute value or 25.0 max =? in relative value), which represents the maximum value of the standard deviation of the regression that is considered acceptable for successful modeling fitting. When the condition max ?? ? fit is fulfilled, the process exits the procedure, otherwise the surface reflectance value exhibiting the largest absolute departure with respect to the model prediction is eliminated, and the fitting process starts anew. This iteration procedure is pursued until an acceptable fit is obtained or the number of surface data points remaining in the time series becomes too low to ensure a reliable retrieval of the geophysical parameters. After the parameters of the BRDF are obtained, values of surface reflectance at observations under non-clear atmosphere are calculated using the illuminating/viewing angles along with the parameters of the fitted BRDF model. The use of the BRDF model has effect on the surface reflectance derived from GOES visible band TOA radiance. As shown in Figure 4-2, the surface reflectance resulting from the BRDF model simulation (represented by black dots) are appreciably different from the surface reflectance that converted from the TOA radiance without BRDF model fitting. In general, the BRDF model removed the extremely large and small reflectance 77 value, which would have significant impacts on resulting surface reflectance from linear interpolation method. Figure 4-2 An examples of the TOA radiance converted surface reflectance (triangles) and the BRDF simulated surface reflectance (crosses) during day 191-198, 2004 at Lost Creek, Wisconsin. The BRDF model used in GOES PAR algorithm is a modified three-parameter linear Ross?Li model, which is used with an assumption that the land surface bi-directional reflecting property keeps stable during the period time used for the BRDF model fitting. The BRDF model requires a sufficient amount of individual satellite observations with different viewing/illuminating geometry for model fitting. The 30 minute observation interval of the GOES provides approximately 20 daytime observations a day, while MODIS and AVHRR typically only have no more than 4 observations at the mid and low 78 latitude. To meet the required amount of observations, it is inevitable to use much longer time series for MODIS and AVHRR, which is in violation of the BRDF model?s underlying assumption: the land surface bi-directional reflecting property keep stable. Therefore the BRDF model is not used in the MODIS and AVHRR algorithms. 4.2.5 Error budget analysis In this GOES PAR algorithm, look-up table entries for every representative illuminating/viewing geometry consist of results from MODTRAN4 simulations under a set of atmospheric conditions. To render the look-up table of manageable size and to meet the strict monotonic change of PAR in regard to GOES-12 band-1 TOA radiance, the look-up tables only contain the entries generated from MODTRAN4 simulation with atmospheric profiles specified by one cloud type (altostratus cloud) and one aerosol type (rural aerosol model). The error generated by using reduced amount of atmospheric profile and linear interpolation was assessed similar to that of MODIS algorithm conducted in chapter 2. The analysis results indicate that the error caused by using simplified atmospheric profiles varies with both solar zenith angle and magnitude of PAR, and the majority of the relative errors are less than 5%. 4.3 Validation 79 To assess the accuracy of the algorithm, the derived PAR values are validated against ground measurements at four FLUXNET sites: Canaan Valley in West Virginia, Lost Creek and Willow Creek in Wisconsin, and Metolius in Oregon. The locations of the four validation sites are summarized in Table 2. To mitigate the adverse effect of the spatial mis-match between satellite observations and ground sites, a 3 by 3 pixel window centered on the pixel that matches the ground site?s latitude and longitude is used. The arithmetic mean of the derived PAR values of the 9 pixels within the window is compared with the PAR value measured at the corresponding ground site. To temporally match the observed and derived PAR values, linear interpolations are used to derive the ground measured PAR value at the time that matches the satellite observation. For example, if the satellite observation is taken at time , which lies between two ground measurements taken at time and , a linear interpolation is used to derive the ground measured PAR value at time based on the measured values at time and . At all four sites, the estimated and measured PAR values are from days 191 to 198, 2004. sat T 1_mea T 2_mea T sat T 1_mea T 2_mea T 80 Figure 4-3 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at Canaan Valley, West Virginia. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. The first validation site, the Cannan Valley site in West Virginia, is located in temperate grassland. The elevation of the site is 1000 meters, and the tower is four meters high. The scatter plot depicting the comparison between the measured and estimated instantaneous PAR values along with fitted linear regression is shown in Figure 4-3. There are some underestimations and overestimations, though most points are within the vicinity of the one-to-one line. The comparison has a RMSE of 190.223 umol/m 2 /s, and bias of 53.639 umol/m 2 /s. 81 Figure 4-4 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at the Lost Creek site in Wisconsin. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. The second site at Lost Creek Wisconsin is situated in an alder-willow wetland. The flux tower is about nine meters above the ground. The underlying terrain is 480 meters above sea level. Figure 4-4 compares the estimated and measured instantaneous PAR values at this site. The comparison resulted in a RMSE of 141.43 umol/m 2 /s, and bias of 11.09 umol/m 2 /s. The liner regression between estimated and measured PAR resulted in a R 2 Value of 0.9227 82 Figure 4-5 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at the Willow Creek site, Wisconsin. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. At the third site, Willow Creek in Wisconsin, the dominant vegetation is alder willow. Figure 4-5 shows the comparisons between estimated and measured PAR value, which has a RMSE of 131.44 umol/m 2 /s, and a bias of 16.56 umol/m 2 /s. Validation at the last site, Metolius in Oregon, resulted in the highest bias of the four validation sites (Figure 4-6). The bias value of -101.54 points to a systematic underestimation of PAR, which is also indicated by the high R 2 value resulted from the linear regression between estimated and measured PAR value. 83 Figure 4-6 Scatter plot between estimated instantaneous PAR (X axis) and measured instantaneous PAR (Y axis) at Metolius site, Oregon. The solid line is the 1:1 line, and the dashed line is the fitted linear regression line. Table 4-1 Summary of site locations and validation statistics at the four validation sites. Metolius Willow Creek Lost Creek Canaan Valley R 2 0.9434 0.9302 0.9227 0.8255 Slope 1.01 1.022 1.0302 0.9173 Intercept 89.39 -39.114 -42.081 16.08 Bias -101.54 16.56 11.09 53.639 RMSE 130.71 131.440 141.443 190.223 Relative Error (%) 14.59 13.97 14.16 22.22 RMSE/Average (%) 9.52 13.01 13.92 24.09 N 456 445 438 436 4.4 Correcting topographic impacts on PAR 84 PAR value estimated using the method described in the previous sections does not take into account the effect of land surface topography. Complex land surface topography has impact on the amount of PAR that is actually available for vegetation?s photosynthesis (Winslow et al., 2001). Many studies have be done to assess and correct topographic effects on surface solar radiation (Corripio, 2003; Duguay, 1995; Kumar et al., 1997; Van Laake and Sanchez-Azofeifa, 2005; Varley et al., 1996). Topographic correciton for surface solar radiation modeled using both polar-orbiting (Dubayah 1992) and geostationary satellites (Dubayah 1997) been reported. In addition, topographic impact on incident PAR is closely related to spatial resolution. For instance, reflected PAR flux from neighboring pixels at finer resolution becomes within-pixel phenomenon as resolution grows coarser. In this section, a set of methods adopted from previous studies are used to evaluate topographic impact on PAR at 1 km nominal resolution in which PAR is retrieved from GOES-12 visible band data. A brief description of the topographic correction method is as follows. 4.4.1 Digital Elevation Dataset When making topographic correction, each pixel is treated as a block of terrain with uniform topographic and reflecting characteristics. The underlying topography is provided by the USGS GTOPO30 DEM dataset. The GTOPO30 has a global coverage with 30 arc second in spatial resolution. To match the grid of GOES PAR, the USGS GTOPO30 DEM dataset covering the study area is re-projected to the same projection 85 system as that of GOES PAR with 1km nominal spatial resolution. Figure 4-7 shows the USGS GTOPO30 DEM covering the continental U.S. and part of southern Canada. 0 1250 2500 3750 5000m Figure 4-7 USGS GTOPO30 DEM projected to Lambert Azimuthal Equal Area projection, with center of projection at 90W/30N (Longitude/Latitude). The study area covers the continental U.S. and part of southern Canada. The three parts are considered in the GOES PAR 1km topographic correction include: correction for direct component of PAR based on angular effects; correction of diffuse component by sky view factor, and correction for the reflected direct and diffuse PAR flux from neighboring terrains. 86 4.4.2 Angular effects on direct PAR Direct PAR received on a given terrain is determined by the angle ( I ) between the solar illumination direction and terrain surface normal (Wang et al., 2005). Angle I is expressed as formula 4-2. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = )cos( )sin()sin( )sin()sin( )cos( )sin()sin( )cos()sin( )cos( S TS TS Z AZ AZ I s s s ss ss 4-2 s Z : Solar zenith angle : Solar azimuth angle s A S : Terrain slope : Terrain aspect s T For an arbitrary pixel X, the topographical impact on direct radiation is calculated as an impact factor ranging from 0 to 1, with 0 representing no direct PAR is received due to terrain?s effect, and 1 representing 100 percent of the available direct PAR is received. In addition, effects of neighboring terrain?s shadows are corrected by a two-step approach: first, calculating the shadow of each pixel cast at the given solar zenith and azimuth angle; and second, testing whether a given pixel is in any other pixel?s shadow. If a pixel is in any other?s pixel?s shadow, it will be assigned a shadow factor of 0, otherwise its shadow factor will be 1. The total correction factor for angular effect on direct component of PAR is the product of the impact factor and shadow factor. 87 Figure 4-8 shows the impact of topography on direct PAR at 19:05 GMT on May 04, 2005. Within the continental U.S., the Rocky Mountain and the Appalachian Mountain range areas reveal the most significant topographical effect on direct PAR. 0% 25% 50% 75% 100 % Figure 4-8 Effect of topography on direct PAR based on the angle between solar illumination direction and the normal to the terrain slope. Dark colors represent pixels where low percentage of incident direct PAR is received. Bright colors represent pixels where high percentage of incident PAR is received. Note that the topographic effect on direct PAR is dependent on time of day and Julian day of year. This image is of 19:05 GMT, May 04, 2004. 88 Figure 4-9 Histogram of values of sky view factor covering the same area as Figures 4-8. The sky view factor ranges from 0.3 to 1.0 with mean value of 0.98 (X axis) and standard deviation of 0.019. The Y axis is the percentage value. 4.4.3 Sky view factor for diffuse PAR In this study, the diffuse PAR is treated as isotropically distributed. Therefore, the topographic impact on diffuse PAR is determined only by the terrain topography. The topographic effect on diffuse PAR is corrected by a sky view factor, which is defined as the ratio of diffuse irradiance at a point to that on an unobstructed horizontal surface (Dubayah and Rich, 1995). The sky view factor is calculated using the angle to the local horizon in for 8 sub-directions of the pixel in question. At each sub-direction, the sub-sky view factor is calculated by the elevation of a blocking pixel and the distance from the blocking pixel to pixel X. In order to find the maximum 89 blocking angle, in each direction, pixels within the distance of the searching radius are examined, and the largest blocking angle is identified. Sky view factor of pixel X then is calculated as the average of 8 sub-sky view factors of the 8 sub-directions. The view factor ranges from 0 to 1, with 0 representing a completely obstructed view, and 1 representing a unobstructed view. Figure 4-9 shows the histogram of view factor values covering the study area. 4.4.4 Reflected direct and diffuse flux from neighboring terrains Previous studies have indicated that there is little justification to use sophisticated methods to correct the reflected radiation from neighboring terrains as the resultant improvement on accuracy is insignificant (Duguay, 1995). In this study a method developed by Dubayah (1992) is used to account for the reflected radiation from neighboring pixels. The reflected flux from neighboring terrains can be expressed as: )( difdirmeantr EECE += ? 4-3 where is the terrain configuration factor: t C skyt S C ?? + = ] 2 cos1 [ 4-4 mean ? is the average reflectance of the 8 neighboring pixels, S is the solar vector, is the direct component of PAR, is the diffuse component of PAR, and dir E dif E sky ? is the sky view factor as described in section 4.3. 90 4.4.5 The results of the topographically corrected PAR The value of the topographically corrected PAR is the sum of the three components described in the previous three sections, and can be expressed as: rskydifdirdirtopo EEfEE +?+?= ? 4-5 where is the topographically corrected PAR, and is the factor of topographic impact on direct PAR. Other symbols are the same as in equations 4-3 and 4-4. topo E dir f Validation of topographically corrected PAR is difficult because of lack of ground measurements. This is because pyranometers are typically horizontally oriented regardless of underlying slope. This study does not carry out validation the topographically corrected PAR. Figure 4-10 shows the topographically corrected instantaneous PAR covering the continental U.S. and part of southern Canada, along with corresponding at-sensor radiance GOES-12 visible band images. It can be observed from the imageries that topographic impact on PAR is much more pronounced in mountainous region than in relatively flat regions. It is also evident that topographic impact is much less significant in areas with heavy cloud cover than areas under clear atmosphere. 91 4.5 Discussion A new algorithm for estimating PAR using GOES visible band data is developed in this paper. Different from the existing GOES PAR algorithms, this new algorithm derives both surface reflectance and atmospheric parameters simultaneously from at-sensor radiance values, without the requirement for external knowledge of atmospheric parameters. Validation against ground measurements indicates that this new algorithm is capable of reaching reasonably high accuracy. Compared with the similar algorithm that was developed to estimate PAR from MODIS visible band data, this GOES algorithm is improved in two important aspects: first, the GOES algorithm counts for the surface bi- direction reflectance using a semi-empirical BRDF model while the MODIS PAR is based on assumption of Lambertian surface. Second, this new GOES PAR algorithm utilizes spatial relation to exclude possible cloud shadowed pixels, while the MODIS method uses only temporal relation for cloud shadow exclusion. However, further improvement is necessary. The topographic correction has not been validated due to the lack of ground measurements. The error associated with the DEM used for topographic correction may introduce extra uncertainty into the topographically corrected PAR. Future studies should consider using topography data with greater accuracy to improve the accuracy of topographic correction. For instance, DEM generated by the shuttle radar topography mission (STRM) at 30m spatial resolution, with global coverage between 60 north and 56 south latitude, is an alternative to the GTOPO30 DEM. 92 93 0 20001000 0 400200 Figure 4-10 GOES visible band data and topographically corrected PAR covering the study area. The 5 images on the left panel are the TOA radiance at 20:00 GMT, day 191, 192, 193, 194, 195, year 2004. The five images on the right hand panel are the topographically corrected PAR at corresponding times. The impact of topographic impact is more pronounced at absence of cloud than at presence of cloud. Mountainous areas experience more topographic impact than relatively flat areas. The white pixel on the TOA radiance images are clouds. The unit of TOA radiance is W/m2/rad/nm. The unit of instantaneous PAR is umol/m 2 /s. 94 The uncertainty involved in validating PAR value from 1km remote sensing data against ground-based point measurements is another issue that needs to be addressed. Concerns have been raised about using ground-based measurements to validate model performance at grid level (Li et al., 2005). More sophisticated validation methods, including inter- comparison among different remote sensing PAR products and comparison against meteorological reanalysis datasets, should be considered for future works. 4.6 Summary A new algorithm for estimating PAR using GOES visible band data is developed in this paper. Different from the existing GOES PAR algorithms, this new algorithm derives both surface reflectance and atmospheric parameters simultaneously from at-sensor radiance values, without the requirement for external knowledge of surface or atmospheric parameters. Validation against ground measurements indicates that this new algorithm is capable of reaching reasonably high accuracy. However, further improvement is necessary. The topographic correction has not been validated due to the lack of ground measurements. The error associated with the DEM used for topographic correction may introduce extra uncertainty into the topographically corrected PAR. Future studies should consider using topography data with greater accuracy in order to improve the accuracy of topographic correction. For instance, DEM generated by the shuttle radar topography mission (STRM) at 30m spatial resolution, 95 with global coverage between 60 north to 56 south latitude, is an alternative to the GTOPO30 DEM. The uncertainty involved in validating PAR value from 1km remote sensing data against ground-based point measurements is another issue that needs to be addressed. Concerns have been raised about using ground-based measurement to validate model performance at grid level (Li et al., 2005). More sophisticated validation methods, including inter- comparison among different remote sensing PAR products and comparison against meteorological reanalysis datasets, will be evaluated in future work. Research described in this chapter has been submitted to the Journal of Applied Meteorology and Climatology (Zheng et al., 2006). 96 Chapter 5 Conclusions PAR is an important variable in terrestrial ecosystem modeling because of its impact on green vegetations? photosynthesis. A complete PAR data product that spans multiples years at regional to continental scale is essential for model initialization, calibration, and validation. To supplement the current remote-sensing PAR products, this study has developed new algorithms to estimate PAR using remotely sensed data acquired by MODIS, AVHRR, and GOES. The new PAR algorithms are different from the existing ones in that they derive both atmospheric and land surface information solely from satellite TOA radiance. A common attribute of the PAR algorithms for MODIS, AVHRR, and GOES is that they are all built upon the near-monotonic change of PAR with changing TOA radiance. In addition, the algorithm for each of the three sensors has its own special features to accommodate each sensor?s unique characteristics. 5.1 Estimating PAR using MODIS data The MODIS sensors aboard NASA?s Terra and Aqua satellites possess many advanced features for the on-board calibration, geo-location, and channel position. Despite the technological superiority made available by the MODIS, no standard PAR product over land has been produced using MODIS data. To fill this gap, this study develops an algorithm to estimate incident instantaneous PAR using MODIS data. What set this new algorithm apart from existing algorithms is that it derives both the atmospheric and land surface information solely from satellite TOA radiance. 97 The PAR estimation using MODIS data consists of two main steps: first, retrieval of broadband surface reflectance using time-series satellite observations; second, retrieval of the atmospheric optical property parameters from MODIS TOA radiance and the subsequent retrieval of PAR. The theoretical basis of the retrieval of PAR from TOA radiance has been thoroughly examined and errors introduced by the simplification of atmosphere parameterization and linear interpolations of the look-up table entries have been analyzed and quantified. The new PAR algorithm has been applied to MODIS data from 2003 to 2004. An overall validation of the retrieved instantaneous PAR against ground measurements at six FLUXNET sites is given in Figure 5-1. Figure 5-1 Validation of the instantaneous PAR derived using MODIS data at the all six FLUXNET sites used in chapter 2. The solid line the is 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression. 98 5.2 Estimating PAR using AVHRR data As the continuation of developing methods to produce large scale PAR product on a robust basis, this study develops a new algorithm to estimate PAR using AVHRR data. Because of its long data record dating back to the late 1970s, AVHRR offers a unique attraction for remote sensing based PAR production: the opportunity to generate about 30-year PAR product. The new AVHRR PAR algorithm shares a similar theoretical basis with its MODIS counterpart: the derivation of surface reflectance and atmospheric parameters from TOA radiance. Results from MODTRAN4 simulations are used to examine the validity of the monotonic relation between PAR and AVHRR band-1 TOA radiance, which is the foundation of the algorithm. Application of the new algorithm to data acquired by AVHRR aboard the NOAA-14 satellite during 1996 and the validations (Figure 5-2) indicate the ability of the algorithm to reach reasonable accuracy in estimating instantaneous PAR. Furthermore, to meet the need for daily integrated PAR, the estimated instantaneous PAR from AVHRR data is scaled to daily PAR using an adjusted sinusoidal model. Although there is no physical mechanism in the adjusted sinusoidal model to account for the impact of abrupt changes in the atmospheric condition on PAR, the performance of the model is reasonably good: the validations against ground measurements indicate that the estimated daily integrated 99 PAR reaches higher accuracy than their instantaneous counterpart at all of the three validation sites. Figure 5-2 Validation of the instantaneous PAR derived using AVHRR data at the three FLUXNET sites used in chapter 3. The solid line is the 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression. 5.3 Estimating PAR using GOES data In scaling from instantaneous to daily integrated PAR, a larger number of instantaneous PAR tends to result in greater scaling accuracy because it better captures the PAR diurnal cycle better. As a result of this consideration, this study extends the MODIS/AVHRR PAR algorithm to the data acquired by the imager aboard on the geostationary satellite GOES. Compared with the MOIDS and AVHRR, the GOES PAR algorithm features two important improvements: excluding cloud-shaded pixels utilizing the spatial relationship and accounting for the surface bidirectional reflectance characteristics by fitting a three- 100 parameter linear Ross-Li model. It is the high temporal resolution of GOES data, which provides satellite observations with varying illuminating/viewing geometry over a daytime, that make it possible to incorporate the BRDF model in the GOES algorithm. Figure 5-3 Validation of the instantaneous PAR derived using GOES data at the four FLUXNET sites used in chapter 4. The solid line is the 1:1 line, and the dashed line is the fitted line using the ordinary least squares regression. Similar to that of MODIS and AVHRR, the theoretical basis of the GOES algorithm has been thoroughly examined and error budget analyses have been conducted. Application of the algorithm to data acquired by GOES-12 in 1999 resulted in estimated PAR covering the continental United States. It is evident from visual inspection that the new algorithm successfully captured the dynamics of PAR cause by the cloud impacts and solar angle changes. The overall validation results, including the ground measurements obtained from four FLUX network, is depicted in Figure 5-3. 101 Furthermore, to meet the need of terrestrial ecosystem modeling, PAR estimated using GOES have been corrected for the topographic impacts. The direct and diffuse components of PAR were corrected for the angular effect, sky viewing factor, and radiation from neighboring terrains. Table 5-1 Statistical summary of validation results of instantaneous PAR estimated using MODIS, AVHRR, and GOES. N is number of samples. MODIS GOES AVHRR R 2 0.9011 0.9131 0.8788 Slope 0.8771 1.0327 0.9272 Intercept 103.9378 -28.0842 25.2575 Bias 32.4876 -6.0201 47.50 RMSE 194.815 161.130 195.478 Relative Error (%) 11.763 16.207 31.85 RMSE/Average (%) 17.494 15.33 20.527 N 1496 1775 1529 5.4 Combining MODIS, AVHRR, and GOES data for daily PAR estimation Because the primary use of PAR products is for the terrestrial ecosystem modeling, and terrestrial models typically require long-term PAR data, it is desirable to generate PAR products with extensive coverage both spatially and temporally. To meet this goal, this study utilizes the complementary characteristics of the remotely sensed data acquired by three different sensors: MODIS, AVHRR, and GOES. 102 The validation results of the instantaneous PAR derived from the three sensors (Table 5- 1) indicate that MODIS derived PAR reaches the higher accuracy (relative error of 11.76%) than AVHRR (31.85%) and GOES (16.21%). In terms of estimation accuracy of instantaneous PAR, all of the three algorithms need to be improved. In addition, what is more relevant to the terrestrial biological modeling is the accuracy of the daily PAR estimation. The integration of the daily integrated PAR from the instantaneous PAR derived using AVHRR data shows that daily PAR estimation reaches a higher accuracy (relative error at 24.22%) than the instantaneous PAR (relative error at 31.85%). One approach for achieving higher accuracy in estimating daily PAR is to combine the instantaneous PAR estimations from all of the three sensors: the different overpass times of the three sensors provide more snapshots of the PAR diurnal variation and therefore presumably provide a better physical basis for the daily PAR integration. However due their differences in spatial, temporal, and radiometric properties, it is challenging to develop daily PAR integration methods that effectively accommodate and utilize each sensor?s characteristics. Furthermore, the different assumptions made in the PAR algorithms for the three sensors also make it difficult to integrate the three instantaneous PAR results. Although more efforts are needed, integrating daily PAR using instantaneous PAR derived from multiple sensors constitutes a promising research direction. In addition, using multiple sensors is necessary to generate PAR product at large spatial scale. For example, to encompass the entire North America continent, GOES alone is not sufficient 103 because its coverage does not go beyond 66 degree latitude. Furthermore, the spatial resolution becomes very coarse as GOES observations approach the high-latitude. One possible alternative, which has been adopted in this study, is to take advantage of polar orbiting satellites. Polar orbiting satellites, such as AVHRR and MODIS, have more frequent overpasses at the higher latitude than at lower latitude, and therefore are able to generate more frequent instantaneous PAR estimates, which enable a better daily PAR integration. 5.5 Issues for future research Through the study reported in this dissertation, several issues have been identified for the future exploration. The first is to improve the instantaneous-to-daily PAR scaling approach. Although the adjusted sinusoidal model achieves a reasonable accuracy in this study, there is need and room for improvements. One possible direction for such improvements could be through using instantaneous PAR estimates from different sensors in integrating daily PAR. Combined, AVHARR and MODIS have more daytime overpasses than when being used alone, thus providing a better chance to capture the PAR diurnal variation. The second issue is the utilization of the sea-viewing wide field of view satellite (SeaWiFS) in estimating PAR. The addition of SeaWiFS will further increase the number of instantaneous PAR estimates per day, providing the possibility to improve the daily PAR integration. 104 The third issue is related to topographic correction of PAR. This study carried out topographic correction on the estimated PAR, but did not validate the results due to the lack of topographically corrected PAR ground measurements. Given the importance of correcting PAR for the topographic impact, greater efforts should be made in the future to collect PAR ground measurements with topographic correction. The new algorithms developed in this study enable the generation of PAR product on the continental scale over extended time period, and therefore provide the platform for examining the spatial and temporal variation of PAR. Future research in this direction may link PAR variation with the natural or anthropogenic activities, which may shed light on the problem of global warming and climate change. Finally, the finer spatial resolution PAR products generated from this study should be aggregated and compared with coarser resolution PAR products to study the scale on which PAR varies. 105 Reference Al-Shooshan, A. A., 1997: Estimation of photosynthetically active radiation under an arid climate. Journal of agricultural engineering research, 66, 9-13. Alados-Arboleda, L., F. Olmo, J, I. Alados, and M. Perez, 2000: Parametric models to estimate photosynthetically active radiation in spain. Agricultural and Forest Meteorology, 101, 187-201. Alados, I., F. J. Foyo-moreno, and L. Alados-Arboleda, 1996: Photosynthetically active radiation: measurement and modelling. Agricultural and Forest Meteorology, 78, 121-131. Alados, I., I. Foyo-Moreno, F. J. Olmo, and L. Alados-Arboledas, 2002: Improved estimation of diffuse photosynthetically active radiation using two spectral models. Agricultural and Forest Meteorology, 111, 1-12. Alados, I., I. Foyo-Moreno, F. J. Olmo, and L. Alados-Arboledas, 2003: Relationship between net radiation and solar radiation for semi-arid shrub-land. Agricultural and Forest Meteorology, 116, 221-227. Asrar, G., 1989: Theory and application of optical remote sensing. John Wiley and Sons, Hoboken, NJ, New York Beringer, F., 1994., 1994: Simulated irradiance and temperature estimates as a possible source of bias in the simulation of photosynthesis. Agricultural and Forest Meteorology, 71, 19-32. Berk, A., L. S. Bernstein, G. P. Anderson, P. K. Acharya, D. C. Robertson, J. H. Chetwynd, and S. M. Adler-Golden, 1998: MODTRAN cloud and multiple 106 scattering upgrades with application to AVIRIS. Remote Sensing of Environment, 65, 367-375. Beyer, H. G., C. Costanzo, and D. Heinemann, 1996: Modifications of the Heliosat procedure for irradiance estimates from satellite images. Solar Energy, 56, 207- 212. Beyer, H. G., A. Hammer, and D. Heinemann, 1997: Satellite-derived irradiance statistics for Africa. Solar Energy, 59, 233-240. Cano, D., J. M. Monger, M. Albuisson, H. Guillard, N. Regas, and L. Wald, 1986: A method for the determination of the global solar radiation from meteorological satellite data. Solar Energy, 37, 31-39. Carder, K., S. Hawes, and R. Chen, 1999: Instantaneous photosynthetically available radiation and absorbed radiation by phytoplankton, version 5. ATBD-MOD-20. Ceballos, J. C., M. J. Bottino, and J. M. de Souza, 2004: A simplified physical model for assessing solar radiation over Brazil using GOES 8 visible imagery. Journal of Geophysical Research-Atmospheres, 109, D02211, doi:10.1029/2003JD003531. Chandrasekhar, S., 1960: Radiative Transfer. Dover Publications, Mineola, NY, Charlock, T. P. and T. L. Alberta, 1996: The CERES/ARM/GEWEX experiment (CAGEX) for the retrieval of radiative fluxes with satellite data. Bulletin of the American Meteorological Society, 77, 2673-2683. Corripio, J. G., 2003: Vectorial algebra algorithms for calculating terrain parameters from DEMs and solar radiation modelling in mountainous terrain. International Journal of Geographical Information Science, 17, 1-23. 107 Dai, Y., R. E. Dickinson, and Y.-P. Wang, 2004: A two-big-leaf model for canopy temperature, photosnythesis and stomatal conductance. Journal of Climate, 17, 2281-2299. Diak, G. R. and C. Gautier, 1983: Improvements to a simple physical model for estimating insolation from Goes data. Journal of Climate and Applied Meteorology, 22, 505-508. Diner, D. J., J. C. Beckert, and e. al, 1998: Multi-angle imaging spectroradiometer (MISR) instrument description and experiment overview. IEEE Geoscience and Remote Sensing, 36, 1072-1087. Dubayah, R., 1992: Estimating net solar radiation using Landsat Thematic Mapper and digital elevation data. Water Resources Research, 28, 2469-2484. Dubayah, R. and P. M. Rich, 1995: Topographic solar-radiation models for GIS. International Journal of Geographical Information Systems, 9, 405-419. Dubayah, R. and S. Loechel, 1997: Modeling topographic solar radiation using GOES data. Journal of Applied Meteorology, 36, 141-154. Duguay, C. R., 1995: An approach to the estimation of surface net-radiation in mountain areas using remote-sensing and digital terrain data. Theoretical and Applied Climatology, 52, 55-68. Dye, D. G., 2004: Spectral composition and quanta-to-energy ratio of diffuse photosynthetically active radiation under diverse cloud conditions. Journal of Geophysical Research, 109, doi:10.1029/2003JD004251. 108 Eck, T. F. and D. G. Dye, 1991: Satellite estimation of incident photosynthetically active radiation using ultraviolet reflectance. Remote Sensing of Environment, 38, 135- 146. Frouin, R. and R. T. Pinker, 1995: Estimating Photosynthetically Active Radiation (PAR) at the Earth's surface from satellite observations. Remote Sensing of Environment, 51, 98-107. Frouin, R., B. Franz, and M. Wang, 2000: Algorithm to estimate PAR from SeaWiFS data, version 1.0 Documentation. Gautier, C., G. Diak, and S. Masse, 1980: A simple physical model to estimate incident solar radiation at the surface from GOES satellite data. Journal of Applied Meteorology, 19, 1005-1012. Goward, S. N. and K. F. Huemmrich, 1992: Vegetation canopy PAR absorptance and the normalized difference vegatation index: an assessment using the SAIL model. Remote Sensing of Environment, 39, 119-140. Goward, S. N., C. J. Tucker, and D. G. Dye, 1985: North American vegetation patterns observed with the NOAA-7 Advanced Very High Resolution Radiomter. Vegetation, 64, 3-14. Gu, J. J. and E. A. Smith, 1997: High-resolution estimates of total solar and PAR surface fluxes over large-scale BOREAS study area from GOES measurements. Journal of Geophysical Research, 102, 29685-29705. Gu, J. J., E. A. Smith, and J. D. Merritt, 1999: Testing energy balance closure with GOES-retrieved net radiation and in situ measured eddy correlation fluxes in BOREAS. Journal of Geophysical Research, 104, 27881-27893. 109 Gu, J. J., E. A. Smith, and H. J. Cooper, 2004: Modeling carbon sequestration over the large-scale Amazon basin, aided by satellite observation. Part I: wet- and dry- season surface radiation budget flux and precipitation variability based on GOES retrievals. Journal of Applied Meteorology, 43, 870-886. Gu, L. H., D. Baldocchi, S. B. Verma, T. A. Black, T. Vesala, E. M. Falge, and P. R. Dowty, 2002: Advantages of diffuse radiation for terrestrial ecosystem productivity. Journal of Geophysical Research-Atmospheres, 107, 4050, doi:10.1029/2001JD001242. Gu, L. H., D. D. Baldocchi, S. C. Wofsy, J. W. Munger, J. J. Michalsky, S. P. Urbanski, and T. A. Boden, 2003: Response of a deciduous forest to the Mount Pinatubo eruption: enhanced photosynthesis. Science, 299, 2035-2038. Gueymard, C., 1989: An atmospheric transmittance model for the calculation of the clear sky beam, diffuse and global photosynthetically active radiation. Agricultural and Forest Meteorology, 45, 215-229. Gueymard, C., 1993: Critical analysis and performance assessment of clear sky irradiance models using theoretical and measured data. Solar Energy, 51, 121-138. Hall, D. K., G. A. Riggs, and V. V. Salomonson, 2001: Algorithm theoretical basis document (ATBD) for the MODIS snow and sea ice-mapping algorithms. Hicks, B. B., J. J. DeLuisi, and D. R. Matt, 1996: The NOAA Integrated Surface Irradiance Study (ISIS) - A new surface radiation monitoring program. Bulletin of the American Meteorological Society, 77, 2857-2864. Iqbal, M., 1983: An introduciton to solar radiation. Academic Press, Orlando, FL, 169- 213. 110 Jacobs, J. M., D. A. Myers, M. C. Anderson, and G. R. Diak, 2002: GOES surface insolation to estimate wetlands evapotranspiration. Journal of Hydrology, 266, 53-65. Kaufman, Y. J., 1993: Aerosol optical-thickness and atmospheric path radiance. Journal of Geophysical Research, 98, 2677-2692. Kobayashi, H. and T. Matsnaga, 2004: Sattellite estimation of photosynthetically active radiation in southeast Asia: Impacts of smoke and cloud cover. Journal of Geophysical Research, 109, doi:10.1029/2003JD003807. Kumar, L., A. K. Skidmore, and E. Knowles, 1997: Modelling topographic variation in solar radiation in a GIS environment. International Journal of Geographical Information Science, 11, 475-497. Landsberg, J. J., S. D. Prince, P. G. Jarvis, R. E. McMurtrie, R. Luxmoore, and B. E. Medlyn, 1996: Energey coversion and use in forests: an analysis of forest production in terms of radiation utilisation efficiency. The use of remote sensing in the modeling of forest productivity, Kluwer Academic Press, Dordrecht, The Netherlands Lavigne, M. B. and M. G. Ryan, 1997: Growth and maintenance respiration rates of aspen, black spruce, and jack pine stems at northern and sourthern BOREAS sites. Tree Physiology, 17, 543-551. Li, Z. Q., M. C. Cribb, and F. L. Chang, 2005: Natural variability and sampling errors in solar radiation measurements for model validation over the atmospheric radiation measurement southern great plains region. Journal of Geophysical Research, 110, D15S19, doi:10.1029/2004JD005028. 111 Liang, S., 2004: Quantitative remote sensing of land surface. John Wiley and Sons, Hoboken, NJ. Liang, S., T. Zheng, R. G. Liu, H. L. Fang, S. C. Tsay, and S. Running, 2006: Estimation of incident photosynthetically active radiation from Moderate Resolution Imaging Spectrometer data. Journal of Geophysical Research, 111, D15208, doi:10.1029/2005JD006730 Lopez, G., M. A. Rubio, M. Martinez, and F. J. Batlles, 2001: Estimation of hourly global photosynthetically active radiation using artifical neural network models. Agricultural and Forest Meteorology, 107, 279-291. Maier, C. A., S. J. Zarnoch, and P. M. Dougherty, 1998: Effects of temperature and tissue nitrogen on dormant season stem and branch maintenance respiration in a young loblolly pine plantation. Tree Physiology, 18, 11-20. Mottus, M. and J. R. Sulev, 2001: Experimental study of ratio of PAR to direct integral solar radiation under cloudless conditions. Agricultural and Forest Meteorology, 109, 161-170. Myneni, R., S. O. Los, and G. Asrar, 1995: Potential gross primary productivity of terrestrial vegetation from 1982-1990. Geophysical Research Letters, 22, 2617- 2620. Myneni, R. B., R. R. Nemani, and S. W. Running, 1997: Estimation of global leaf area index and absorbed par using radiative transfer models. IEEE Transactions on Geoscience and Remote Sensing, 35, 1380-1393. 112 Ohmura, A., 1998: Baseline surface radiation network (BSRN/WCRP): new precision radiometry for climate research. Bulletin of American Meterological Society, 2115-2136. Oktin, J. A., M. C. Anderson, J. R. Mecikalski, and G. R. Diak, 2005: Validation of GOES-based insolation estimates using data from the US Climate Reference Network. Jounal of Hydrometeorology, 6, 460-475. Pinker, R. T. and I. Laszlo, 1992: Global distribution of photosynthetically active radiation as observed from satellites. Journal of Climate, 5, 56-65. Pinker, R. T., I. Laszlo, C. H. Whitlock, and T. P. Charlock, 1992: Modeling surface solar irradiance for satellite applications on global scale. Journal of Applied Meteorology, 31, 194-211. Pinker, R. T., R. Frouin, and Z. LI, 1995: A review of satellite methods to derive surface shortwave irradiance. Remote Sensing of Environment, 51, 108-124. Prince, S. D., 1991: A model of regional primary production for use with coarse resolution satellite data. International Journal of Remote Sensing, 12, 1313-1330. Prince, S. D. and S. N. Goward, 1995: Global primary production: a remote sensing approach. Journal of Biogeography, 22, 815-835. Ramsay, B. H., 1998: The interactive multisensor snow and ice mapping system. Hydrological Processes, 12, 1537-1546. Reich, P. B., B. D. Kloeppel, D. S. Ellsworth, and M. B. Walters, 1995: Different photsynthesis-nitrogen relations in deciduous hardwood and evergreen coniferous tree species. Oecologia, 104, 24-30. 113 Ross, J. and J. R. Sulev, 2000: Sources of errors in measurements of PAR. Agricultural and Forest Meteorology, 100, 103-125. Rubio, M. A., G. Lopez, J. Tovar, D. Pozo, and F. J. Battles, 2004: The use of satellite measurements to estimate photosynthetically active radiation. Physics and Chemistry of the Earth, 30, 159-164. Running, S. W., R. R. Nemani, F. A. Heinsch, M. S. Zhao, M. Reeves, and H. Hashimoto, 2004: A continuous satellite-derived measure of global terrestrial primary production. Bioscience, 54, 547-560. Schwarz, P. A., T. J. Fahey, and T. E. Dawson, 1997: Seasonal air and soil temperature effects on photosynthesis in red spruce saplings. Tree Physiology, 17, 187-194. Sellers, P. J., 1985: Canopy reflectance, photosynthesis, and transpiration. International Journal of Remote Sensing, 6, 1335-1371. Sellers, P. J. and D. S. Schimel, 1993: Remote sensing of the land biosphere and biogeochemistry in the EOS era: Science priorities, methods and implementation- EOS land bioshpere and biogeochemical cycles panels. Global and Planetary Change, 7, 279-297. Simpson, J. J. and J. R. Stitt, 1998: A procedure for the detection and removal of cloud shadow from AVHRR data over land. IEEE Transactions on Geoscience and Remote Sensing, 36, 880-897. Simpson, J. J., Z. H. Jin, and J. R. Stitt, 2000: Cloud shadow detection under arbitrary viewing and illumination conditions. IEEE Transactions on Geoscience and Remote Sensing, 38, 972-976. 114 Stewart, J. B., C. J. Watts, J. C. Rodriguez, H. A. R. De Bruin, A. R. van den Berg, and J. Garatuza-Payan, 1999: Use of satellite data to estimate radiation and evaporation for northwest Mexico. Agricultural Water Management, 38, 181-193. Van Laake, P. E. and G. A. Sanchez-Azofeifa, 2004: Simplified atmospheric radiative transfer modelling for estimating incident PAR using MODIS atmosphere products. Remote Sensing of Environment, 91, 98-113. Van Laake, P. E. and G. A. Sanchez-Azofeifa, 2005: Mapping PAR using MODIS atmosphere products. Remote Sensing of Environment, 94, 554-563. Varley, M. J., K. J. Beven, and H. R. Oliver, 1996: Modelling solar radiation in steeply sloping terrain. International Journal of Climatology, 16, 93-104. Vermote, E. F., J. Tanre, and e. al, 1997: Second simulation of the satellite signal in the solar spectrum, 6s: An overview. IEEE Geoscience and Remote Sensing, 35, 675- 686. Wang, D. D., S. L. Liang, and T. Zheng, 2007: Estimation of Daily-intergrated PAR from MODIS data. Letter of IEEE Transactions on Geoscience and Remote Sensing, Submitted. Wang, K., X. Zhou, J. Liu, and M. Sparrows, 2005: Estimating surface solar radiation over complex terrain using moderate-resolution satellite sensor data. International Journal of Remote Sensing, 26, 47-58. Will, R. E. and R. O. Teskey, 1997: Effect of elevated carbon dioxide concentration and root restriction in net photosynthesis, water relations and foliar carbohydrate status of loblolly pine seedling. Tree Physiology, 17, 655-661. 115 Winslow, J. C., E. R. Hunt, and S. C. Piper, 2001: A globally applicable model of daily solar irradiance estimated from air temperature and precipitation data. Ecological Modelling, 143, 227-243. Zheng, T., S. Liang, and K. C. Wang, 2006: Estimation of incident Photosynthetically Active Radiation from GOES visible imagery. Journal of Applied Meteorology and Climatology. In revision.