An Efficient Method of Estimating Downward Solar Radiation Based on the MODIS Observations for the Use of Land Surface Modeling

Solar radiation is a critical variable in global change sciences. While most of the current global datasets provide only the total downward solar radiation, we aim to develop a method to estimate the downward global land surface solar radiation and its partitioned direct and diffuse components, which provide the necessary key meteorological inputs for most land surface models. We developed a simple satellite-based computing scheme to enable fast and reliable estimation of these variables. The global Moderate Resolution Imaging Spectroradiometer (MODIS) products at 1° spatial resolution for the period 2003–2011 were used as the forcing data. Evaluations at Baseline Surface Radiation Network (BSRN) sites show good agreement between the estimated radiation and ground-based observations. At all the 48 BSRN sites, the RMSE between the observations and estimations are 34.59, 41.98 and 28.06 W·m for total, direct and diffuse solar radiation, respectively. Our estimations tend to slightly overestimate the total and diffuse but underestimate the direct solar radiation. The errors may be related to the simple model structure and error of the input data. Our estimation is also comparable to the Clouds and Earth’s Radiant Energy System (CERES) data while shows notable improvement over the widely used National Centers for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) Reanalysis data. Using our MODIS-based datasets of total solar radiation OPEN ACCESS Remote Sens. 2014, 6 7137 and its partitioned components to drive land surface models should improve simulations of global dynamics of water, carbon and climate.


Introduction
Solar radiation (E) is the principal energy source that drives all land surface processes such as the hydrological cycle, heat transfer and greenhouse gas exchange.It is required by many land surface models and ecological models as the input to simulate land-atmosphere interactions and therefore recognized as a critical variable for the quantitative analysis of climatic change questions [1].At the time of writing, many global solar radiation datasets have been developed for various applications.For example, reanalysis products such as the National Centers for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) Reanalysis [2] (hereafter called NCEP data), the European Center for Medium Range Weather Forecasts (ECMWF) ERA-15 [3], ERA-40 [4] and ERA-Interim [5] provide long-term and multiple time-scale global solar radiation data, but are usually associated with coarse spatial resolution (e.g., >1°) and biased estimation [6,7].In addition, a number of global datasets of solar radiation based on various remotely sensed data sources and algorithms have been available at different spatio-temporal resolutions such as Earth Radiation Budget Experiment (ERBE) [8], International Satellite Cloud Climatology Project based on the NASA Goddard Institute for Space Studies radiative transfer model (ISCCP/GISS) [9], University of Maryland/MODIS (UMD/MODIS) [10], Clouds and Earth's Radiant Energy System (CERES) [11], Climate Monitoring-Satellite Application Facility (CM-SAF) [12], Land Surface Analysis-Satellite Application Facility (LSA-SAF) [13,14] and Global Land Surface Satellite (GLASS) products [15], etc.
Solar radiation propagates the atmosphere and is absorbed or scattered by atmospheric particles and molecules before reaching the land surface.Many well-parameterized atmospheric radiative transfer models have been established to simulate solar radiation transfer, with the input of atmospheric profile parameters including cloud cover and optical depth, aerosol properties, water vapor amount, and ozone column amount.Compared to General Circulation Model (GCM) estimates, satellite measurement of these parameters is considered to be more accurate and has higher spatial resolutions [16].In comparison with other approaches, the satellite-based global solar radiation datasets shall provide better estimation of downward solar radiation on land surface by using direct measurements of atmospheric profiles, and especially the incorporation of more adequate information of cloud properties [7].Benefiting from their remarkable accuracy, the satellite-based global solar radiation datasets have been frequently used in the studies on global land surface energy budget.However, these datasets have not been widely applied for land surface modeling, mainly due to the limited temporal coverage of satellite-based data, but another possible reason is that these datasets often only provide the total solar radiation (E g ) without the direct and diffuse components (E b and E d , respectively), with some exceptions such as CM-SAF, Solar radiation Data (SoDa) [17], Satel-light (http://www.satel-light.com/)and GLASS.However, growing number of studies support that these two components of solar radiation are significantly important to land surface processes and can potentially influence the greenhouse gas exchange between terrestrial ecosystems and the atmosphere [18,19].Considering the importance of direct vs. diffuse solar radiation, most popular land surface models and ecological models such as CLM 4.0 [20], JULES [21] and IBIS [22], are designed to require the direct and diffuse components of solar radiation as the inputs.Among the widely used global solar radiation datasets, to our knowledge, only the reanalysis NCEP data have both direct and diffuse solar radiation data for the input of land surface models.Hence, in many cases, running such land surface models only uses E g as input and has to partition E g into E b and E d with either empirical or semi-empirical relationships (e.g., [23,24]), which may bias the results at different locations [25].No doubt that more global solar radiation datasets with direct and diffuse components are necessary to reduce the uncertainty of land-surface modeling experiments due to uncertain input of solar radiation.
The instrument MODerate-resolution Imaging Spectroradiometer (MODIS), on-board the National Aeronautics and Space Administration (NASA) Terra and Aqua satellites, has provided large-scale, continuous, and detailed measurements of global atmospheric profile parameters [26], which makes direct estimation of the total and partitioned solar radiation possible.Up to date, there have been a number of attempts of using the MODIS products to estimate downward solar radiation [10,[27][28][29][30], but many of them were limited to cloud-free conditions [31,32], and none of them has provided the direct and diffuse products of the solar radiation.In this study, we aim to develop a new method to estimate the global-scale total, direct, and diffuse solar radiation over the MODIS era, using MODIS retrieved atmosphere and land parameters over the global terrestrial ecosystem area under all-sky conditions for the period 2003-2011.Hereafter we refer to our developed data as MODIS E g , MODIS E b , and MODIS E d for total solar radiation, direct radiation and diffuse radiation, respectively.

Overview of the Method
We use a combination of a clear-sky solar radiation model and a cloud transmittance algorithm to estimate incident downward solar radiation under all-sky conditions, utilizing the atmospheric composition data from the MODIS atmosphere products.The method simplifies the complicated radiative transfer processes by assuming atmosphere as homogeneous and plane-parallel layers without three-dimensional effects.Therefore, the one-dimensional transmittance of each layer is calculated separately and combined to estimate the surface solar radiation.The atmosphere is treated as a single layer under clear-sky conditions so that the surface solar radiation is only affected by the gaseous and particle components.Under cloudy-sky conditions, the atmosphere is divided into two layers: the layer above clouds and the cloud layer.For the purpose of simplification, the lower part under the cloud is considered to be clear.Interactions between clouds and atmospheric gaseous and particle components (e.g., aerosol) are not considered but land surface reflectance is included in the presented method.The framework of the method is summarized in the flow chart as shown in Figure 1.The following sections are organized to describe the algorithms used under clear-sky and all-sky conditions.All the key terms in the equations within these sections is documented in Table 1 [33][34][35][36].

Algorithms for Estimating the Incident Solar Radiation under Clear-Sky Conditions
We use an existing high-performance clear-sky model to estimate the incident solar radiation under the clear-sky conditions.The model, REST2, is a two-broadband model that has been well parameterized and validated [34,37].Band 1 covers the UV and visible (VIS, 0.29~0.70μm) and Band 2 covers near-infrared (NIR, 0.7~4 μm).REST2 is designed in a simple radiative-transfer manner to estimate both direct beam and diffuse radiation under cloudless condition, which are calculated as functions of individual transmittance and scattering coefficients of major atmospheric constituents (well-mixed gases, ozone, nitrogen dioxide, water vapor and aerosols).The model uses the latest solar constant value at the mean solar-earth distance as the extraterrestrial solar radiation and calculate the incident direct beam radiation under clear sky (E bi,clear ) as for each of the two band i as: The incident diffuse irradiance under clear sky (E dpi,clear ) is calculated as: For each term in Equations ( 1) and ( 2), the detailed calculation procedures follow the equations in [34], so they are not repeated here but can be found in the original documents.

Algorithms for Estimating the Incident Solar Radiation under Cloudy-Sky Conditions
Under cloudy conditions, cloud is assumed to be an additional single-plane homogeneous layer.The two-stream approximation method developed by [38] is then used to calculate cloud reflectance (ρ Ci ) and transmittance (T Ci ) for incident direct beam irradiance at each band i.
For band 1: For band 2: where    where τ C is cloud optical thickness, ω 0 is the single-scattering albedo, β 1 and β 2 are the backscattered fraction of mono-directional incident radiation at the solar zenith angle Z for band 1 and band 2. ω 0 , β 1 and β 2 are linearly interpolated using the look-up table from Table 1 of [38].A δ-function adjustment is applied to cloud optical thickness τ C and single scattering albedo (ω 0 ) to incorporate the forward peak contribution in multiple scattering [39].
where and are the adjusted cloud optical thickness and single scattering albedo, respectively.g is the cloud asymmetric factor, which uses a typical value 0.85 [35,36].The adjusted cloud optical thickness and single scattering albedo are then used in Equations ( 3)-( 8) to calculate the cloud transmittances of the direct solar radiation.The cloud transmittances of diffuse radiation (T C1,d and T C2,d ) are calculated by averaging β 1 and β 2 over the hemispheric incidence angle, as the incoming diffuse radiation is assumed to be isotropic distributed.
Therefore, in the presence of clouds, direct beam radiation on the earth surface is calculated as the cloud-fraction (f C ) weighted sum for each band: and the diffuse radiation on an ideally horizontal black surface is calculated as the weighted sum of cloud-induced diffuse radiation and diffuse radiation transmitted from clear sky for each band: , represent the diffuse radiation generated from direct beam and diffuse radiation while transmitting through clouds, respectively.m a is the relative air mass estimated by cloud top air pressure and Z.
Considering the interactions between sky and ground, the sky albedo with clouds (ρ si ) is first calculated as the weighted average of the clear sky albedo ρ si,clear (calculated in the clear-sky algorithm based on the aerosol characteristics, see [34]) and the cloud albedo ρ Ci (assuming the cloud top and bottom have the same albedo): Using the ground albedo ρ gi , therefore the backscattered diffuse radiation E ddi is: Finally the total diffuse radiation is The global irradiance is the sum of direct beam and diffuse irradiance:

MODIS data for Estimation of Global Land Surface Solar Radiation
MODIS atmosphere products continuously provide instantaneous measurements of global atmospheric profile parameters that are needed in our computing scheme.The level-3 products (MOD08 and MYD08 for daily, 8-day and monthly) have been released for global application purpose, which are gridded at 1 by 1 degree resolution spanning 1-day/8-day/1-month interval and contain all the MODIS derived atmospheric parameters.Since the spatial resolution of MODIS albedo product (MCD43C3) is 0.05 by 0.05 degree, we averaged the overlaid 20 × 20 MCD43 pixels in each of the 1 by 1 degree grid to match the resolution of the atmosphere products.
As discussed above, this method needs atmospheric inputs of ozone column, water vapor column, nitrogen dioxide column, aerosol Ångström coefficients, cloud top pressure (required by REST2), cloud fraction and cloud optical depth (required by the cloud radiative transfer scheme).The MODIS atmosphere product directly provides all the inputs except the nitrogen dioxide column and aerosol Ångström coefficients (α and β, representing the size distribution of aerosol particles and the general haziness of the atmosphere).The nitrogen dioxide has minor effect on solar radiation and is assumed to a constant (0.0002 cm).For the Ångström coefficients, MODIS products only provide the Ångström exponent α, so we use the Ångström law [40] to derive β, and then use both α and β to calculate aerosol scattering effects following REST2 algorithms.Specifically, the REST2 needs band-average optical depth to calculate the aerosol broadband transmittances, therefore it uses α and β, the same formalism as in the original Ångström law, and an estimated effective wavelength for the whole band to transform the MODIS provided spectral AOD into broadband values.Rather than using the 8-day and monthly products, we use the daily product since aerosol lifetime is normally shorter than a week [41] and the atmospheric conditions could change considerably in a week or longer.
The local overpass solar time for MODIS Terra and Aqua is around 10:30 am and 1:30 pm, respectively.Having a morning and an afternoon observation is important when integrating radiation from sunrise to sunset, because often the weather shows distinct patterns in the morning and in the afternoon [32].Therefore we use both MODIS Terra and Aqua products (MOD08_D3 and MYD08_D3, collection 051) to obtain solar radiation.The white-sky (bi-hemispherical) albedos of global land surface from the 16-day MCD43C3 products are used for calculating multiple scattering between the sky and land surface since this process mainly produces diffuse radiation.We assume that the land surface albedo remains constant during each 16-day period.The detailed information of MODIS parameters used in this study is listed in Table 2. Missing values of each parameter are filled using a gap-filling algorithm based on discrete cosine transforms [42,43] which has been validated and applied for filling the gaps of global satellite products [44].Specifically, since aerosol properties over land have a large number of missing values, we combine the most information we have together before gap filling.For example, the parameter, deep-blue Ångström exponent, is used to fill the gaps over bright land targets where the standard MODIS retrieval algorithm may fail [45].The Ångström exponent over ocean is then combined with that over land.These operations allow maximizing the spatial observation sampling size and quality and therefore are able to promote better gap-filling results.The gap-filled MODIS products is then used as the inputs for the computing scheme to estimate incident downward solar radiation for each 1 by 1 degree grid over the globe land surface area (excluding the Antarctic).

Daily Integrations of Solar Radiation
The daily-integrated solar radiation is calculated by averaging instantaneous values in hourly intervals for the 24 h in each calendar day.The MODIS Terra observation is assumed to be representative of morning atmospheric conditions, while the MODIS Aqua observation is used for the afternoon atmosphere.After setting up the computing scheme, we use MODIS data to calculate the instantaneous solar radiation, including the direct and diffuse parts with varying solar zenith angles for each day over the period of 2003-2011.The diurnal cycle of incident fluxes are determined by the extraterrestrial solar radiation flux, which is determined by the instantaneous solar zenith angle Z.

Evaluating the MODIS Estimates
For validation purpose, we collect ground measurements from the Baseline Surface Radiation Network (BSRN) sites [46].BSRN is a project of the radiation panel of the Global Energy and Water Cycle Experiment (GEWEX) and its member sites whose latitudes range from 80°N to 90°S provide the most comprehensive coverage of the major geographical regions in the world.The components of solar radiation (total, direct normal and diffuse) are continuously measured at these sites at 1 to 3 min intervals.To examine spatial characteristics of the method performance, we choose 48 sites that have available data during the period 2003-2011 and divide them into the following groups according to their geographic locations: (1) Arctic: EUR, NYA, TIK; (2) Islands: KWA, LAU, NAU, BER, CLH, LER, MAN, MNM, SAP, FUA, ISH, COC, DWN, DAR; (3) Continental sites: BIL, BON, BOS, BOU, DRA, E13, FPE, GCR, REG, PSU, SXF, LIN, CAB, CAM, PAL, PAY, CAR, CNR, ASP, DAA, ILO, SBO, TAM, TAT, TOR, XIA, PTR, BRB, SMS.The spatial distribution of these sites is shown in Figure 2.More details of the BSRN sites can be found at http://www.bsrn.awi.de/.Daily mean values of BSRN data are calculated for comparison to reduce large uncertainty caused by varying cloud conditions.Only the days that have above 90% valid observed values are used.
The NCEP data is the only global dataset that has the partitioned solar radiation among the widely used datasets.The data has coarse spatial resolution (T62 Gaussian grid with 192 × 94 points) provides daily values of surface downward solar radiation, the visible and near-infrared band direct and diffuse radiation from 1948 to present, which has been widely used for various applications [47,48].The NCEP data is a reanalysis product based on the National Meteorological Center (NMC) model [49].The NCEP algorithm of solar radiation uses a rather complicated radiative transfer scheme [50], and clouds are simulated from model derived relative humidity, vertical motion and convective rain.For purposes of comparison, we collect the following variables of NCEP data: "Downward solar radiation flux" (dswrf), "Visible beam downward solar flux" (vbdsf), "Visible diffuse downward solar flux" (vddsf), "Near IR beam downward solar flux" (nbdsf) and "Near IR diffuse downward solar flux" (nddsf).Therefore NCEP provided total, direct and diffuse radiation fluxes can be calculated as E g = dswrf; E b = vbdsf + nbdsf; and E d = vddsf + nddsf, with neglecting the UV radiation.The NCEP data is available at http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html.
In addition, the satellite-based solar radiation data from the CERES are also used for comparison.CERES instruments were launched aboard Tropical Rainfall Measuring Mission (TRMM) in 1997 and on Terra satellite in 1999.It is part of NASA's Earth Observation System (EOS) and provides long-term, accurate estimation of the radiative fluxes within the Earth atmosphere.The algorithms of CERES mainly follow the previous project ERBE and uses the instrument-derived multiple satellite-based parameters, but combine more-frequent satellite data and more-detailed cloud information from other EOS instruments such as MODIS [51].Specifically, we collect the recently released the global daily downward solar radiation at 1-degree resolution (SYN1deg) from CERES to compare with our estimates.The CERES data was downloaded from its official website (http://ceres.larc.nasa.gov/).
The daily integrated MODIS estimates and the NCEP, CERES data are interpolated at the BSRN sites using two-dimensional spline interpolation method for comparison.We notice that NCEP has coarser spatial resolution than the MODIS and CERES data, but we assume that this issue could be ignored through spatial interpolation and we don't spend discussion on it in the following analysis.Standard statistics of R 2 , Root Mean Square Error (RMSE) and mean and bias are calculated to provide quantitative assessment of the evaluation.We notice that NCEP has coarser spatial resolution than the MODIS and CERES data, but we assume that this issue could be ignored through spatial interpolation and we don't spend discussion on it in the following analysis.

Comparison at the BSRN Sites
Figure 3 shows the scatter plots between the BSRN observations and data from MODIS, NCEP and CERES.Our estimates have comparable performance with the CERES data, which are closer to BSRN observations with lower biases and RMSEs, and higher R 2 values (Table 3).The CERES data is more accurate as expected, since CERES has adapted high-frequency geostationary satellite cloud properties in addition to the MODIS cloud products, therefore has utilized more detailed atmospheric profile information, while with our method only uses twice-per-day MODIS atmosphere products.However, our estimates are consistently better than NCEP data at BSRN sites with higher R 2 values and lower RMSE values (Table 3).As shown in Figure 3, NCEP data have relatively reasonable estimates of E g, but its E b and E d have large errors.NCEP data has been evaluated by a few previous evaluation studies that indicate NCEP data systematically underestimates cloudiness [6,52,53], therefore less solar radiation is reflected, scattered or absorbed and the land surface could receive more direct beam radiation.This explains the dramatic overestimates of NCEP E b even when BSRN observes zero direct beam radiation (Figure 3b,e,h), and the overestimates of E g (Figure 3a,d,g) but underestimates of E d (Figure 3c,f,i).In addition, since UV radiation is not included in the NCEP data, which is about 3%-5% of total solar radiation (calculated from data in American Society for Testing and Materials G173-08 Reference Spectra http://rredc.nrel.gov/solar/spectra/am1.5/), the actual bias of NCEP E g after considering UV radiation should be even higher than 38.11W•m −2 as shown in Table 3.
Overall our estimates of solar radiation using MODIS products show good agreements with the observed data.At all the 48 BSRN sites and under the all-sky conditions, the statistics of MODIS data and BSRN observations are: R 2 = 0.88 and RMSE = 34.59W•m −2 for E g ; R 2 = 0.79 and RMSE = 41.98 W•m −2 for E b ; and R 2 = 0.64 and RMSE = 28.06W•m −2 for E d .These statistics are close to Wang and Pinker's results when they compare the UMD/MODIS products and 31 BSRN sites [10].Despite the overall good performance, errors of E b and E d estimations are relatively higher and the MODIS data generally underestimates E b and overestimates E d over all the geographical regions.These errors may be caused by multiple reasons.The measuring errors due to systematic problems such as the measurement "cosine response error" of pyranometers, which are widely used at BSRN sites [46], exists at all the BSRN sites and varies under different solar zenith angles.Beside the error induced by interpolation, the simplified radiative transfer processes and the assumptions of the calculating method, the spatial-temporal resolution and the uncertainty of the MODIS atmospheric products can all contribute to the errors of the estimated solar radiation, which is discussed in the below sections.
Under the clear-sky conditions, the results at all the BSRN sites show similar statistics as the cloudy-sky results but generally lower RMSEs.Note that the clear-sky algorithms are exactly the same as that are described in Gueymard's work in [34], but the accuracy of our estimates are notably lower than that reported in that document.The inconsistent accuracy of the algorithm can be caused by quality of the input data, including MODIS derived atmospheric composition parameters and land surface albedo, as well as the neglected topography at the BSRN sites, which can be rather heterogeneous in a 1-degree grid.The bias of NCEP clear-sky values is largely reduced comparing to cloudy-sky conditions, which further supports the suggestion that the cloudiness characterization is a major error source of the NCEP data.Consistent with the other two datasets, the MODIS data generally has the lowest accuracy at the Island sites (Table 2).All the three datasets have the lowest R 2 and highest RMSE at the Island sites.These results are consistent with previous findings, most likely due to the island effect and mixed grid of land and ocean along the coasts [10], which could substantially affect the accuracy of the land surface albedo and therefore the multiple reflection process.Comparing to the other continental sites, sites located in the arctic zone are often with bright land surface and high solar zenith angle, which make the MODIS retrieval of clouds and other atmospheric compositions more difficult.In addition, our method is less reliable when the solar zenith angle is high, since it is a 1D model that cannot deal with the geometrical effects of scattering.Both reasons may contribute to larger bias of estimated solar radiation.

Global-Scale Comparison
Over the study period, our method estimates the global land surface mean total solar radiation to be 182.5 W•m −2 , which is 33 W•m −2 lower than the NCEP data, but close to other satellite-based estimates such as the International Satellite Cloud Climatology Project (ISCCP) data [9] and the CERES data [51,54] which are 188.8 and 180.1 W•m −2 , respectively.For E b and E d , comparison shows that our estimated E b and E d are 39.2 W•m −2 lower and 6.6 W•m −2 higher than NCEP data respectively, which is consistent with site-level conclusions.
Figure 4 shows the spatial patterns of our estimated land surface solar radiations, which are comparable to the NCEP data and very similar to the CERES data.Over the study period, all the data show that most of the high surface downward solar radiation takes place at between 30°S to 30°N, especially in the arid regions where are lack of clouds and precipitation.High-latitude regions have the lowest solar radiation because of the averagely low solar zenith angle.Southern China experiences relatively low solar radiation comparing with the other regions at the same latitude, probably because of the substantial influence of the East Asian Monsoon as well as heavy aerosol loadings caused by industrial air pollutions.Figure 4h shows that the MODIS solar radiation has the similar latitudinal pattern to the CERES and NCEP data.The MODIS data match the CERES data very well in most latitudes, but slightly lower than CERES data in northern high latitudes, probably due to underestimation of direct radiation as suggested in Figure 3b.The MODIS and NCEP data have major discrepancies in both of the magnitudes and spatial patterns of the direct and diffuse solar radiation (Figure 4b,c,e,f,i,j).The largest discrepancies happen in the northern middle and low latitudes, where NCEP provides higher direct radiation but lower diffuse radiation.This region, including the majority of southeastern Asia and part of central Africa as well as Amazon tropical forests, plays a significant role in global climate change.However, the average biases between 0-30°N of the NCEP and MODIS (NCEP-MODIS) direct and diffuse solar radiation are 40.9 and −15.6 W•m −2 , respectively, equivalent to 104% and 236% of the global mean biases.
Comparisons show the MODIS data well captured the seasonal change as well as the inter-annual change of solar radiation (Figure 5).In each image of Figure 5, each row of pixels shows the global zonal mean values across the study period (from 2003 to 2011).The highest and lowest solar radiation estimates alternately occur in high-latitude regions because of polar days and nights.Zonal means of solar radiations are low in equatorial regions due to the small area of land.Consistent with the above findings, MODIS E g matches the CERES E g very well at both spatial and temporal domain while exhibiting magnitude difference comparing with NCEP data.MODIS-NCEP difference of E g is mostly located in middle-latitude and low-latitude regions in the middle of a year, mainly attributed to the differences of E b in these regions and the compensation from the differences of E d .These findings are consistent with previous comparisons of NCEP data with ISCCP/GISS data [6].

Discussion
Solar radiation is a key driver of the land surface models, which usually run at large scales.The objective of this study is to develop a computationally efficient, and easy-to-use method for estimating the partitioned downward solar radiation at the global scale so as to be used as the input the land surface models.Therefore, the calculation scheme of this method has a rather simpler structure and simplified many radiative transferring processes comparing to the elaborate models such as MODTRAN and HITRAN.For instance, the method doesn't distinguish the aerosol and cloud types, which however have considerable impacts on distribution of the asymmetrically scattered radiation.The plane-parallel, homogeneous "two-layer" assumption of the cloudy sky may not be very realistic under most conditions, particularly when the cloud over is the discontinuous cumulus-type.Interactions between ground and atmosphere has been considered, but interactions between the clear-sky layer and the cloud layer is not included in the method, although the multiple scattering could be significant with heavy aerosol load.For instance, the multiple scattering between the clear-sky and the cloud layers may add more downward path radiation and result in larger estimation of the solar radiation.However, the overall accuracy of the MODIS data indicates that these simplifications are acceptable and estimates by this method show significant improvement over the widely used NCEP dataset, and therefore are more accurate and suitable for driving the large-scale land surface models to better predict land surface processes [55].
The uncertainty and error of the MODIS derived atmospheric parameters essentially influences the accuracy of the estimated solar radiation.The retrieval algorithms of MODIS atmosphere products allow relatively high accuracy of estimating aerosol optical depth but less accurate aerosol size parameters.Retrieval of aerosol properties, which is directly related to the aerosol attenuating effects, is also highly influenced by the large solar zenith angle, cloud and bright surface.The relatively large errors as shown in the evaluation results at the arctic sites, for instance, may be attributed to this reason.In addition, using the MODIS products from Terra and Aqua, which provide only two observations per day, is not adequate to represent high-frequency diurnal variations of the atmospheric conditions.This shortage should be more significant in the tropical regions or during the summers, since the atmospheric dynamics are more active and causing rapid change of the atmospheric profile.
The assumption of constant atmosphere conditions in the morning and afternoon is therefore over-simplified for these situations and may be problematic in the daily integration and could be a large error source.Another important source of error is the heterogeneity of the atmosphere compositions and land surface within the grid.A 1° × 1° grid roughly covers 10 4 km 2 area, in which the atmospheric compositions usually are not homogeneously distributed as assumed.Particularly, if the land surface within a grid has rough terrain, the surface-received solar radiation can be dramatically affected by the surface topography [56][57][58].For example, the bottom of a valley in a mountainous area may not be able to receive direct-beam solar radiation with large solar zenith angle and only part of the diffuse solar radiation because of the shading effect of the nearby mountain peaks.However, the spatial-resolution issue of MODIS atmospheric products may be not problematic for the application of the method at large scale, but could contribute to explain the inconsistency between our estimates and BSRN observations.Another important aspect of the error is the spatial gaps in the MODIS atmosphere products, especially the aerosol property data.The sun-synchronous and near-polar orbit of the Terra and Aqua satellite results in relatively sparse overpasses in the low-latitude area and therefore causes spatial gaps of atmospheric composition retrievals; the complex retrieval of aerosol always fails in the bright-surface regions, resulting in large gaps of aerosol-related parameters over the deserts and the high-latitude area.In this case, the gap-filling method is applied to estimate missing data and introduces errors into the calculation of solar radiation.There are many existing gap-filling methods (e.g., Kridging, linear/nonlinear spatial interpolation, etc.), most of which are based on the principle that uses the information of nearby pixels to approximate the gaps with the different definitions of the spatial correlations.The gap-filling method used in our computing scheme is one of them, but incorporates the temporal dimension rather than only the spatial domain, thus provides robust gap filling and has been applied successfully in previous studies [44].However, the accuracy of any gap-filling method depends on the useful nearby data of the gap.Once the gap is too large, e.g., the gaps of aerosol properties in the arctic regions, the gap-filling method may not be able to work accurately and therefore results in significant errors.For example, the REST2 model is sensitive to the Ångström exponent α.According to the results in Gueymard 2008 [34], in a case with low solar zenith angle (i.e., Z = 0), medium aerosol haziness (β = 0.5) and a fixed reference atmosphere conditions (clear sky, p = 1013.25 mb, water vapor column = 1.5 cm, ozone column = 0.35 cm, NO 2 column = 0.0002 cm), a unit error of α caused by the gap-filling process could roughly result in 110 W•m −2 error of the model estimated E g .This error is smaller under higher solar zenith angles but could be even larger when experiencing heavier aerosol loads (i.e., larger β).This could be used to explain the large discrepancies between the zonal mean solar radiation between the CERES and MODIS data in Figure 4h.
As indicated above, the overall accuracy of solar radiation estimated by our method is comparable with another well-known satellite-based data, the CERES data.Compared to many satellite-based algorithms, our method is simple and easy to use with the publicly available MODIS products to rapidly retrieve the most recent global solar radiation.Furthermore, the advantage of the method is that one can partition direct and diffuse components in both the visible and near-infrared broadbands, which are necessary for most land surface models, and therefore promotes better estimation of land surface energy balance, hydrological cycling, as well as regional and global climate [59,60].This method can also contribute to explore the relative importance of atmospheric compositions in affecting the land surface processes.The direct radiative effect of each atmospheric composition, e.g., aerosol and clouds, can be easily examined by turning on/off the correlated terms in the computing scheme [61].

Conclusions
This study developed a simple but efficient method to estimate two-band downward solar radiation, including the direct and diffuse components over the global land surface.Ground-based measured data from 48 BSRN sites distributed over the globe was used to evaluate our products.Evaluations show that our MODIS-based estimates are in good agreement with the observed values (R 2 = 0.88, bias = 5.72 W•m −2 and RMSE = 34.59W•m −2 for the total solar radiation; R 2 = 0.79, bias = −4.12W•m −2 and RMSE = 41.98 W•m −2 for the direct beam radiation; and R 2 = 0.64, bias = 9.84 W•m −2 and RMSE = 28.06W•m −2 for the diffuse radiation), and have lower but comparable accuracy with the satellite-based CERES data.Our estimated global land surface mean total solar radiation over the study period is 182.5 W•m −2 , which is close to other satellite-based ISCCP and CERES data.Comparisons indicate that this study provides substantially better estimates of total, direct and diffuse solar radiations than the widely used NCEP reanalysis data, therefore can provide a reliable alternative solar radiation dataset for the scientific use.
The method introduced in this study intends to keep the simplicity of the modeling framework, which only requires the inputs of the time, and the publicly available MODIS products.While the simplified processes may bring considerable errors into the calculation, it however enables users to efficiently estimate the solar radiation for various purposes in the land surface modeling studies.

Figure 1 .
Figure 1.Schematic diagram of the method used in this study to estimate each component of the downward solar radiation.The definition of each variable in the flowchart can be found in Table1.

Figure 2 .
Figure 2. Map of the selected Baseline Surface Radiation Network (BSRN) sites used in this study.

Figure 3 .
Figure3.Scatterplots of the MODIS, NCEP and CERES daily total, direct and diffuse solar radiation with BSRN observations under all-sky conditions for different geographic regions (a-i).x-axis is BSRN observations and y-axis is the data from MODIS, NCEP and CERES.The red, green and blue dots represent the NCEP, MODIS and CERES data, respectively.

Figure 4 .
Figure 4. Comparison of the spatial pattern of the MODIS (a-c), NCEP (d-f) and CERES (g) temporal mean total, direct and diffuse solar radiation over the period 2003-2011.(h), (i) and (j) show the zonal mean of the solar radiation over the global land surface.

Figure 5 .
Figure 5.Comparison of the temporal variation of the MODIS (a-c), NCEP (d-f) and CERES (g) zonal mean solar radiation received by global land surface over the period 2003-2011.The x-axis is the time series of days from 2003 to 2011; the y-axis is the latitude.

Table 1 .
1. Key variables in the computing scheme.
[34]lar zenith angle Local time; longitude and latitude; algorithms from [33] P Cloud top pressure MODIS cloud top pressure data E 0,i Extraterrestrial solar constant for each of the two broad bandsValues by[34]

Table 2 .
MODIS products used in the computing scheme.

Table 3 .
Statistics of the MODIS, NCEP and CERES daily total, direct and diffuse solar radiation with BSRN observations for different geographic regions.CERES data only provides E g hence no statistics for E b and E d .n is the number of valid observed.