A New Set of MODIS Land Products (MCD18): Downward Shortwave Radiation and Photosynthetically Active Radiation

: Surface downward shortwave radiation (DSR) and photosynthetically active radiation (PAR), its visible component, are key parameters needed for many land process models and terrestrial applications. Most existing DSR and PAR products were developed for climate studies and therefore have coarse spatial resolutions, which cannot satisfy the requirements of many applications. This paper introduces a new global high-resolution product of DSR (MCD18A1) and PAR (MCD18A2) over land surfaces using the MODIS data. The current version is Collection 6.0 at the spatial resolution of 5 km and two temporal resolutions (instantaneous and three-hour). A look-up table (LUT) based retrieval approach was chosen as the main operational algorithm so as to generate the products from the MODIS top-of-atmosphere (TOA) reﬂectance and other ancillary data sets. The new MCD18 products are archived and distributed via NASA’s Land Processes Distributed Active Archive Center (LP DAAC). The products have been validated based on one year of ground radiation measurements at 33 Baseline Surface Radiation Network (BSRN) and 25 AmeriFlux stations. The instantaneous DSR has a bias of − 15.4 W / m 2 and root mean square error (RMSE) of 101.0 W / m 2 , while the instantaneous PAR has a bias of − 0.6 W / m 2 and RMSE of 45.7 W / m 2 . RMSE of daily DSR is 32.3 W / m 2 , and that of the daily PAR is 13.1 W / m 2 . The accuracy of the new MODIS daily DSR data is higher than the GLASS product and lower than the CERES product, while the latter incorporates additional geostationary data with better capturing DSR diurnal variability. MCD18 products are currently under reprocessing and the new version (Collection 6.1) will provide improved spatial resolution (1 km) and accuracy.


Background
Solar radiation, mainly in the shortwave spectrum, is the primary energy resource of the Earth system and, consequently, affects many processes of the system [1]. Downward shortwave radiation (DSR) that reaches the land surface is of particular importance to the terrestrial ecosystem and human society not only because of its role in surface radiation and energy budget but also because it is the driving force of many land processes [2]. It also affects the weather and climate processes through regulating exchanges of energy, carbon, water, other matter, and momentum between the surface and the atmosphere [3]. The visible component of DSR (with wavelengths between 400 and 700 nm), known as photosynthetically active radiation (PAR), is the amount of energy available for photosynthesis, and plays a significant role in ecosystem productivity and carbon assimilation [4]. DSR also directly determines the amount of energy that can be used as clean renewable energy.
Extraterrestrial shortwave radiation at the top of atmosphere (TOA) varies mainly with latitude and day of year since the two parameters determine the Earth-to-Sun distance, solar zenith angle, and daytime length. Factors affecting DSR and PAR at the Earth's surface are more complicated. Other than those two parameters, atmospheric conditions strongly regulate DSR and PAR at the surface because of scattering and absorption of various atmospheric constituents [5]. Moreover, surface reflectance is also a factor affecting DSR and PAR due to the multiscattering between the atmosphere and the surface [6].
DSR and PAR can be estimated from satellite data using several different strategies. One approach is to calculate DSR and PAR from data of atmospheric and surface parameters using atmospheric radiative transfer models or simplified parameterization schemes derived from radiative transfer simulations [7]. The other strategy is to directly retrieve DSR and PAR from satellite TOA signatures. The retrieval algorithms can be either based on an empirical relationship between TOA reflectance and surface radiative fluxes [8] or predefined look-up table that is established from radiative transfer simulation [9] or other more complicated inversion algorithms such as optimization [10].
Several global data sets of DSR and/or PAR have been generated using some of the abovementioned algorithms from various sources of satellite data or products. The products were mainly developed for the weather and climate community and bear rather coarse spatial resolutions, which cannot satisfy the requirement of the terrestrial research and applications on spatial details [11]. Additionally, many of the existing products have stopped updates and covered only historical periods. An operational global high resolution product of land surface shortwave radiative fluxes is still limited.
NASA's MODIS land science team has generated 16 global products of key land surface parameters, such as surface albedo, leaf area index, land surface temperature, and so on [12]. These products have been extensively used to address a wide range of application and scientific issues from regional to global scales. Since 2017, NASA has started to produce the new set of MODIS DSR and PAR products (MCD18) over land surfaces. This paper summarizes the retrieval algorithm in theoretical and practical perspectives and presents an initial assessment of the new product suite. Following this background, Section 2 briefly introduces the algorithm. Section 3 summarizes the data processing system and describes the data products. Section 4 presents the validation results and Section 5 concludes the paper with discussions and outlook.

Algorithm Description
The look-up table (LUT) based shortwave irradiance retrieval algorithm is used to generate the MODIS DSR and PAR products. The LUT approach takes the satellite TOA signature as the main information source to estimate atmospheric radiative properties and does not rely on additional products of aerosol or cloud optical parameters, which typically contain a substantial amount of missing values. By directly using TOA data, the LUT approach has the advantage of performing at the native resolution of the satellite data and producing fewer data gaps. Through more than one decade of development, the LUT approach has become mature and robust.
The LUT approach was initially developed by Liang et al. [9] to estimate instantaneous PAR from MODIS data. A temporal scaling scheme was added in a following study so that it can estimate temporally averaged or integrated quantities [13]. The absorption effects of water vapor in shortwave broadband radiation was later incorporated to extend the LUT approach capability in estimating DSR. In its initial version, surface reflectance was estimated from time series of MODIS TOA reflectance data with minimal blue reflectance approach. This approach was later replaced by the direct use of surface reflectance product for simplicity and also to avoid the issue of cloud shadow contamination [14]. Other than MODIS data, the LUT approach has also been adapted to other polar-orbiting data such as AVHRR and geostationary satellite data such as GOES and MTSAT [15]. Combining MODIS and geostationary data, Zhang et al. [16] generated a global data set of DSR and PAR for three years (2008-2010) using the LUT-based approach. In their study, the data products were generated on a yearly basis because a time-series based filter was needed to smooth the annual surface reflectance data. The near real time MODIS surface albedo product is used in this paper to reduce the latency of data production. In Zhang's study, a post correction of water vapor absorption was applied. Here, the effects of water vapor on both TOA spectral reflectance and surface fluxes are fully considered through adding water vapor as an additional dimension of LUT. The theoretical basis of the refined the LUT approach is briefly summarized as follows.
Two parameterization schemes are used to simplify the interactions of electromagnetic radiation between the atmosphere and the surface with the improved efficiency of modeling the atmospheric radiative transfer process [17]. The first formula is to parameterize TOA spectral reflectance. TOA spectral reflectance R(λ) for a given wavelength λ can be obtained from surface spectral reflectance r(λ) and three parameters related to atmospheric radiative transfer: where θ s is solar zenith angle (SZA). R 0 (λ), ρ(λ), and γ(λ) are path reflectance, atmospheric spherical albedo, and transmittance for the spectral band, respectively. The second parameterization scheme is to formulate surface incident broadband radiation fluxes (i.e., DSR for the total shortwave spectrum and PAR for its visible component). Broadband downward irradiance in the LUT approach is calculated using the following formula: where r is surface broadband reflectance, E 0 is extraterrestrial solar broadband irradiance, θ s is SZA, F 0 is path irradiance, ρ is atmospheric spherical albedo, and γ is atmospheric transmittance. The parameters related to the atmospheric radiative transfer in both equations are simulated offline with the numerical atmospheric radiative transfer models. The parameters in Equation (1) used for TOA reflectance calculation are saved in the reflectance LUT and the counterpart for flux retrieval in Equation (2) are saved in the flux LUT.
The radiative transfer parameters are functions of observation geometry and atmospheric conditions. Observation geometry in Equation (1) refers to SZA, view zenith angle (VZA), and relative azimuth angle (RAA), while that in Equation (2) includes only SZA. The factors that affect radiation scattering and absorption in the atmosphere are complex, including the quantity, type, and vertical distribution of all kinds of atmospheric constituents. In the LUT approach, Rayleigh scattering and gas absorption are determined mainly by the air pressure (air mass), which is simply represented by surface elevation. The data of total precipitable water is used in the LUT approach to quantify absorption of water vapor. Other than the abovementioned variables, aerosol and cloud are the major factors that affect TOA radiance received by satellite sensors and incident irradiance reaching the land surface. Atmospheric optical depth, which refers to aerosol optical depth for the clear-sky cases and cloud optical depth for the cloudy-sky cases, is used as a general variable to calculate atmospheric scattering and absorption.
The LUT approach consists of two major steps: estimation of atmospheric optical depth and calculation of DSR and PAR. A simple optimization algorithm is used in the first step to estimate atmospheric optical depth. With surface spectral reflectance and water vapor known, TOA spectral reflectance of a given value of atmospheric optical depth can be predicted from Equation (1). The predicted reflectance will be compared with the observed TOA reflectance. The value of atmospheric optical depth that produces the best match between the predicted and the observed reflectance is chosen. In the second step, the LUT approach uses Equation (2) to calculate surface irradiance. Equation (2) contains five parameters. Among them, view geometry and surface reflectance are obtained from the input data. The atmospheric optical depth from the first step, together with ancillary input data of water vapor and surface elevation, is used to select appropriate entries in the flux LUT and get the remaining three radiative transfer parameters in Equation (2).

Data Processing and Product Generation
The LUT approach uses MODIS L1b TOA reflectance, geolocation, MCD43A3 surface albedo, MERRA2 water vapor, GTOPO30 DEM, and surface reflectance climatology data as the input data (Table 1). MCD43A3 500 m data is aggregated to 1 km to match with the spatial resolution of TOA input data. The current version of the LUT approach assumes Lambertian surface and surface albedo data is used as alternative of surface reflectance. Surface reflectance climatology is generated from the historical MODIS surface albedo product and calculated as the mean of 10 years' albedo. The climatology value will be used as the backup value of surface reflectance when the MCD43A3 data is not available. The LUT approach is executed at the daily time step. The L1b reflectance and geolocation data of one given day are first gridded into the MODIS sinusoidal map projection and all the swaths covering one given tile are stacked together and organized as one tiled input TOA file (MOD18IPGD for Terra and MYD18IPGD for Aqua). Both Terra and Aqua MODIS data of one given day are used to increase the number of observations. Each MODIS TOA instantaneous observation is fed into the LUT approach as the input of the TOA information to estimate instantaneous atmospheric optical depth. The input of the surface reflectance is from the MCD43A3 product or the surface reflectance climatology when the valid MCD43A3 data is not available. The reflectance and flux LUTs are multidimensional arrays ( Table 2). Data of view geometry, elevation, and total precipitable water are known and used to reduce the dimensions of LUTs. The LUTs are linearly interpolated according to the actual values of these variables. After the dimension reduction, the single dimension of the LUTs is atmospheric optical depth, which is the only unknown to be estimated in the first step. The results of atmospheric optical depth derived from the first step are the instantaneous values corresponding to the associated MODIS overpass time. The instantaneous optical depth data are then fed into the second step of the LUT approach to estimate instantaneous DSR and PAR values. For the convenience of the users, a nearest neighborhood interpolation approach is used to produce a 3-h time series of optical depth, which is fed again into the second step of the LUT approach to generate 3-h DSR and PAR values. The MCD18 product suite includes two data products: DSR (MCD18A1) and PAR (MCD18A2). The products are archived in Hierarchical Data Format V4-Earth Observing System (HDF-EOS) format files. MCD18A1 and MCD18A2 are daily gridded Level 3 products in the MODIS sinusoidal map projection. Each MCD18 file contains two major types of scientific data sets: instantaneous DSR or PAR arrays for each individual MODIS overpass and 3-h radiation arrays at 00:00, 03:00, 06:00, 09:00, 12:00, 15:00, 18:00, 21:00 in Coordinated Universal Time (UTC). Besides, an array of quality flag is used to indicate the input source of surface reflectance data. The MCD18 Collection 6.0 products are produced at the spatial resolution of 5 km, while the new Collection 6.1 products will be available at an improved 1 km spatial resolution. In the Collection 6.1, additional Climate Modeling Grid (CMG) products (MCD18C1 for DSR and MCD18C2 for PAR) are generated at the spatial grid of 0.05 • in the latitude/longitude map projection. An example of global maps of 3-hourly DSR product is shown in Figure 1. It is clearly seen that large amounts of solar radiation move from the east to the west with the Earth rotation.

Field Data and Processing
Field measurements of DSR in 2018 collected at 33 Baseline Surface Radiation Network (BSRN) stations and those of PAR from 25 AmeriFlux stations are used to quantitatively validate the MCD18 products. The Baseline Surface Radiation Network (BSRN) is an international collaboration network established by World Climate Research Programme (WCRP) to collect, archive, and distribute high quality measurements of surface irradiance [18]. DSR is one of BSRN basic measurements, which all stations are obligated to monitor. AmeriFlux is a regional partner of the global FLUXNET network. At the AmeriFlux sites, flux towers are set to measure the exchange rates of matter, radiation, and energy between the atmosphere and the biosphere, as well as meteorological and other ancillary data [19]. Different from that from BSRN, PAR is measured as the number of photosynthetic photons in the unit of umol/m 2 /s at the AmeriFlux sites. Division by a coefficient of 4.56 is used to convert the AmeriFlux PAR measurements from umol/m 2 /s to W/m 2 [20], which the MCD18 product has.
BSRN measures surface radiation fluxes at high sampling frequency and in situ data are available at an interval of 1 min. The station data within a 60-min difference of the satellite acquaintance time are averaged to compare with the instantaneous satellite estimates to account for the spatiotemporal mismatch between the satellite products and the ground measurements [21]. To handle the issue of missing values, the daily mean of ground measurements is calculated in two steps. The 30-min mean values are first calculated if there is at least one valid measurement during the 30-min time window. The daily mean values are then computed as the average of the 30-min means of a given day only when there is no missing 30-min mean data in that day. AmeriFlux data are available at coarser temporal resolution with a half hour interval and a similar procedure is used to obtain instantaneous and daily values that match with the satellite products. R 2 , root mean square error (RMSE), bias, and relative RMSE are calculated as the validation metrics.

Validation of Instantaneous Product
The comparison results of all the MCD18A1 instantaneous values against the field measurements over the BSRN stations are shown in Figure 2. The satellite estimated DSR shows strong correlation with the in situ data with an R 2 value of 0.860. A slight underestimation with a negative bias value of 18.2 W/m 2 is observed with this validation data set. The overall RMSE is 119.2 W/m 2 . The accuracy is comparable to another MODIS radiation product which uses atmosphere and land products as inputs [22]. Previous studies found spatial aggregation of high resolution satellite data can improve its agreement with the in situ data [23]. Being averaged to a 50-by-50 km grid, the MCD18A1 product does see RMSE reduces to 101.0 W/m 2 and R 2 increases to 0.889 (Figure 3). The results reported in this study refer to those aggregated to a 50-by-50 km grid unless otherwise stated.
To quantify the impacts of surface reflectance on retrieval results, the MCD18A1 data are separated by the availability of MODIS surface reflectance (Figure 4). The quality of DSR retrievals is higher when there is valid MODIS surface reflectance input, with a RMSE of 89.6 W/m 2 , bias of -7.7 W/m 2 . The DSR accuracy is greatly degraded when no actual surface reflectance product is available and climatological values are used. In such cases, R 2 decreases to 0.806 and RMSE increases to 131.7 W/m 2 . Besides, DSR retrieved with surface reflectance climatology is substantially underestimated with a negative bias of 43.7 W/m 2 . A further analysis discloses the underestimation mainly occurs over the snow-covered conditions, where a lower value of surface reflectance is obtained from climatology and the clear-sky pixels are mistakenly treated as clouds because of the high TOA reflectance which is actually the result of snow cover. However, the use of surface reflectance climatology does significantly improve the frequency of valid DSR estimation. For the validation data used in this study, the number of valid DSR values is improved by 18.3% from with only valid MCD43A3 data to with the optional climatology data.    Table 3 lists the summary statistics of the DSR validation results at the individual BSRN sites. Over most of the 33 sites, the MODIS DSR product demonstrates high accuracy when compared with the field measurements. At the GOB site, RMSE is only 53.9 W/m 2 , and the relative error is as low as 6.8%. A total of 22 out of 33 sites have relative RMSE smaller than 20%. There are five sites with relative RMSE greater than 30%, which are located mainly at high latitudes or high altitudes with frequent snow coverage. This indicates the single band information used in the current algorithm is limited in discriminating the bright targets between cloud and snow.

Validation of 3-h and Daily Products
The interpolated 3-h DSR values are also validated against the field measurements at the BSRN stations. R 2 of the 3-h data is 0.934, higher than that of the instantaneous values ( Figure 5). The bias and RMSE are −8.3 and 73.7 W/m 2 , respectively, smaller than those of their instantaneous counterparts. However, the relative RMSE of the 3-h DSR is 38.6%, much greater than the relative RMSE of the instantaneous DSR (19.4%). The higher relative error of 3-h data is not surprising, because the instantaneous data are directly retrieved from the MODIS observations while the temporal interpolation procedure of generating 3-h data will introduce uncertainties especially for the cases where atmospheric conditions are highly dynamic.
MCD18A1 does not directly provide the daily DSR value. However, the eight 3-h MCD18A1 DSR values during a day can be conveniently averaged to obtain daily mean DSR. While it is a simple approach, such an averaging method fails to account for nonlinearity in DSR time series data and can introduce some levels of uncertainties in estimating daily DSR. Using the field measurements, we calculated the daily DSR values by averaging the eight 3-h BSRN measurements and compared it with the daily DSR reference value, which uses all the 1-min BSRN measurements ( Figure 6). It shows the simple averaging approach introduces an error of 12.8 W/m 2 over the BSRN sites.   The comparison with in situ data over all the stations shows the MODIS daily mean DSR has a negative bias of −7.8 W/m 2 and RMSE of 32.3 W/m 2 (Figure 7), which represents a relative RMSE of 17.0%. The relative error of the daily DSR is slightly smaller than that of the instantaneous values and much smaller than that of the 3-h interpolated values. That indicates the temporal averaging process and the use of multiple observations of one day reduce the random noises of estimating DSR and result in higher accuracy of calculating daily mean DSR.
At the individual stations, the performance of the daily DSR is highly correlated with the quality of the instantaneous DSR (Table 3). The stations at high latitudes or altitudes with larger errors in instantaneous DSR typically have larger relative RMSE in estimating daily DSR. Besides, GAN, a tropical station with latitude of 23.11 • N, has a relative RMSE of 43.0%. One main contributor to the difference in data quality between the instantaneous values and the daily values is the dynamic cloud coverage in the area. Another important factor is the number of MODIS overpasses available in one day. Due to the polar-orbit characteristics, MODIS samples the lower latitudes with reduced frequency.
To assess the dependency of product performance on time of year, the validation results over Northern hemisphere sites are categorized by months ( Table 4). The MODIS DSR product shows good agreement with the field measurements for all the seasons. Winter (December, January, and February) finds the smallest absolute RMSE but the highest relative error. This high relative RMSE in winter can be attributed to difficulty in distinguishing snow-covered surfaces with clouds and high solar zenith angles in this season.

Figure 7.
Validation results of the daily MCD18A1 DSR product using the BSRN 2018 DSR measurements. To further evaluate the impacts of daily overpass counts on the accuracy of estimating daily DSR, the validation results are grouped by their number of MODIS overpasses within a day. The overpass counts have little impact on the bias of daily DSR. However, the importance of multiple observations in terms of reducing estimation errors is significant (Figure 8). If there is only one MODIS overpass, RMSE of estimating daily DSR is as high as 35.4 W/m 2 . RMSE drops to 33.0 W/m 2 , when the number of overpasses increases to two. RMSE can further decrease to 25.8 W/m 2 if the daily overpass count reaches four. Not surprisingly, the R 2 value increases with the number of daily MODIS overpasses. This finding is in line with our previous studies on temporal scaling of daily surface shortwave net radiation. As the overpass time of MODIS typically is close to noon when solar radiation reaches the maximum, two MODIS observations a day can generate daily estimates with rather high accuracy [24].

Comparison with Other DSR Products
Two other global DSR products (Clouds and the Earth's Radiant Energy System (CERES) and Global LAnd Surface Satellite (GLASS)) were also validated using the same data set of ground measurements and the results were compared with the MCD18A1 product. The gridded CERES DSR product is available at 1 • spatial resolution and several levels of temporal resolution [25]. The latest GLASS DSR product has a spatial resolution of 5 km and is only available at the daily temporal resolution [26]. The validation results of daily CERES and GLASS products are shown in Figure 9. CERES has the highest accuracy with RMSE of 28.8 W/m 2 and R 2 of 0.920. GLASS has a larger negative bias of 7.6 W/m 2 and greater RMSE of 35.9 W/m 2 . The accuracy of the MCD18A1 product is between them, with RMSE 3.5 W/m 2 higher than CERES and 3.6 W/m 2 lower than GLASS. However, it needs to be noted that the CERES product incorporates global geostationary satellite data better capturing the diurnal DSR variability with a much coarser spatial resolution, while MCD18 and GLASS radiation products use only the polar-orbiting MODIS observations.

Validation of PAR Product
The MCD18A2 PAR product was validated with one year data of PAR measurements at 25 AmeriFlux stations. PAR and DSR values are highly correlated. The earlier conclusions on DSR about the effects of surface reflectance and daily overpass count are also applicable to PAR. The PAR validation results are only summarized in brief here ( Figure 10). The MODIS instantaneous PAR product shows good agreement with the measured PAR with R 2 of 0.881, slightly higher than that of the DSR data. RMSE of PAR data is 45.7 W/m 2 , only half of that of the DSR data. Considering the absolute values of PAR and DSR, the two products have similar relative errors. Compared to DSR, PAR has a much smaller bias. The negative bias of DSR is mitigated mainly by the slight overestimation of PAR for high value part (clear-sky cases with small solar zenith angle). RMSE of daily PAR is 13.1 W/m 2 or 14.8% in a relative value, smaller than 17.6% that of the daily PAR product.

Summary and Future Plans
A new global product suite (MCD18) of land surface incident shortwave irradiance and its visible part PAR has been generated from MODIS data using NASA's MODIS Adaptive Processing System (MODAPS). The data products are archived and distributed through NASA's Land Processes Distributed Active Archive Center (LP DAAC). The current Collection 6.0 MCD18 products are produced at 5 km resolution and the forthcoming Collection 6.1 data sets will be available at the improved resolution of 1 km. The instantaneous values at the MODIS overpass time as well as the interpolated 3-h data of DSR and PAR are stacked in daily MODIS sinusoidal tile files. The products were validated with one year of ground measurements of surface shortwave irradiance collected at 33 BSRN and 25 AmeriFlux stations. The daily DSR has a bias of −7.8 W/m 2 and RMSE of 32.3 W/m 2 , while the daily averaged PAR has a bias of −1.5 W/m 2 and RMSE of 13.1 W/m 2 .
At the spatial resolution of km, the 3D effects of cloud coverage will cause issues in product validation because of the inconsistency between space-based satellite view geometry and perspective of ground-based pyranometers. As suggested by [2], temporal averaging of the field measurements and spatial averaging of the satellite products is a common approach to mitigate such 3D mismatch. The validation in the study employed both strategies. The results show RMSE of instantaneous DSR is reduced from 119.2 to 101.0 W/m 2 when the data was spatially aggregated to 50 km. Future study will include a comprehensive assessment of the impacts of 3D cloud effects on products with various spatial and temporal resolutions in different seasons and locations.
The current Collection 6.0 MCD18 products are the first release of the MODIS-based operational global high resolution DSR and PAR products. They were assessed in this study against a collection of in situ data with limited data length, climate type, and surface type variations. A comprehensive evaluation with more representative locations and longer time periods which can account for intra-annual variability and long-term trend is needed to provide users with more complete information on the product accuracy and meanwhile supply the product developers with diagnosis information which is essential for algorithm tuning and refining.
MODIS data onboard both the Terra and Aqua satellites are currently used in the MCD18 data production. While combining the twin sensors enables a global coverage at least once a day over most of the Earth's surface, significant coverage gaps remain unmapped at the low latitudes. Therefore, the daily global maps display data gaps. Our validation suggests the number of MODIS overpasses affects the accuracy of daily radiation estimation. The addition of geostationary data is crucial in filling the data gaps as well as enhancing the retrieval quality, particularly in the low latitudes where the MODIS overpass count per day is limited. The LUT approach used in the MCD18 generation is flexible in terms of ingesting multiple sources of satellite data. Additional satellite data will be used in the same manner as combining MODIS twin sensors. Each individual satellite observation, no matter on which satellite it is derived, will be used as an independent source of information to estimate separate instantaneous atmospheric optical properties. This represents a long range development plan to combine geostationary and polar-obit satellite data to deliver improved operational global shortwave irradiance products.
The LUT approach is used in the current MCD18 operational data production because of its robustness, maturity, and efficiency. However, the LUT approach also has disadvantages. The major drawback of the LUT approach is that it fails to take advantage of multispectral information of the MODIS sensor. As a result, its performance is compromised when the contrast between TOA reflectance and surface reflectance is limited. It is ideal when used to combine the LUT approach with other DSR estimation approaches to take advantage of their respective strengths. The next-generation DSR and PAR products should be an integrated data set, which may utilize multiple sources of satellite data as well as different retrieval algorithms.