Satellite-Based Estimation of Daily Average Net Radiation under Clear-Sky Conditions

Net radiation (R _n) is the balance of downward and upward short- and longwave radiation at the Earth's surface, and drives the process of evaporation, photosynthesis, and heating of air and soil (Bisht and Bras, 2010). It is a critical variable for estimation of surface energy budget and is used for weather prediction, climate monitoring, and agricultural meteorology (Bisht et al., 2005). DANR is referred to as the average accumulated net radiation during local time with a positive R _n value. Compared with instantaneous net radiation (INR), DANR or the diurnal cycle of net radiation has more applications, especially in models attempting to estimate evapotranspiration from remotely sensed data (Bastiaanssen, 2000; Jacobs et al., 2000; Bisht et al., 2005; Allen et al., 2006; Gao et al., 2008; Long et al., 2010).

So far, many efforts have been made to estimate net radiation under clear sky conditions, including those that estimate INR (Bastiaanssen et al., 1998; Roerink et al., 2000; Jiang and Islam, 2001; Su, 2002; Ma et al., 2002; Nishida et al., 2003; Norman et al., 2003; Bisht et al., 2005; Sobrino et al., 2007; Ryu et al., 2008), and daily net radiation (Dong et al., 1992; Hurtado and Sobrino, 2001; Irmak et al., 2003; Kjaersgaard et al., 2007; Samani et al., 2007; Wang and Liang, 2009a). However, compared with INR, there have been relatively few studies of DANR. According to the data used, the methodologies to obtain DANR can be classified into three categories: (2) ground-only measurements; (3) a combination of remotely sensed data and ground measurements; and (4) remotely sensed data only. DANR can be obtained directly from weather stations or field measurements (Jegede, 1997). However, the limitation of weather stations and field data is that they can only provide point observations of DANR, or the components of DANR, which often hinders its applications at regional scales (Irmak et al., 2003; Su et al., 2005; Samani et al., 2007; Long et al., 2010). Field observations and remotely sensed data have been combined to estimate DANR in several studies (Allen et al., 1998; Bastiaanssen, 2000; Jacobs et al., 2000; Gao et al., 2008; Long et al., 2010); however, these methods have been dependent on ground measurements. Studies that have estimated DANR using remotely sensed data only are rare. (Bisht et al., 2005) first proposed a scheme to estimate DANR for clear-sky days using remotely sensed observations only, and one of the strengths of the scheme is the use of a sinusoidal model to directly retrieve DANR with the INR estimate obtained from the satellite. This avoids having to define each component of DANR. However, the DANR estimates are easily affected by the INR estimates.

In the present paper, we propose a method that could overcome the dependence of DANR estimates upon ground observations by only using remotely sensed data to capture the spatial distribution of DANR estimates over large heterogeneous areas under clear sky conditions. The study uses data from remote sensing platforms, MODIS, a key instrument onboard the Earth Observing System (EOS) Terra and Aqua satellites. Terra MODIS and Aqua MODIS view the entire Earth's surface every one to two days, acquiring the attributes of the atmosphere and surface in 36 spectral bands spanning the visible, near-infrared and thermal regions of the electromagnetic spectrum, thereby offering an opportunity to improve global monitoring of terrestrial ecosystems. All of the components of DANR are estimated using MODIS land products (land surface albedo, land surface temperature, and land surface emissivity) and atmospheric products (air temperature, water vapor, ozone burden, and aerosol optical depth).

A strength of the proposed methodology is its ability to provide a spatially consistent and distributed DANR map for large heterogeneous areas under clear sky conditions. The key objectives for the paper are to develop a method to estimate the four components of DANR using remotely sensed data only as model input, and eliminating the need for ground measurements, and to validate the method using ground measurements from the Songnen Plain in Northeast China.

The study area is located southwest of the Songnen Plain, in Northeast China, covering the area (43^#cod#x000b0;02'-46^#cod#x000b0;06'N, 122^#cod#x000b0;21'-125^#cod#x000b0;36'E) (represented by the red rectangle in Fig. 1). The elevation of the area ranges from 64 m to 488 m and the area is characterized by various land cover types classified by the 500-m MODIS Land Cover Type product (MCD12Q1). The ten flux sites used here are shown in Fig. 1. Outside the study area, the eight flux sites (green points in Fig. 1) are part of a coordinated enhanced observation project for arid and semi-arid regions in northern China and represent different vegetation/land cover types: temperate grassland, cropland and deciduous broadleaf forest. JZ (Jingzhou), YK (Yingke) and GT (Guantao) are cropland sites, planted with maize with regular irrigation. MY (Miyun) is a deciduous broadleaf forest site. For grasslands, AR (Arou) is a sub-alpine meadow site, YZ (Yuzhong) represents typical steppe, DS (Dongsu) represents desert steppe, and NM (Naiman) is a sandy grassland site. Intensive calibration and maintenance were carried out weeks before the coordinated enhanced observation period (July to September) to ensure instrument performance and data quality. The data over these three months were used to study the longwave component of DANR.

Figure 1.  Study region (red rectangle) and the location of Tongyu-crop (red point), Tongyu-grass (black point), and the other eight flux sites (green points).

In the study area, TYC (Tongyu-crop; red point in Fig. 1) is a cropland site planted with maize/sunflowers (44^#cod#x000b0;35'N, 122^#cod#x000b0;52'E) and TYG (Tongyu-grass; black point in Fig. 1) is a degraded meadow steppe site (44^#cod#x000b0;34'N, 122^#cod#x000b0;55'E) (Wang et al., 2010). Both TYC and TYG are reference sites for the CEOP (coordinated enhanced observing period), initially proposed in 1997 during the Global Energy and Water Cycle Experiment (GEWEX). This experiment focused on the measurement, understanding and modeling of water and energy cycles within the climate systems. Quality and error checking procedures were performed regularly to maintain continuous, high-quality measurements during the planned CEOP periods. At the TYC and TYG sites (between 2003 and 2004), the four components of net radiation were measured at a height of 4 m above ground using a Kipp and Zonen CNR-1 radiometer (Kipp and Zonen, The Netherlands), which acquired data at intervals of several seconds and averaged them every 15 minutes. The data for 30 clear sky days at TYC and 27 clear-sky days at TYG were used to validate the DANR estimates in the study.

Satellite daily time series were applied to estimate the four radiation components of DANR. The MODIS datasets used in this study included MOD11_L2 and MCD43B3 land data products (collection 5) provided by the MODIS Land Team (http://modis-land.gsfc.nasa.gov/). MOD03 geolocation data, MOD04, MOD05 and MOD/MYD07 atmosphere data products were archived by the MODIS Adaptive Processing System (MODAPS) in the Level 1 and Atmospheres Archive and Distribution System (http://ladsweb.nascom.nasa.gov/), which are shown in Table 1. The pixels for this study are classified as clear-sky pixels if a 1-km LST estimate within a MODIS overpass is available from MOD11_L2. To select the MODIS datasets for days with clear skies, we examined the fifth byte of the quality assurance bit field (i.e., the number of 1-km² clear pixels within a 25-km² pixel) in MOD/MYD07_L2, and selected those data for perfectly clear days (i.e., the number of clear pixels is 25). MOD03 and MCD43B3 products can be used for clear sky and cloudy sky conditions, while MOD04, MOD/MYD07 and MOD11_L2 can only be used for clear sky conditions. MOD05_L2 contains column water-vapor amounts over clear land areas of the globe and above clouds over both land and ocean. These products are assumed to be homogenous over the corresponding grid, and the MODIS overpass data used in this paper were re-projected to the geographic projection at a spatial resolution of 0.009^#cod#x000b0;. The dimensions of the study area were 344 rows by 373 columns with a pixel resolution of 1 km. Table 2 provides a description of the study period and satellite overpass time over the study region. The fifth and sixth column of Table 2 respectively contain the times of sunrise and sunset (local time), as obtained from the US Naval Observatory website, Astronomical Application Department (http://aa.usno.navy.mil/), which were used in the estimation of the diurnal cycle of downward shortwave radiation.

DANR can be expressed as

where R _n,q is the DANR (W m^-2), q is the period between about 1 h after local sunrise time and 1 h before local sunset time [these two values correspond to the local time at which R _n becomes positive and negative, respectively (Bisht and Bras, 2010)]; R _S,q^#cod#8595; and R _S,q^#cod#8593; are the q-h average downward and upward shortwave radiation (W m^-2); and R _L,q^#cod#8595; and R _L,q^#cod#8593; are the q-h average downward and upward longwave radiation (W m^-2).

3.1.1 Instantaneous downward shortwave radiation estimates

Various parameterizations for downward shortwave radiation have been presented in the scientific literature (Niemela et al., 2001), in which it shows that Iqbal's (1983) parameterization scheme outperforms the others. This scheme considers both diffuse and direct beam radiation. The instantaneous solar irradiance (R _s) can be expressed as

where R_dir and R _dif are direct and diffuse beam radiation \hbox(W m^-2), respectively; z is the solar zenith angle; I _sc is the solar constant (about 1367 W m^-2); #cod#964; _r is the transmittance of Rayleigh scattering of air molecules; #cod#964; _g is the transmittance of absorptance of uniformly mixed gases; #cod#964; _w is the transmittance of water vapor absorptance; #cod#964; _a is the transmittance of aerosol absorptance and scattering; #cod#964; _o is the transmittance of ozone absorptance; and R _dr, R _da and R _dm represent the broadband irradiance on the ground arising from Rayleigh scattering, aerosol scattering, and multiple reflection between the Earth's surface and atmosphere, respectively. Because there are several parameters in the estimations of #cod#964; _w and #cod#964; _a that are not readily available at the regional scale, we used three schemes——air mass at standard pressure (m _r); absorptance of water vapor (a _w); and transmittance of aerosol absorptance and scattering (#cod#964; _a)——to substitute for the original scheme (Appendices A-C). The specific schemes for other sub-terms can be found in (Iqbal, 1983).

3.1.2 Diurnal cycle and q-hour average of downward shortwave radiation

(Monteith, 1973) proposed a framework to obtain the diurnal cycle of solar irradiance during the day length (the time period between sunrise and sunset) under cloudless days by modeling it as a sinusoidal function, and the irradiance at t hours after sunrise can be expressed as

where R _s is instantaneous solar irradiance; R _{s, max} is the maximum irradiance at solar noon; and N is day length in hours.

(Jackson et al., 1983) demonstrated that the value of N obtained from the sinusoidal fit of solar irradiance in Eq. (6) was smaller than day length, and that these two values of sunrise and sunset time in the sinusoidal function of Eq. (6) virtually correspond to a time after the local sunrise time and a time before the local sunrise time [shown in Fig. 1 of (Jackson et al., 1983)]. We argue that the day length value in Eq. (6) can cause an increased diurnal cycle retrieval value of downward shortwave radiation, and thus we propose an improved sinusoidal model for estimating the diurnal cycle of R _s for clear-sky conditions, as expressed below:

Here, h is a time between h _rise and h _set, and these two values are approximated as 30 min after local sunrise time and 30 min before local sunset time——discussed in more detail in subsequent sections.

The diurnal cycle of downward shortwave radiation at the satellite overpass time (t _overpass) is shown in Fig. 2. For a given study day, R _s,max and R _S,q^#cod#8595; can be expressed as

and

where D=h _set-h _rise; and t _rise and t _set correspond to the local time at which R _n becomes positive and negative (Fig. 2), respectively. These two values were approximated as 1 h after local sunrise time and 1 h before local sunset time, respectively (Bisht and Bras, 2010).

Figure 2.  The improved sinusoidal model with MODIS overpass.

The q-h average upward shortwave radiation (R _S,q^#cod#8593;) can be estimated as

where #cod#945; is land surface albedo.

The MODIS BRDF/albedo algorithm combines a Ross-Thick-LiSparse-Reciprocal (RTLSR) kernel-driven BRDF model with multi-date, cloud-free and multi-spectral data to characterize the anisotropic reflectivity of the land surface (Ross, 1981; Li and Strahler, 1992; Roujean et al., 1992; Wanner et al., 1995; Schaaf et al., 2002). For the present paper, we selected the MCD43B3 product, a 16-day composite consisting of black- and white-sky albedos for seven spectral bands (band 1-7) and the three broadbands (0.30-0.7 #cod#956;m, 0.7-3.0 #cod#956;m and 0.3-5.0 #cod#956;m). Black-sky albedo (directional-hemispherical reflectance) is referred as albedo in the absence of a diffuse irradiance component and is a function of the solar zenith angle. White-sky albedo (bi-hemispherical reflectance) is albedo in the absence of a direct irradiance component when the diffuse component is isotropic. Black-sky albedo and white-sky albedo mark the extreme cases of completely direct and completely diffuse illumination. For a waveband (#cod#923;), actual albedo, #cod#945;(#cod#952;,#cod#923;), is a value interpolated between these two depending on the aerosol optical depth (#cod#915;). The equation to compute actual albedo was given by (Lucht et al., 2000) and is expressed as follows:

Here, #cod#945; _bs(#cod#952;,#cod#923;) is black-sky albedo; #cod#945; _WS(#cod#952;,#cod#923;) is white-sky albedo; s(#cod#952;,#cod#915;(#cod#923;)) is the fraction of diffuse sky light assumed to be isotropic (obtained by a look-up table available from the MODIS albedo products homepage); and shortwave albedo was estimated using seven MODIS spectral albedos and weights of converting surface MODIS narrowband albedos to broadband albedos (Liang et al., 1999).

Presently, R _L,q^#cod#8595; is not directly estimated due to limitations of the remotely sensed data: the mean atmospheric emissivity and the mean air temperature during a positive R _n time were not measured. The observations show that the 24-h curve of downward longwave radiation (R _L^#cod#8595;) is gentle compared to the other components of R _n. The observed R _L^#cod#8595; from 643 clear days were used to compare the daily average downward longwave radiation (R _L,24^#cod#8595;) with R _L,q^#cod#8595;, obtained from the 10 flux sites used in this study. The comparison between R _L,24^#cod#8595; and R _L,q^#cod#8595; is shown in Fig. 3, which demonstrates good agreement [bias, RMSE and correlation (r²) of 4.59 W m^-2, 8.21W m^-2 and 0.99, respectively]. This level of error is acceptable, and thus we use the value of R _L,24^#cod#8595; as a replacement for R _L,q^#cod#8595;, and can be estimated as: Here, #cod#949; _a is the air emissivity (Swinbank, 1963); #cod#963; is the Steffan-Boltzmann constant (5.67#cod#215; 10^-8 W m^-2 K^-4); T _a and are the air temperature (K) and daily average air temperature at screen level, respectively; and is the daily average air emissivity, assumed to be approximately equal to the value obtained from Eq. (13) using .

Figure 3.  Comparison between daily average downward longwave radiation and q-h average downward longwave radiation.

The Terra MODIS instrument acquires data twice daily at 1030UTC and 2230UTC, as does the Aqua MODIS (0130 UTC and 1330 UTC). For a given area, the four daily MODIS observations provide an opportunity to track the diurnal cycle of air temperature or LST, and are a great help in obtaining the or daily average LST at regional scales. The MOD/MYD07_L2 product provides the air and dew-point temperature at 20 different atmospheric pressure levels ranging from 5 to 1000 hPa, retrieved using a statistical regression algorithm (Seemann et al., 2003). (Tang and Li, 2008) suggested that air temperatures at screen level can be used to interpolate temperatures at 1000 hPa when assuming a hydrostatic atmosphere. These values at 1000 hPa were extrapolated to obtain temperature values at screen level using Eq. (8) given by (Bisht and Bras, 2010). (Long et al., 2010) estimated the daily average LST by fitting a three-order polynomial function with four observations of LST obtained from the MOD/MYD11A1 data product. In the present paper, we fit a three-order polynomial function using four observations of air temperature from MOD/MYD07_L2 data to estimate . The fitting of a three-order polynomial function can be expressed as where b₀, b₁, b₂ and b₃ are regression coefficients. According to the four observations from MOD/MOD07_L2 and four acquisition times, the value can be obtained after the regression coefficients are computed.

The R _L,q^#cod#8593; from the land surface can be given as

where #cod#949; _s is the land surface emissivity and is the q-h average land surface temperature (K).

The MOD11_L2 product provided the daily LST values retrieved using the generalized split-window algorithm (Wan and Dozier, 1996), and emissivity in bands 31 and 32 estimated using the classification-based emissivity method (Snyder and Wan, 1998). A study to estimate from remote sensing sensors has not been reported in the literature thus far. In the present study, we used one daytime observation of LST from the MOD11_L2 as a surrogate for . The #cod#949; _s value was computed using the following nonlinear formula (Liang, 2004).

where #cod#949;₃₁ and #cod#949;₃₂ denote the emissivity in MODIS bands 31 and 32, respectively.

The air temperature values obtained from the MOD/ MYD07_L2 product at the 1000 hPa level were underestimated, with a bias, RMSE and r² of -1.13^#cod#x000b0;C, 3.21^#cod#x000b0;C and 0.89 at the TYC site, and -1.33^#cod#x000b0;C, 3.08^#cod#x000b0;C and 0.90 at the TYG site, respectively, (Fig. 4a); this was perhaps created by the differences in screen level and 1000 hPa height. When these temperature values were extrapolated to obtain values at screen level assuming a hydrostatic atmosphere, these extrapolated temperature values at both sites had slightly larger errors, with biases and RMSEs of -1.23^#cod#x000b0;C and 3.31^#cod#x000b0;C, and -1.45^#cod#x000b0;C and 3.20^#cod#x000b0;C, respectively. This could be due to the air pressure being lower than 1000 hPa at the study sites. Air temperature values obtained using Eq. (8) by (Bisht and Bras, 2010) will be lower than MOD/MYD07 temperature values at the 1000 hPa level. In addition, the decision to perform the MOD/MYD07 atmospheric and surface parameter retrievals depended upon the validity of the MODIS cloud mask algorithm (Ackerman et al., 1998), which determined if a given pixel was clear by combining the results of several spectral threshold tests. It required that at least five of the 25 pixels in a 5#cod#215; 5 field-of-view area be assigned a 95% or greater confidence of clear by the cloud mask. Cloud contamination occurred if the cloud mask failed to detect a cloud; this affected the retrieval accuracies of MOD/MYD07 air temperature, and caused the differences between air temperature measured at screen level and that derived by MODIS, which can cause an error in the downward longwave component of DANR.

Figure 4.  Comparison of (a) instantaneous air temperature obtained from the MOD/MYD07_L2 product at the 1000 hpa level and the measured air temperature, and (b) the estimated daily average air temperature and the daily average accumulated measurement.

Compared to the MODIS-derived instantaneous air temperature, the daily average air temperature demonstrated a stronger correlation with smaller errors: the bias, RMSE and r² for were -0.18^#cod#x000b0;C, 1.85^#cod#x000b0;C and 0.95 at the TYC site, and 0.17^#cod#x000b0;C, 2.01^#cod#x000b0;C and 0.94 at the TYG site, respectively (Fig. 4b). Sensitivity analysis demonstrated that errors up to #cod#177;2 K in resulted in errors of about #cod#177; 5-#cod#177; 12 W m^-2 in the estimation of R _L,24^#cod#8595;, within a daily average air temperature range of 240-315 K.

Direct measurements of LST were not available, and thus observations of the upward longwave radiation flux were converted to obtain surrogate measurements of LST using the Steffan-Boltzmann equation. Surface emissivity was assumed to be equal to the corresponding #cod#949; _s value obtained using Eq. (16). The average accumulated LST during the time with a positive R _n was used as the observed value of . The bias, RMSE and r² for were 1.44 K, 2.68 K and 0.93 at the TYC site, and 0.63 K, 2.46 K and 0.92 at the TYG site (not shown), respectively. Sensitivity analysis demonstrated that errors up to #cod#177;3 K in resulted in errors of about #cod#177;14-#cod#177;22 W m^-2 in the estimation of R#cod#8593; , within a q-h mean LST range of 275-320 K.

4.2.1 Instantaneous downward shortwave radiation estimates

The statistical information about all radiation flux is summarized in Table 3. The comparison between the observed and modeled R _s is shown in Fig. 5a, revealing a bias and RMSE of 25 W m^-2 and 31 W m^-2 at the TYC site, and 22 W m^-2 and 28 W m^-2 at the TYG site, respectively. The mean relative error in the estimation of R _s at both sites was approximately of the order of 4% of the ground observations. The RMSE value obtained from the modified (Iqbal, 1983) parameterization scheme was comparable to the value reported by (Houborg and Soegaard, 2004), which used Iqbal's (1983) scheme and obtained an overall RMSE of 47.7 W m^-2. Using Zillman's (1972) scheme, which used the screen-level vapor pressure and solar zenith angle as input variables, (Bisht et al., 2005) reported the bias, RMSE and r² for shortwave radiation downward as 41 W m^-2, 51 W m^-2 and 0.97, respectively. This suggests that the modified (Iqbal, 1983) scheme can be easily used to estimate the downward shortwave radiation at the regional scale with reasonable estimation accuracy.

Figure 5.  Comparison at the cropland and grassland sites of estimated and measured (a) instantaneous downward shortwave radiation; (b, c) downward and upward q-h average shortwave radiation; (d, e) downward and upward q-h average longwave radiation; and (e) daily average net radiation.

Solar irradiance is attenuated during its passage through the atmosphere. Water vapor, aerosols, and the scattering of air molecules are three main attenuators under clear-sky conditions (Iqbal, 1983) that affect the atmospheric transmittance. During study days, minimum attenuation was affected by the uniformly mixed gases with a transmittance range of 0.98-0.99, while ozone attenuation was slightly greater, ranging from 0.97-0.99. The transmittance of water vapor ranged from 0.85-0.93; aerosols obtained by the replaced scheme ranged from 0.49-0.99; and Rayleigh scattering ranged from 0.82-0.92. This indicated that the transmission of water vapor, aerosols, and Rayleigh scattering are three main factors affecting atmospheric transmittance, calculated as the product of above five attenuating atmospheric constituents. (Liu and Jordan, 1960) stated that atmospheric transmittance for clear-sky days is between 0.45 and 0.75, while (Gates, 1980) stated a range of between 0.60 and 0.70, and for extremely clear-sky conditions reported atmospheric transmittance reaching around 0.75. In the present study, atmospheric transmittance was measured between 0.46 and 0.77, with 13 values greater than 0.75 and 32 values greater than 0.70 at both sites. This suggests that several atmospheric transmittance values were overestimated, resulting in the overestimation of R _s.

In addition, during conditions of very low aerosol loading presence for clear-sky days, retrieval accuracy of aerosol optical depth (AOD) from MODIS affects estimation of downward shortwave radiation. The 31 values from the MODIS-derived AOD at both sites were between -0.05 and 0.05. Levy et al. (2009) stated that MODIS is not sensitive enough over land to retrieve aerosol to better than #cod#177;0.05, which means that in very clean conditions the algorithm cannot determine #cod#915;=0.00 from 0.05 or -0.05. Thus, a value in this range may be over- or underestimated; if the value is underestimated, it can cause an overestimation of aerosol transmittance values, contributing to the overestimation of atmospheric transmittance, and ultimately causing an overestimation of R _s.

4.2.2 Diurnal cycle of downward shortwave radiation estimates

The diurnal cycle of downward shortwave radiation estimated using Monteith's (1973) sinusoidal model (green line with triangles) and the improved model (red line) are shown in Fig. 6. The first two cases demonstrate the overestimation and underestimation of the instantaneous R _s value, which shows that the error in R _s affects the retrieval of the diurnal cycle. The third case shows that the estimation of h _rise and h _set as 30 min after sunrise and before sunset might not be accurate, and highlights the joint impact of the sinusoidal model time period, and of the error in R _s on the diurnal cycle retrieval. The fourth case shows that the diurnal cycle was retrieved accurately. In all the cases, the area enclosed by the green sinusoidal curve is larger than the red sinusoidal curve. Compared to the diurnal cycle of observations, Monteith's (1973) sinusoidal model better matches with the observations during a time period between about 1 h before sunset time and sunset time. However, the red line fits better than the green line during the period between the time of sunrise and about 1 h before sunset time, as shown in Figs. 6a, b and d. Overall, the improved sinusoidal model better captures the diurnal variation in downward shortwave radiation for days with clear sky.

The diurnal cycle retrieval of downward shortwave radiation is affected by the time period (D) between h _rise and h _set, as is the R#cod#8595; . Because h _rise and h_\mboxset are not directly observed, an accurate time period was not obtained, but one can be calculated approximately using Eq. (8), as described by (Jackson et al., 1983). The TYC and TYG sites are roughly located at the same latitude, and the values of D calculated by (Jackson et al., 1983) were the same for the same dates. Thus, we used one TYC site as a case study. The estimation of h _rise and h _set as 20, 25, 30 or 35 min after sunrise and before sunset produced the different time periods (D₁, D₂, D₃ and D₄), which affected the R _S,q^#cod#8595; value (Table 4). For the given study period, the difference between D and D₁ was a negative value every day, and the mean bias when estimating the R _S,q^#cod#8595; value was 35.73 W m^-2, while the difference between D and D₄ was a positive value, and the mean bias was 19.56 W m^-2. Similarly, the mean biases obtained using D₂ and D₃ were 30.55 W m^-2 and 25.16 W m^-2, respectively. This indicates that the estimation of h _rise and h _set as 20 or 35 min after sunrise and before sunset resulted in the overestimation and/or underestimation of the sinusoidal model time period and was likely to produce a systematic bias. Both values of D₂ and D₃ approach the value of D. The RMSE and bias between D and D₂ are 8 min and -6 min respectively, while the RMSE and bias between D and D₃ are 6 min and 4 min respectively, indicating that D₃ is closer to D than D₂ . Thus, the estimations of h _rise and h _set as 30 min after sunrise and before sunset are rational and reliable, and the effects of these estimates on the R _S,q^#cod#8595; value are negligible.

Figure 6.  Diurnal cycle of downward shortwave radiation when instantaneous R _s was (a) overestimated and (b) underestimated. (c) Diurnal cycle of downward shortwave radiation, showing the joint impact of the time period of the sinusoidal model and the overestimated R _s. (d) Diurnal cycle of downward shortwave radiation retrieved accurately during the Terra-MODIS overpass.

4.2.3 q-hour average downward shortwave radiation estimates

Compared to the instantaneous R _s value, the R _S,q^#cod#8595; value had a lower magnitude of error, with a bias and RMSE of 25 W m^-2 and 28 W m^-2 at the TYC site, and 21 W m^-2 and 27 W m^-2 at the TYG site, respectively, (Fig. 5b). The estimated R _s value was substituted for the observed value, and the bias, RMSE and r² for R _S,q^#cod#8595; were 1 W m^-2, 10 \hboxW m^-2 and 0.99, respectively. This indicated that an overestimation of R _S,q^#cod#8595; was mainly due to an overestimation of the R _s. The decreasing trend for R _S,q^#cod#8595; values from summer to autumn in 2003 and 2004 was almost the same as the observed at both sites, and trends of these values decreased with dates (Figs. 7a and b). This seasonal change of R _S,q^#cod#8595; values could be closely related to the variation of the solar zenith angle, caused mainly by variation in solar declination. However, this decreasing trend of R _S,q^#cod#8595; values on adjacent days of the year was not obvious, due to the different extinction of solar radiation and the variation of the solar altitude angle at satellite overpass time. It can be concluded that the improved sinusoidal model is able to reliably retrieve R _S,q^#cod#8595; from the instantaneous R _s value under clear-sky conditions, using remotely sensed data only, at a regional scale.

Figure 7.  Trend comparison of the estimated and observed daily average net radiation with four components at the (a) TYC and (b) TYG sites.

The R _S,q^#cod#8593; value was underestimated with a bias and RMSE of -10 W m^-2 and 13 W m^-2 at the TYC site, and -12 \hboxW m^-2 and 14 W m^-2 at the TYG site, respectively (Fig. 5c). The MCD43B3 albedo product in this study was underestimated with a bias, RMSE and r² of -0.02, 0.027 and 0.31 at the TYC site, and -0.027, 0.031 and 0.32 at the TYG site, respectively, as compared to the observed albedo, which was obtained by the ratio of the upward to downward shortwave radiation flux using field measurements. The bias and RMSE between the observed albedo at the TYC site and that at the TYG site were -0.023 and 0.032 respectively, indicating that there was an obvious difference between the two surface albedos during the study days. This difference affected the R _S,q^#cod#8593; value (Figs. 7a and b). When the observed albedo value was used instead of the MODIS-derived albedo, the bias, RMSE and r² for R _S,q^#cod#8593; were -1 W m^-2, 4 W m^-2, and 0.93 at the TYC site, and 1 W m^-2, 4 W m^-2 and 0.94 at the TYG site, respectively. This indicated that the underestimated R _S,q^#cod#8593; value was caused mainly by the underestimation of the albedo.

The upward shortwave radiation was affected by land cover heterogeneity (Jin et al., 2003; Salomon et al., 2006; Ryu et al., 2008) and temporal heterogeneity (Ryu et al., 2008). At the TYC site, the source area for the upward radiation measurements from the beginning of August to the end of September in 2003 and 2004 was mainly covered by maize, and the land surface of the MODIS 1-km pixel from the surrounding site was covered by a mixture of crops, bare soil, and trees. During this period, between the MODIS-derived albedo and the observed albedo, the RMSE and bias were 0.024 and -0.02, respectively. When the source area was completely covered with bare soil, the RMSE and bias were 0.034 and -0.03, respectively, from the beginning of October to the end of November. At the TYG site, the RMSE and bias were 0.028 and -0.026 from the beginning of August to the end of September in 2003 and 2004, and 0.035 and -0.031 from October to November, respectively. A change in color of the grass during this time occurred, suggesting that the albedo-induced underestimation of R _S,q^#cod#8593; was partly caused by land cover heterogeneity. It was assumed that the vegetation cover did not change during this 16-d period; if change occurred, the MODIS-derived surface albedo estimates likely caused error in net radiation estimates (Bisht et al., 2005), and created error in R _S,q^#cod#8593; values. In addition, the large discrepancy of R _S,q^#cod#8593; values at both the sites may have been affected by temporal heterogeneity. The RMSE and bias between the MODIS 16-d albedo product and q-h averaged ground albedo were 0.038 and -0.032 at the TYC site, and 0.039 and -0.035 at the TYG site, respectively, showing that the 16-d composite MODIS albedo was underestimated, and that the mismatch of temporal scales between MODIS-derived albedo and the observed albedo contributed to the error in R _S,q^#cod#8593; values.

The R _L,q^#cod#8595; value obtained from the R _L,24^#cod#8595; value was slightly overestimated, with a bias and RMSE of 4 W m^-2 and 12 W m^-2 at the TYC site, and 2 W m^-2 and 12 W m^-2 at the TYG site, respectively (Fig. 5d), probably caused by the slight overestimation of . The bias, RMSE and r² for were 0.011, 0.025 and 0.63 at the TYC site, and 0.010, 0.025 and 0.66 at the TYG site, respectively. These values were calculated from the average accumulated atmospheric emissivity during a 24-h period, and measurements of downward longwave radiation and air temperature at intervals of 15 min that were converted to obtain surrogate observations of atmospheric emissivity. The error range of the various methods was about 10-32 W m^-2 for instantaneous clear-sky R#cod#8595; values (Ellingson, 1995; Zhou et al., 2007; Tang and Li, 2008; Ryu et al., 2008; Wang and Liang, 2009b). Although the above mentioned error was obtained from the instantaneous values, it can also be regarded as a reference for studying R#cod#8595; . It suggested that the R _L,24^#cod#8595; value obtained from Eqs. (12) and (13) approached the observed R _L,q^#cod#8595;, and thus using the R _L,24^#cod#8595; value as a replacement for R _L,q^#cod#8595; was a feasible method.

Estimation of R _L,q^#cod#8595; is a challenge at regional scales due to complex atmospheric constituents and the limitations of remote sensing sensors. It was difficult to acquire high spatial and temporal resolution information for air temperature and atmospheric emissivity, which control the q-h downward longwave radiation. Thus, we could not directly obtain the R _L,q^#cod#8595; value using remote sensing data. In the present study, we used the R _L,24^#cod#8595; value as an approximate surrogate for the R _L,q^#cod#8595; value. For clear-sky conditions, the difference between R _L,24^#cod#8595; and R _L,q^#cod#8595; was affected by elevation and topography. The RMSE value between R _L,24^#cod#8595; and R _L,q^#cod#8595; at the AR site (located in the northeastern Tibetan Plateau) was the largest among all the sites (Table 5), followed by the YK site (located in the artificial oasis to the south of Zhangye City (38^#cod#x000b0;51'N, 100^#cod#x000b0;24'E), where the main crop is maize with row structure and regular irrigation. It has been argued that the daily temperature range is proportional to the elevation over large-scale terrains (i.e., the higher elevation, the larger the daily range in air temperature) (Hu and Chen, 2006). This can cause a large discrepancy between R _L,24^#cod#8595; and R _L,q^#cod#8595;. Compared to the YK site, the YZ site (located at the top of Cuiying Mountain) had a lower RMSE value due to the smaller daily range of air temperature.

The R _L,q^#cod#8593; value was overestimated at both the sites with a bias and RMSE of 8 W m^-2 and 16 W m^-2 at the TYC site, and 3 W m^-2 and 15 W m^-2 at the TYG site, respectively (Fig. 5e), caused mainly by an overestimation of . In this paper, we used the LST in the MOD11_L2 acquired around 1030 UTC for Terra-MODIS as a surrogate estimate of . The difference in values was possibly due to the heterogeneity of the land cover/use, as was the R _L,q^#cod#8593; value. The error statistics of different land cover types for the value and the R _L,q^#cod#8593; value from July to September 2008 are shown in Table 6. The overall RMSE and mean bias of R _L,q^#cod#8593; and from 10 sites were 18.54 W m^-2 and 8.08 W m^-2, and 3.02^#cod#x000b0;C and 1.28^#cod#x000b0;C, respectively. The RMSE of at the AR site was the largest among all the sites, perhaps due to the complex diurnal variations of LST for large-scale terrain areas. The RMSE of at the JZ site was the smallest, perhaps due to the homogenous land cover (maize) within the MODIS LST extent (1 km). The sites that were densely vegetated during the growing season had a relatively smaller RMSE for the values and the R _L,q^#cod#8593; values, as compared with those of the sites that were sparsely vegetated. The values of were overestimated and resulted in the overestimation of the R _L,q^#cod#8593; value, which could have been due to the fact that daytime land surface temperature in MOD11_L2 does not truly represent the .

We acknowledge that the biases computed from the nine sites only may not have been appropriate for other land cover types, such as urban areas and desert, and did not involve other factors that affected the diurnal variations of LST and may have caused large errors at certain sites: surface air temperature, humidity, wind speed and soil moisture. Various methods were used to estimate clear-sky upward longwave radiation, with errors ranging from about 12-31 W m^-2 for instantaneous upward longwave radiation (Ryu et al., 2008; Wu et al., 2012). This suggests that the error in R _L,q^#cod#8593; over the densely vegetated areas under clear-sky conditions using daytime LST values in the MOD11_L2 product as a replacement for the value [using Eq. (15)] was acceptable, but the method may cause large errors over sparsely vegetated areas.

DANR was estimated using only the Terra-MODIS and Aqua-MODIS data, and the spatial datasets of DANR were obtained using the scheme discussed in section 3. An example of the DANR map for the study domain on 14 August 2004 is shown in Fig. 8. A comparison between the observed and simulated DANR values indicated that there was overestimation of values at both sites (Fig. 5f). This overestimation was likely due to the overestimation of the net shortwave component of DANR, which in turn was caused by the overestimation of the R _S,q^#cod#8595; value and underestimation of the R _S,q^#cod#8593; value. The error in DANR was approximately of the order of 12% and 14% of the ground observations at the TYC and TYG sites, respectively. The trends of the observed and the modeled DANR values for clear-sky days from summer to autumn were generally in agreement at both sites (Figs. 7a and b). (Bisht et al., 2005) performed a study of estimating net radiation using only Terra-MODIS data and they reported that the bias, RMSE and r² for DANR were 50 W m^-2, 60 \hboxW m^-2 and 0.85, while our results revealed a bias, RMSE and r² of 31 W m^-2, 37 \hboxW m^-2 and 0.95 at the TYC site, and 33 \hboxW m^-2, 40 W m^-2 and 0.93 at the TYG site, respectively. When the observed albedo value was used instead of the MODIS-derived albedo, the bias, RMSE and r² for DANR were 21 W m^-2, 28 W m^-2 and 0.96 at the TYC site, and 20 W m^-2, 30 W m^-2 and 0.94 at the TYG site, respectively. This suggests that more accurate information about land surface albedo can improve the proposed methodology.

Figure 8.  A typical daily average net radiation map for 14 Aug 2004. The white region on the map refers to where no information was available.

In the present study, we developed an approach to estimate daily average net radiation with four radiation components at a regional scale using only remotely sensed data for clear-sky conditions. The method attempted to overcome the need for ground observations as model input by exploiting the land and atmospheric data products from Terra-MODIS and Aqua-MODIS. The estimates were validated by in-situ data at cropland and grassland sites located in the Songnen Plain of Northeast China during a period of 30 clear-sky days during 2003-04. The main results of the study can be summarized as follows:

(2) Although Iqbal's (1983) scheme has definite physical meaning with reasonable accuracy, it requires atmospheric transmittance information about Rayleigh scattering, uniformly mixed gases, water vapor, aerosols and ozone. Such information is not readily available, and thus limits its application at regional scales. By substituting for the schemes of air mass at standard pressure, aerosols and water vapor, we easily estimated the downward shortwave radiation at the regional scale with a small level of error.

(3) Compared to Monteith's (1973) sinusoidal model, the improved model, using two values of 30 min after local sunrise time and 30 min before local sunset time as model input, was able to better capture the diurnal variation of downward shortwave radiation from a single instantaneous downward shortwave radiation value obtained from the satellite.

(4) The comparison of errors in R _S,q^#cod#8595; obtained using downward shortwave radiation estimates and ground measurements suggested that R _S,q^#cod#8595; was overestimated mainly due to the overestimation of downward shortwave radiation.

(5) The R _S,q^#cod#8593; value was underestimated due to the underestimation of the albedo. The MCD43B3 albedo value for both sites was underestimated compared to ground observations, probably caused by land cover heterogeneity and a mismatch of the temporal scale between the q-hour averaged ground albedo and 16-day composite MODIS albedo.

(6) The comparison between the R _L,q^#cod#8595; value obtained from the R _L,24^#cod#8595; value and the observed R _L,q^#cod#8595; value suggested that, under conditions of R _L,q^#cod#8595; not being available, it is possible to use remote sensing data to obtain an R _L,24^#cod#8595; value and use it as a replacement for the R _L,q^#cod#8595; value.

(7) The R _L,q^#cod#8593; value was overestimated mainly due to the overestimation of the value. Comparing the errors of the R _L,q^#cod#8593; from those sites suggested that the error in the R _L,q^#cod#8593; value for the densely vegetated areas was acceptable by using daytime LST in MOD11_L2 as a replacement for the value, but the method may cause large errors over sparsely vegetated areas.

(8) The DANR was estimated using remotely sensed data only, and the results showed an overestimation at both sites, largely due to the overestimation of the R _S,q^#cod#8595; value and underestimation of the R _S,q^#cod#8593; value. The error in the DANR estimate from this study was comparable to that obtained from the current method, which also used remotely sensed data.

Although our proposed method does not involve a complex radiative transfer model, it can estimate regional DANR together with four components for clear-sky conditions using only remotely sensed data (without the need for ground observations). Our satellite-derived DANR datasets are spatially consistent and reasonably distributed, and could be used to estimate evapotranspiration, which is an important parameter used in regional climate models.

In Iqbal's (1983) scheme, the air mass at actual pressure can be expressed as

and

where m_a (dimensionless) is the air mass at actual pressure; p (hPa) is the actual pressure; p₀ is the standard atmospheric pressure (1013.25 hPa); and m _r is computed according to (Kasten and Young, 1989).

Absorptance of water vapor is computed by correlation according to (Lacis and Hansen, 1974):

where u₁ is the pressure-corrected relative optical path length of precipitable water; u₁=wm _a; w is the water vapor amount obtained from the MODIS near-infrared total precipitable water product (MOD 05); and 1-km near-infrared daylight only is used.

The transmittance of aerosol absorptance and scattering is computed by an equation from the REST model:

where #cod#915; is aerosol optical depth, which can be obtained from the MODIS aerosol product (MOD04) and that at 0.550 #cod#956;m is used for the computation of transmittance.