Parameterization of clear-sky surface irradiance and its implications for estimation of aerosol direct radiative effect and aerosol optical depth

Aerosols impact clear-sky surface irradiance () through the effects of scattering and absorption. Linear or nonlinear relationships between aerosol optical depth (τa) and have been established to describe the aerosol direct radiative effect on (ADRE). However, considerable uncertainties remain associated with ADRE due to the incorrect estimation of (τa in the absence of aerosols). Based on data from the Aerosol Robotic Network, the effects of τa, water vapor content (w) and the cosine of the solar zenith angle (μ) on are thoroughly considered, leading to an effective parameterization of as a nonlinear function of these three quantities. The parameterization is proven able to estimate with a mean bias error of 0.32 W m−2, which is one order of magnitude smaller than that derived using earlier linear or nonlinear functions. Applications of this new parameterization to estimate τa from , or vice versa, show that the root-mean-square errors were 0.08 and 10.0 Wm−2, respectively. Therefore, this study establishes a straightforward method to derive from τa or estimate τa from measurements if water vapor measurements are available.

( ↓ E g in the absence of aerosols). ADRE is widely reported in the literature, based on a combination of measurements of ↓ E g and aerosol optical depth (τ a ) [4][5][6][7] . ↓ E g 0 cannot be obtained straightforwardly from observations since the atmosphere almost has aerosols present. One of the difficulties relating to the derivation of ADRE stems from the need to accurately estimate ↓ E g 0 7,8 . Radiative transfer model or single-layer clear-sky solar radiation model can be used to calculate ↓ E g 0 9-12 and thereby ADRE is obtained, however, this method is sensitive to the calibration uncertainties of pyranometer (independent of model calculation) and dependent on model assumptions about the atmospheric parameters 13 . ↓ E g 0 can be estimated from observations, this method avoids dependence on a model, furthermore, ADRE estimation should not be very sensitive to the calibration errors since ↓ E g 0 is derived from observations 13 . A linear relationship of ↓ E g to τ a has been popularly assumed and ↓ E g 0 is then derived by linear regression analysis of ↓ E g and τ a 4, 5,13 . Some studies suggested that an exponential decay of ↓ E g with an increase in τ a would be expected according to Beer-Lambert Law, which is especially true for cases with τ a values larger than 0.5 7,8,13 . However, contrary to expectation, a better estimation of ↓ E g 0 was derived from the linear regression than using the exponential relationship. This was thought to be because the systematic underestimation of ↓ E g 0 by the linear regression was compensated by the positive correlation between τ a and water vapor content (w) 8 .
In this study, we show that these previous methods produce a systematic bias in the derivation of ↓ E g 0 and thereby result in an overestimation of ADRE. A new parameterization of ↓ E g is developed based on global Aerosol Robotic Network (AERONET) data, in which the relationships of ↓ E g to the cosine of the solar zenith angle (μ), τ a and w have been established by using a combination of nonlinear equations. The results show that the mean bias error (MBE) of the estimations of ↓ E g 0 decreases from 4-7 W m −2 to 0.32 W m −2 . The same improvement in the estimation of ADRE would be expected when using the new method. Furthermore, one of the important advantages of this parameterization is that a straightforward method to derive τ a from ↓ E g , or vice versa, has been established. Specifically, we find it is possible to derive τ a from ↓ E g with a root-mean-square error (RMSE) of 0.08, and vice versa with an RMSE of 10.0 W m −2 .

Results
The solid dots in Fig. 1 represent the scatter between τ a at 550 nm (τ a hereafter) and ↓ E g at two narrow μ and w ranges. The analysis is firstly performed on data points with a very narrow range of μ (~0.2°) and w (10%w), to isolate the effect of τ a on ↓ E g . The mean ↓ E g 0 amounts, using the AERONET calculation and fitting ↓ E g 0 values using regression analysis on the basis of Eqs (1)-(4) in the Method section, are also presented. The performance of these equations is evaluated by the agreement in instantaneous ↓ E g 0 between the AERONET calculations and the regression analysis results. We can see that the best performance is achieved using Eq. (4), the new parameterization proposed in this study, which produces a difference in ↓ E g 0 between the mean AERONET calculation and the regression analysis result (∆ ↓ E g 0 ) of nearly zero. This nearly zero ∆ ↓ E g 0 is in fact always derived using Eq. (4) for the full μ and w ranges (not shown). In contrast, ∆ ↓ E g 0 , when using Eqs (1), (2) and (3), varies from a few to tens of W m −2 in these two cases. The fact that Eq. (3) occasionally produces unrealistic results that indicate the poor performance of the nonlinear regression analysis 8 , we eliminate it hereafter. Poorer performance of Eqs (1) and (2) than Eq. (4) is further shown by the histogram of ∆ ↓ E g 0 for Eqs (1) and (2)   extrapolate ↓ E g to zero τ a rather than linear regression. This conclusion is also supported by the fact that smaller residuals of the regression analysis are derived from Eq. (2) than from Eq. (1). This is expected because attenuation of ↓ E g by aerosols shows nonlinear decay, as implied by Beer-Lambert Law. The difference between Eq. (2) and Eq. (4) is the introduction of a new parameter, C, into Eq. (4). The value of C is nearly always lower than 1.0 (the exact value of Eq. (3)), which is one of the most important reasons for the better performance of Eq. (4). Eq. (3) is somewhat similar to Beer-Lambert Law, which depicts the attenuation of solar direct radiation by aerosols; however, some part of the attenuation of solar direct radiation is backscattered to the surface, which is certain to enhance ↓ E g . It is therefore expected that C should be lower than 1.0. In addition, the effect of aerosols on ↓ E g should not be independent from μ, since μ governs the transfer path of photons. This implies that C should vary with μ, which is reflected in the following analysis. Figure 3 presents the dependence of ↓ E g 0 (A in Eq. (4)), on w and μ. Variation of ↓ E g 0 is mainly governed by μ, which is somewhat modulated by w at the same value of μ. Therefore, parameter A for a given amount of w is firstly simulated using a power law function of μ (Eq. (5) in the Methods section). The first and most important reason for the selection of a power law function is because it models the simple physics of the situation with only two parameters 14,15 . Parameter a 1 of Eq. (5) represents expected measurements of ↓ E g 0 for a μ of 1. Parameter a 2 governs the ↓ E g 0 variation with μ. The second reason for selection of a power law function is that it provides a faithful approximation to the data.
Since the dependence of ↓ E g 0 on μ is depicted by parameters a 1 and a 2 of Eq. (5), the w effect on ↓ E g 0 is further simulated through the parameterization of a 1 and a 2 as a function of w. Figure 4 shows the relationships of a 1 and a 2 to w. We can see that the effect of water vapor per one unit of w on ↓ E g 0 decreases as w increases, which is expected since the relative increase in water vapor absorption gradually decreases as w increases 16 . These relationships are simulated using an exponential equation. The RMSEs of the regression analysis for a 1 and a 2 of Eq. (5) are 1.37 and 0.0004, respectively, indicating a faithful approximation. A parameterization of ↓ E g 0 to μ and w can then be established through a combination of Eqs (5)-(7) that leads to Eq. (8). Therefore, it is straightforward to calculate ↓ E g 0 from Eq. (8) if w is available, since μ can be calculated from location and time very accurately.
Further analysis of parameters B and C of Eq. (4) shows that both parameters are moderately related to μ and w, which thereby leads to a parameterization of ↓ E g as a function of τ a , μ and w. As shown in Fig. 5, parameters B and C are approximated well using Eqs (9) and (10). In terms of the variability of both parameters, 99.8% is explained by the regression analysis. Since the relationship between instantaneous ↓ E g to τ a as well as w for a specified value of μ is established, therefore, a straightforward method is developed that can be used to derive ↓ E g if τ a and w are available from another source, such as satellite remote sensing. On the other hand, it can also be used to derive τ a if ↓ E g and w are available from sources such as the Baseline Surface Radiation Network (BSRN). The proposed method is evaluated by using the 20% of validating data and BSRN data at Xianghe. Figure 6a shows the comparison of instantaneous AERONET ↓ E g 0 values of validating data and calculations from Eq. (8) based on validating AERONET τ a and w. The MBE is 0.33 W m −2 , one order magnitude smaller than the results from Eq. (1) and (2), even though the latter is derived from the training data. Since ↓ E g 0 is simulated well by Eq. (5), the ADRE derivation based on Eq. (5) should be very close to that derived from the AERONET model calculations. This expectation is supported by Fig. 6b, and parameters (a) a 1 and (b) a 2 ofEq. 5. The curve represents the regression result using Eqs (6) and (7). The figure was produced using MATLAB. values and compared with AERONET τ a products. Figure 7a shows that ↓ E g can be estimated with an MBE and RMSE of 0.02 and 10.0 W m −2 , respectively. This certainly relies on the fact that both τ a and w are available. On the contrary, if w and ↓ E g are available and τ a is not known, τ a can be retrieved from ↓ E g and w using this parameterization. The MBE and RMSE values of τ a retrievals are 0.0005 and 0.08, respectively (Fig. 7b). Measurements of ↓ E g and τ a at Xianghe 7 , a BSRN and AERONET station in China are used to further evaluate the effectiveness of the parameterization of ↓ E g . The results are shown in Fig. 8. The estimations of ↓ E g from AERONET τ a and w products using the proposed parameterization method agree with the BSRN measurements very well, with an MBE and RMSE of − 3.9 and 12.5 W m −2 , respectively. On the other hand, the retrievals of τ a from the measurements of ↓ E g and w are compared with AERONET τ a products and the MBE and RMSE values are − 0.03 and 0.08, respectively. These results once again proved the reliability of the proposed parameterization method.

Uncertainty analysis.
In the parameterization of ↓ E g 0 (Eq. 8 of the Method section), surface albedo effect was excluded that likely produced bias in the estimation of ↓ E g 0 . Figure 9 shows the scatter-plot of  surface albedo and ∆ ↓ E g 0 , from which we can see a significant negative correlation between both quantities. Uncertainty of 0.1 in surface albedo may produce 1-3 W m −2 bias in ↓ E g 0 . τ a is the dominant aerosol optical property driving the variation of ↓ E g and therefore ADRE. However, aerosol absorption also plays an important role in ADRE 17 , which shows a wide range of variations and thereby may lead to uncertainties in the parameterization of ↓ E g , since it was not accounted for. Remarkable impact of aerosol single scattering albedo at 550 nm (ω 550nm ) on the estimation of ↓ E g is presented in Fig. 10. ↓ E g changes as a result of uncertainty of ω 550nm (0.03) was estimated to be a few Wm −2 that depends on optical path and τ a . Significant impacts of ω 550nm on estimation of ↓ E g from τ a are further evidenced in Fig. 11 in which ∆ ↓ E g 0 shows a significant correlation to ω 550nm . The best estimation is achieved for ω 550nm of ~0.90 that is close to the median value of ω 550nm of AERONET data points.
In the above error analysis of ↓ E g 0 , water vapor is assumed to be known without any uncertainty. This is, of course, not realistic, since water vapor products from AERONET, satellite measurements are not free of uncertainty. By differentiating Eq. (8) with respect to w, the uncertainty of ↓ E g 0 is estimated to be < 3 W m −2 that is slightly dependent on optical path if the uncertainty of w is assumed to be < 10% of w. The uncertainty of ↓ E g estimation was estimated to < 2 W m −2 if AERONET τ a products with uncer- tainty of 0.01 ~ 0.02 are used. However, this may reach 10 W m −2 if satellite τ a products are used since their uncertainty was estimate to be 20% of τ a over land 18 . The uncertainly of BSRN ↓ E g measurements is estimate to be 2% 19 , which may lead to the uncertainty of τ a < 0.02.

Discussion
↓ E g is one of the key parameters governing a large number of diverse surface processes, and therefore accurate measurement or estimation of ↓ E g is significant 1-3 . Surface ↓ E g networks are still limited in  Figure 10. Scatter-plot of ω 550nm to ↓ E g for three τ a values and two solar zenith angles. The figure was produced using MATLAB.
Scientific RepoRts | 5:14376 | DOi: 10.1038/srep14376 spatial coverage; therefore, satellite remote sensing is a promising method for the accurate estimation of clear sky ↓ E g . Some highly complex algorithms have been developed to estimate ↓ E g from satellite remote sensing data 20 . Given that τ a values have been inverted from a few spaceborne radiometers since 2000 with good quality (e.g. the Moderate Resolution Imaging Spectraradiometer 18 ), establishment of the parameterizations in this study provide a straightforward method to calculate clear sky ↓ E g from such satellite aerosol products.
Broadband pyranometer measurements have been used to retrieve τ a , which is e expected to be a promising method to build a long-term dataset of aerosol loading, since early pyranometer measurements can be tracked to the beginning of the last century 1 . Broadband direct solar radiation is widely used in these previous studies [21][22][23] . The method proposed in this study is based on global solar radiation measurements that are available more often than direct solar radiation. For example, there are only dozens of stations with direct solar radiation measurements; however, global solar radiation is measured at more than 100 stations in China. Certainly, it should be noted that measurement of global solar radiation is impacted by contamination of the upward facing glass dome, leveling of the instrument and cosine response of the pyranometer. The disadvantage of using direct solar radiation measurements is that it is occasionally disturbed by solar tracker malfunctions. Furthermore, both methods are impacted by calibration uncertainties.
The reason for only considering τ a in the proposed method is that the availability of aerosol absorption is very limited, especially from satellite remote sensing. Similar analysis can be performed for a specified area characterized by a special aerosol type, e.g. dust aerosol in desert regions or biomass-burning aerosol in tropical forest regions. In this case, better performance of the parameterization is expected since aerosol absorption shows much less variation for the same aerosol type 24 . Furthermore, lower variation of surface albedo, ozone amount and surface elevation is also expected to reduce the random error of the parameterization.

Conclusion
Solar zenith angle, aerosol and water vapor are the three most important physical quantities governing the variability of ↓ E g . Based on a large quantity of AERONET τ a , w, and ↓ E g products, the effects of these quantities on ↓ E g are fully considered, leading to an effective parameterization of ↓ E g as a nonlinear function of these three quantities. The first advantage is that an accurate estimation of ↓ E g 0 is achieved, which ultimately results in a significant improvement of ADRE estimation compared to previous  Figure 11. Scatter-plot of ω 550nm to the difference in ↓ E g 0 between AERONET product and estimation using the proposed method for six solar zenith angle ranges. The figure was produced using MATLAB.
Scientific RepoRts | 5:14376 | DOi: 10.1038/srep14376 methods. The second is that a straightforward method has been established to estimate ↓ E g from τ a , or vice versa, if w is available. It is expected that potential applications of this new parameterization in the estimation of ↓ E g and τ a will arise in the near future.

Methods
I used ↓ E g , ↓ E g 0 , τ a and w products from those Aerosol Robotic Network (AERONET) sites with an elevation of less than 0.8 km (to eliminate the Rayleigh scattering effect on the analysis) (http://aeronet. gsfc.nasa.gov). The AERONET is a federation of ground-based remote sensing aerosol networks that is composed of more than 700 stations across the world (see Supplementary Fig. S1) 25 . The AERONET products were used because they cover different aerosol types (0 < τ a < 3.0; 0.65 < ω 550nm < 1.0; − 0.2 < α 440_870nm < 2.5, see Supplementary Fig. S2) and thereby realistically represent the aerosol direct effect on ↓ E g . Furthermore, the availability of ↓ E g 0 data provides a benchmark for the evaluation of the parameterizations. Uncertainty of τ a was estimated to be 0.01-0.02 26 . ↓ E g and ↓ E g 0 were calculated using the discrete ordinates radiative transfer model with and without aerosols 27 . The ↓ E g values agree with pyranometer measurements, with the relative difference varying from 0.98 to 1.02 28,29 . Of the 950,000 AERONET data points with surface albedo at 440 nm less than 0.25 (to reduce surface albedo effect on the analysis), I randomly select 80% of data points to develop the parameterization and the remaining 20% were used as test data. The analysis flow chart was presented in Supplementary  (Fig. S3) that was described as follows.
To isolate the dependence of ↓ E g on τ a , the AERONET data were firstly divided into subgroups according to θ s and w. The range of θ s was 0.2°. The range of w was 10% of w (the measurement uncertainty 25 ), respectively. The amount of ↓ E g was normalized for the average Earth-Sun distance and cosine correction of ↓ E g was performed within ranges to its midpoints. Three equations used in the literature were considered to represent the dependence of ↓ E g on τ a : These equations were compared with the following new equation proposed in this paper: The performance of these methods was evaluated using the ↓ E g 0 difference between the mean AERONET model calculations and the regression analysis result (∆ ↓ E g 0 ). Given that Eq. (3) occasionally produces unrealistic results, we eliminated it in the comparison.
To derive ↓ E g 0 for varying μ and w, ↓ E g 0 was further parameterized as follows: where a 1 and a 2 was found to relate to w as follows: Therefore, the parameterization of ↓ E g 0 was finally developed through a combination of Eqs (5), (6) and (7). It was found that parameters B and C of Eq. (4) show moderate variation with μ and w, which was then simulated by the following equations: The parameterization of ↓ E g to μ, w and τ a was finally established. This parameterization can be used to estimate ↓ E g by using a combination of Eq. 4 and 8-10 if τ a and w are available, and conversely, τ a can be directly calculated from ↓ E g and w. μ can be accurately calculated from location and time.