Retrieval of gridded aerosol direct radiative forcing based on multiplatform datasets

Atmospheric aerosols play a crucial role in regional radiative budgets. Previous studies on clear-sky aerosol direct radiative forcing (ADRF) have mainly been limited to site-scale observations or model simulations for short-term cases, and long-term distributions of ADRF in China have not been portrayed yet. In this study, an accurate fine-resolution ADRF estimate at the surface was proposed. Multiplatform datasets, including satellite (MODIS aboard Terra and Aqua) and reanalysis datasets, served as inputs to the Santa Barbara Discrete Atmospheric Radiative Transfer (SBDART) model for ADRF simulation with consideration of the aerosol vertical profile over eastern China during 2000–2016. Specifically, single-scattering albedo (SSA) from the Modern-Era Retrospective Analysis for Research and Application, Version 2 (MERRA-2) was validated with sun photometers over eastern China. The gridded asymmetry parameter (ASY) was then simulated by matching the calculated top-of-atmosphere (TOA) radiative fluxes from the radiative transfer model with satellite observations (Clouds and the Earth's Radiant Energy System, CERES). The high correlation and small discrepancy (6–8 W m⁻²) between simulated and observed radiative fluxes at three sites (Baoshan, Fuzhou, and Yong'an) indicated that ADRF retrieval is feasible and has high accuracy over eastern China. Then this method was applied in each grid of eastern China, and the overall picture of ADRF distributions over eastern China during 2000–2016 was displayed. ADRF ranges from −220 to −20 W m⁻², and annual mean ADRF is −100.21 W m⁻², implying that aerosols have a strong cooling effect at the surface in eastern China. With the economic development and rapid urbanization, the spatiotemporal changes of ADRF during the past 17 years are mainly attributed to the changes of anthropogenic emissions in eastern China. Our method provides the long-term ADRF distribution over eastern China for the first time, highlighting the importance of aerosol radiative impact under climate change.

1 Introduction

Atmospheric aerosols play a significant role in air quality, regional–global climate, and human health (Wang et al., 2018, 2019). Aerosols can directly absorb and scatter solar radiation, and they indirectly affect cloud formation and precipitation by acting as cloud condensation nuclei or ice nuclei (Twomey, 1977; Rosenfeld, 1999). Large amounts of scattering aerosols can generally attenuate incoming solar radiation. This reduction in surface radiation significantly impacts the surface temperature, crop growth, and solar energy availability (Chameides, 1999; Liao et al., 2015). On the other hand, highly absorbing aerosols, such as black carbon, can warm the atmosphere, alter regional atmospheric stability, and even influence the large-scale circulation and hydrologic cycle with significant regional climate effects (Menon et al., 2002; J. Wang et al., 2009). Aerosol direct radiative forcing (ADRF) is a good metric for evaluating the impact of aerosols on radiation by absorption and scattering and is defined as the difference between the net radiative flux of earth–atmosphere systems with and without aerosols. Anthropogenic aerosols produce a global mean negative direct radiative forcing of $-\mathrm{0.35}\pm \mathrm{0.5}$ W m⁻² of ADRF, which has dampened the warming effect of greenhouse gases (IPCC, 2013). However, the current assessment of ADRF remains highly uncertain. This uncertainty mainly results from the large variations in aerosol concentrations, chemical compositions, optical properties, mixing states, and vertical profiles (Haywood and Boucher, 2000; Tian et al., 2018a). Therefore, an accurate and feasible method for ADRF retrieval is greatly needed.

Reduction in these uncertainties requires the integration of different techniques and datasets (e.g., surface measurement, model simulation, and satellite remote sensing) (Yu et al., 2006). To better understand aerosol optical properties and their radiative effect, several ground-based networks have been established worldwide, such as the AEROsol Robotic Network (AERONET) (Holben et al., 2001), Global Atmosphere Watch precision-filter radiometer (GAW-PFR) network (Nyeki et al., 2015), China Aerosol Remote Sensing Network (CARSNET) (Che et al., 2009), and Chinese Sun Hazemeter Network (CSHNET) (Xin et al., 2007). Moreover, intensive field experiments have been carried out over China, in Beijing, Xianghe, Taihu, Wuhan, Shanghai, and Lanzhou (Li et al., 2003; He et al., 2012a; Wang et al., 2014; Yu et al., 2016a; Gong et al., 2017; Zhang et al., 2018). Such measurements are conducive to gaining wider knowledge of aerosol properties, which is helpful for improving the performance of satellite and model simulations through synthesis. Nevertheless, available measurements are usually restricted in terms of spatial and temporal coverage. In addition to surface measurements, model simulations play an indispensable role in the estimation of the aerosol radiative effect at the global scale and excel in predicting past or future trends of ADRF (Chang and Liao, 2009; Qiu et al., 2017). Meanwhile, model simulations are subject to large uncertainties in terms of emissions, transport, and physical and chemical parametrization schemes (Ruiz-Arias et al., 2013).

Figure 1 The map of research area, topography, major lakes, and mountains in eastern China. The red circles denote the locations of three pyranometers (Baoshan, Fuzhou, and Yong'an). This figure was generated by ArcGIS, version 10.2. Map source: Map World (National Platform for Common Geospatial Information Services, http://www.tianditu.gov.cn/, last access: 1 November 2019).

Compared to the above methods, satellite remote sensing has an outstanding advantage of delivering aerosol information with higher spatial resolution and larger spatial coverage. Using solely satellite data or a combination of model simulations and observations, many methods have been developed to retrieve global and regional ADRF estimates (e.g., Yu et al., 2004; Bellouin et al., 2005; De Graaf et al., 2013). However, these studies have mainly concentrated on the top-of-atmosphere (TOA) radiation budget. Thus far, long-term estimates of the surface ADRF distribution have rarely been addressed and few studies gave a full picture of surface ADRF over land (e.g., Thomas et al., 2013; Chung et al., 2016). This lack of research is because satellites are unable to measure surface-level radiative fluxes directly. Furthermore, the retrieval of aerosol microphysical parameters remains challenging, including single-scattering albedo (SSA; see Table S1 in the Supplement for the abbreviations) and asymmetry parameter (ASY). Many attempts have been made to solve this key problem. For instance, Thomas et al. (2013) adopted prescribed aerosol properties from the literature to estimate surface ADRF. Fu et al. (2017) took aerosol optical parameters from some AERONET sites as representative of the entire region to conduct grid-cell ADRF simulations. Undoubtedly, additional uncertainty was introduced by the assumption of aerosol optical representativeness in the temporal and spatial dimensions. Some studies also nudged global model simulations towards AERONET SSA to obtain the aerosol parameters (Chung et al., 2016). With the rapid development of satellite technology, more satellites are providing more detailed aerosol optical products via instruments such as the Polarization and Directionality of the Earth's Reflectance instrument (POLDER) and the Ozone Monitoring Instrument (OMI) (Levelt, et al., 2006; Tilstra and Stammes, et al., 2007). However, the accuracy of the SSA and ASY products over China still needs to be improved (Oikawa et al., 2013; Dubovik, et al., 2019). Recently, using satellite and observational data assimilated into the Goddard Earth Observing System, version 5 (GEOS-5), the National Aeronautics and Space Administration (NASA) has extended the Modern-Era Retrospective Analysis for Research and Application, Version 2 (MERRA-2). Compared with its predecessor (MERRA-1), MERRA-2 offers important improvements in aerosol assimilations (Gelaro et al., 2017). The new dataset has the potential to provide improved estimates of aerosol microphysical parameters, such as SSA, and can be further used in the ADRF estimation. After SSA is determined, ASY, the only unknown model input, can be retrieved by matching the simulated radiative fluxes with satellite measurements from Clouds and the Earth's Radiant Energy System (CERES). Overall, based on the satellite and reanalysis datasets, including MERRA-2, the MODerate Resolution Imaging Spectroradiometer (MODIS) and CERES, the objective of this study is to provide quantitative estimates of fine-resolution ADRF distributions under clear skies using a radiative transfer model over eastern China (114–124^∘ E, 24–38^∘ N, shown in the Fig. 1). Additionally, the aerosol vertical profiles in each grid, which were not considered in previous studies, are used to obtain more accurate ADRF. In our study, aerosol vertical profiles are determined by the Weather Research and Forecasting Model (WRF, version 3.2.1) and the National Centers for Environmental Prediction Final Operational Global Analysis (NCEP FNL). The detailed algorithm of aerosol profiles can be found in Sect. 2. Other data acquisition is also presented in Sect. 2, and Sect. 3 introduces the method of ADRF simulations. Section 4 includes the retrieval of aerosol optical properties, validation of surface radiative fluxes with pyranometers, and detailed discussion of the error sources. Then this method is applied in each grid of eastern China during 2000–2016, and the uncertainty in the retrieval method is also discussed in Sect. 4. The conclusion is presented in Sect. 5.

2 Data

To acquire ADRF, the inputs (aerosol optical depth (AOD), SSA, ASY, albedo, etc.) to the radiative transfer model were determined from a combination of satellite and reanalysis datasets. AOD was derived from Collection 6 (C6) of MODIS Level 2 products over land (10 km resolution at the nadir) from the Terra satellite (Levy et al., 2013). MODIS AOD retrieval primarily employs three spectral channels, centered at 0.47, 0.66, and 2.1 µm and is interpolated at 0.55 µm (Kaufman et al., 1997). Li et al. (2003) demonstrated that the MODIS AOD Level 2 product is appropriate in eastern China and exhibits high precision. Compared with C5, MODIS C6 mainly updated the cloud mask to allow heavy smoke retrievals and fine-tuned the assignments for aerosol types as a function of season and location over the land. Levy et al. (2013) made a comparison between MODIS C5, C6, and AERONET and found that the correlation coefficient of C6–AERONET increases slightly, and the slope and offset of the regression curve only changed slightly compared with C5–AERONET. In addition, He et al. (2010) found that MODIS AOD was highly correlated with sun photometer (CE318) measurements at seven sites in the Yangtze River Delta (YRD) region (118–123^∘ E, 29–33^∘ N), with a correlation coefficient of 0.85 and with 90 % of cases falling in the range of ΔAOD = ± 0.05 ± 0.20 AOD (Chu et al., 2002). Thus, the uncertainty in the AOD is regarded as 20 % in this study.

The hourly SSA product was provided by MERRA-2. MERRA-2 combines GEOS-5 and the three-dimensional variational data assimilation (3D-Var) Gridpoint Statistical Interpolation (GSI) analysis system. GEOS-5 is coupled to the Goddard Chemistry Aerosol Radiation and Transport (GOCART) aerosol module, which includes five particulate species (sulfate, dust, sea salt, and organic and black carbon) (Colarco et al., 2010). The optical properties of these aerosols are primarily from the Optical Properties of Aerosols and Clouds (OPAC) dataset, in which aerosol optical parameters are calculated based on microphysical data (size distribution and spectral refractive index) under the assumption of spherical particles, and they are given for up to 61 wavelengths between 0.25 and 40 µm (Hess et al., 1998). MERRA-2 provides SSA data at 0.55 µm. These data are calculated by the ratio of total aerosol scattering aerosol optical thickness (AOT) to total aerosol extinction AOT at 0.55 µm, and these two are the outputs of the GOCART model (Colarco et al., 2010). More details of the aerosol module in MERRA-2 can be found in Randles et al. (2017) and Buchard et al. (2017). The new dataset has been used in many recent studies and is appropriate for environmental research (Song et al., 2018). The input SSA was interpolated to other wavelengths in SBDART, which will be discussed in detail in the Methodology (Sect. 3).

The upward radiative flux at TOA was used to constrain and determine the ASY. The shortwave (SW, 0.3–5 µm) TOA flux was acquired by the CERES Single Scanner Footprint (SSF) level 2 product from the Terra satellite. CERES SSF measures the instantaneous reflected SW radiance under clear-sky conditions. To convert from radiance to flux, angular distribution models (ADMs) were used in the CERES SSF product (Loeb et al., 2003). The CERES file contains 1 h of data, and the CERES SSF footprint nadir resolution is approximately 20 km. According to Su et al. (2015), the uncertainty of TOA SW flux is 1.6 % over clear land.

Another important parameter for ADRF simulations is the surface albedo, and it was derived from the daily MODIS MCD43C3 black-sky albedo product (C6). The surface albedo product includes seven narrow bands and three broadbands (visible (0.3–0.7 µm), near-infrared (0.7–5.0 µm), and SW (0.3–5 µm)). Here, albedo product in the SW band was used in our study. Each file contains 16 d of combined level 3 data from the satellites Aqua and Terra, with a spatial resolution of 0.05^∘. It also contains the data quality information, that is, the proportion of inversion retrieval information in each pixel. For example, a data quality index of 0 represents the best quality (100 % with full inversion and no fill values); this index increases with the decrease in the proportion of inversion retrieval pixel, and 4 represents 50 % or less of fill values. Notably, to ensure accuracy, only the albedo values with a high quality index (0–4) were used. The uncertainty in the high-quality MODIS albedo is less than 5 % (Cescatti et al., 2012).

The total column ozone, total column water vapor, and atmospheric profile data were from ERA-Interim (European Center for Medium-Range Weather Forecast (ECMWF) Interim Reanalysis). Specifically, the atmospheric profile includes the altitude, temperature, water vapor density, and ozone density at 37 pressure levels (1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100 to 250 at 25 hPa intervals; 300 to 750 at 50 hPa intervals; and 775 to 1000 at 25 hPa intervals). The data quality of the ERA-Interim reanalysis data can be found in Dee et al. (2011).

The aerosol vertical profile plays a non-negligible role in aerosol radiative forcing. In SBDART, the aerosol vertical profile is shaped by aerosol density and the according altitude. The aerosol density is a proportion of AOD in different altitudes, and the overall profile is scaled by AOD. The aerosol density is set to fall exponentially between two altitudes by default. In our study, the aerosol vertical profile in SBDART was derived from a two-layer aerosol vertical distribution model, which is proposed by He et al. (2008). In this two-layer aerosol model (Fig. S1 in the Supplement), aerosol extinction coefficient is assumed to decrease exponentially with altitude above the top of the planetary boundary layer (PBL), and the extinction coefficient keeps uniform below the PBL. Based on this aerosol model, two inputs of aerosol vertical profile need to be determined, PBL and aerosol layer height (ALH). ALH is defined as the level where the aerosol extinction coefficient decreases to 1/e (scaling height) of that at the top of the PBL. PBL and ALH input to SBDART along with the according aerosol density. In this study, PBL was simulated using a three-domain, two-way nested simulation of the WRF Model (version 3.2.1). ALH can be influenced by the transport of air mass and the convective dispersion of aerosols, both of which are usually associated with large-scale weather systems. Based on the different meteorological conditions, an automated workflow algorithm of ALH was constructed, and ALH was estimated by the meteorological parameters (relative humidity, temperature, wind speed, and wind direction) from NCEP FNL (Kalnay et al., 1996). The detailed algorithm and the according calculations of PBL and ALH retrieval can be found in He et al. (2016). The aerosol profiles were utilized to calculate the surface-level visibility from AOD, and the long-term spatial comparison with surface measurements over eastern China showed that 90 % of the samples exhibited correlation coefficients greater than 0.6 and that 68 % of the samples exhibited correlation coefficients greater than 0.7 (He et al., 2016).

Table 1 Satellite and reanalysis datasets used in the study.

All of these multiplatform datasets with their spatial and temporal resolutions were summarized in Table 1. In this study, bilinear interpolation was used in these datasets, and these datasets were interpolated to a spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ to collocate with the MODIS AOD data. The ADRF simulation was also performed in each $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grid over eastern China. For temporal resolution, AOD and TOA radiation fluxes were from the MODIS and CERES sensors aboard the Terra satellite, respectively, and they are available once per day. The temporal resolution of SSA and ERA-Interim is hourly, and surface albedo data are daily means. The ADRF simulations were only performed at the passing over of the Terra satellite under clear skies. The temporal coverage is from 2000 to 2016. The research area and surface measurement sites for validation are shown in Fig. 1.

Figure 2 A schematic diagram to simulate ADRF based on satellite and reanalysis datasets.

3 Methodology

Clear-sky ADRF in the SW (0.25–4 µm) spectral region was simulated by the Santa Barbara Discrete Atmospheric Radiative Transfer (SBDART) model (Ricchiazzi et al., 1998). This model has been widely adopted for the estimation of aerosol radiative forcing and validated with high accuracy (Li et al., 2010). In this study, the SBDART model was used to estimate broadband SW (0.25–4 µm) surface irradiances and ADRF over eastern China. It is based on the DISORT radiative transfer model, the low-resolution band models developed for LOWTRAN 7 atmospheric transmission, and the Mie scattering results for light scattering by water droplets and ice crystals (Ricchiazzi et al., 1998). Here, the LOWTRAN 7 (Low Resolution Atmospheric Transmittance 7) solar spectrum was adopted in SBDART. This radiative transfer model also includes the standard aerosol models derived from Shettle and Fenn (1975), in which aerosol optical parameters are wavelength dependent and the scattering parameters depend on the surface relative humidity. Users can also define different aerosol parameters in different wavelengths. The default of the corresponding spectral information is interpolated or extrapolated to all wavelengths using linear fitting on SSA or ASY and using Ångström coefficients on AOD. According to P. Wang et al. (2009), the input of aerosol parameters has a very minor effect on the accuracy of irradiance simulation when using spectrally averaged values compared with detailed spectral information. Therefore, aerosol parameters (AOD, SSA, ASY) at 0.55 µm were used in the radiative transfer model. As for surface albedo, it is simply assumed that angular distribution of surface-reflected radiation is completely isotropic in the model. In our study, the MODIS SW MCD43C3 (0.3–5 µm) product is used as albedo input, and it is nearly consistent with wavelength coverage (0.25–4 µm) of the output surface irradiances in SBDART.

As shown in Fig. 2, the main inputs of the SBDART model include aerosol properties (AOD from MODIS; SSA from MERRA-2; ASY from the retrieval (Sect. 4.2)), surface albedo (from MODIS), aerosol vertical profile (from NCEP), atmospheric profiles (from ECMWF), total column ozone, and water vapor (from ECMWF). The main outputs are radiative fluxes at the surface and TOA with and without aerosols. ADRF is defined as the difference in net radiative flux (downward minus upward) between aerosol and no-aerosol conditions. Here, we mainly concentrated on ADRF at the surface:

$$\begin{array}{}\text{(1)}& {\text{ADRF}}_{\text{sur}}=\left(F\downarrow -F\uparrow \right)-\left(F_{\mathrm{0}}\downarrow -F_{\mathrm{0}}\uparrow \right),\end{array}$$

where F and F₀ represent radiative fluxes with and without the aerosol at the surface, respectively. The upward and downward arrows denote the directions of the radiative fluxes, which can be obtained by the outputs of SBDART. For simplicity, the upward radiative fluxes at the TOA are called F_u_toa, and the downward/upward radiative fluxes at the surface are called F_d_sur and F_u_sur, respectively (see Table A1 for the acronyms).

In addition, Mann–Kendall (MK) test (Mann, 1945; Kendall, 1975) was used to calculate the trend of ADRF time series and its significance level (above 90 %) in our study. It identifies whether monotonic trends exist in a time series and is widely employed for trend analysis of aerosol data. The detailed analysis can be found in Li et al. (2014). Prior to trend analysis, ADRF data were deseasonalized by subtracting the monthly mean during 2000–2016 to eliminate the influence of annual and seasonal cycles.

Figure 3 (a) The location of six sun photometer sites over eastern China. (b) The scatter plots of SSA between MERRA-2 and the sun photometer in Xuzhou, Shouxian, and Hefei. Orange dots represent Xuzhou samples, and the orange line is the fitting curve of Xuzhou samples while green represents Shouxian and black represents Hefei. Dashed lines are the range of ±10 % relative error. (c) The scatter plots of SSA between MERRA-2 and the sun photometer in Taihu, Pudong, and Hangzhou. Red dots represent Taihu samples, and the red line is the fitting curve of Taihu samples while purple represents Pudong and yellow represents Hangzhou. Dashed lines are the range of ±10 % relative error.

Table 2 The geographical characteristics, observing period, and sample number of sun photometer sites. The fitted regression equations between MERRA-2 and sun photometer SSA are also shown here. In the equation, x represents the SSA sample, and y represents the fitted value of SSA.

4 Results and discussion

4.1 Retrieval of aerosol properties

Before ADRF simulation, the accuracy of the MERRA-2 SSA product was evaluated firstly. In eastern China, six sun photometer sites, Xuzhou (117.14^∘ E, 34.22^∘ N), Shouxian (116.78^∘ E, 32.56^∘ N), Hefei (117.16^∘ E, 31.91^∘ N), Taihu (120.22^∘ E, 31.42^∘ N), Pudong (121.79^∘ E, 31.05^∘ N), and Hangzhou (120.16^∘ E, 30.29^∘ N) (Fig. 3a), were chosen for comparison with MERRA-2 SSA data. The location of the sun photometers was shown in Fig. 3a, and their geographical characteristics, observing periods, sample numbers, and the fitted regression equation between MERRA-2 and sun photometer SSA were presented in Table 2. Five sites (Xuzhou, Shouxian, Hefei, Taihu, and Hangzhou) are AERONET sites, and level 1.5 inversion data of AERONET were used. The uncertainty of AERONET products can be found in Dubovik and King (2000). Another sun photometer (CE318, Cimel Electronique, France) in Pudong was calibrated annually and maintained routinely, and a detailed description of calibration was presented in Cheng et al. (2015). The sun photometer spectral products are available at wavelengths of 440, 675, 870, and 1020 nm, and they were interpolated at 0.55 µm to match MERRA-2 SSA. The collection time was constrained from 09:00 to 14:00 (local time), covering the overpass time of the Terra satellite. Meanwhile, the relatively high solar zenith in this period avoids possible inversion errors and improves the data accuracy (Tian et al., 2018b). Additionally, the specific MERRA-2 grid cell containing the sun photometer was selected, and sun photometer SSA was hourly averaged to match the MERRA-2 SSA product. The detailed comparisons at Xuzhou, Shouxian, and Hefei are shown in Fig. 3b. Figure 3c displays the comparison results at Taihu, Pudong, and Hangzhou. As shown in Fig. 3, dashed lines are the range of ±10 % relative error; all samples in Taihu, Pudong, Hefei, 94 % of samples in Xuzhou, 93 % in Shouxian, and 98 % in Hangzhou fall within the ±10 % error. This finding suggests that MERRA-2 SSA agrees well with the sun photometer data, even though a few SSA samples are beyond the error range. Furthermore, the slopes of the linear fitting curve are less than 1 at all sites except Shouxian (Table 2), and this reveals that MERRA-2 SSA has systematic biases in most areas of eastern China. The primary reason for the discrepancy is the simple aerosol model assumption in MERRA-2 (Buchard et al., 2017). Only five aerosol types (sulfate, dust, sea salt, and organic and black carbon) are involved; the lack of nitrate aerosols, which are highly scattering aerosols, may result in the underestimation of MERRA-2 SSA. In addition, the calibration errors among these instruments should be considered. Generally, the evaluation results in six sites showing that the accuracy of the MERRA-2 SSA product is acceptable in eastern China, with ±10 % uncertainty.

Figure 4 A detailed workflow of the binary search used in ASY retrieval.

After SSA was determined, ASY is the only unknown input parameter. ASY is the key to portraying the scattering direction of aerosols. ASY = 1 denotes completely forward scattering, and ASY = 0 is symmetric (Rayleigh) scattering. Here, gridded ASY was simulated by matching observed F_u_toa (from CERES) with simulated F_u_toa (from SBDART). The sensitivity test indicates that F_u_toa, just similar to F_u_sur (shown in Fig. S3b), is a monotonically increasing function of ASY with other fixed inputs. Consequently, only one F_u_toa can be obtained with one specific ASY. With this premise, a binary search was applied to approximate ASY to improve calculation efficiency (Chang, 2013). The goal of the binary search is to find the ASY when the simulated F_u_toa is close to the observed F_u_toa. To accomplish this, the ranges of F_u_toa are repeatedly diminished by taking the middle ASY as one of the boundary values, and when the difference between the F_u_toa observed by CERES and calculated by SBDART is less than 1, the corresponding approximation of ASY is finally obtained. A detailed scheme is illustrated in Fig. 4. First, the value for ASY is initially assumed in the reasonable range of 0.1–0.9, and the upper and lower boundaries of ASY, along with other parameters, are input to SBDART to yield the initial range of calculated F_u_toa_a and F_u_toa_b. Then, this range is checked to determine whether it includes the F_u_toa (observed by CERES) by multiplying ((F_u_toa_a-F_u_toa)^∗( F_u_toa_b-F_u_toa)). If the multiplication result is negative, meaning that ASY falls within this range (ASYa, ASYb), the average of F_u_toa_a and F_u_toa_b is set as a new boundary (F_u_toa_c). Otherwise, this case is discarded, and the retrieval is not continued (ASY = NaN), perhaps due to inappropriate inputs. Next, for cases in which the multiplication result is negative, the multiplication process is applied to the new boundary ((F_u_toa_a-F_u_toa)^∗(F_u_toa_c-F_u_toa)). If this multiplication result is negative, the ASY falls within this range (ASYa, ASYc). Then, ASYc is set to represent ASYb. Otherwise, ASYc is set to represent ASYa. This process represents the scope-narrowing of the ASY boundary discussed above. With several iterations of narrowing the scope, the boundaries of the simulated F_u_toa become close to the true value of F_u_toa (observed by CERES). When the difference between the simulated F_u_toa boundary and the observed F_u_toa is less than 1, the corresponding ASY is considered to be one approximation. In this process, the input parameters, including AOD (from MODIS), SSA (from MERRA-2), surface albedo (from MODIS), aerosol vertical profile (from NCEP), atmospheric profiles (from ECMWF), total column ozone, and water vapor (from ECMWF), were input into the SBDART together in every iteration. All these inputs from 2000 to 2016 were used to simulate ADRF in each grid of eastern China. All calculations were performed on the Linux system. Following this method, ASY was retrieved in each grid cell over eastern China. The range of retrieved ASY is 0.50–0.80, and the mean ASY is 0.63, which is consistent with the observation site (Taihu) in eastern China (Xia et al., 2007). According to Mie theory, ASY is determined by the size distribution and the complex refractive index of aerosols. Therefore, the difference of ASY in eastern China can be partly related to the difference of fine-mode radius. Xia et al. (2007) has reported that the fine-mode volume median radius at the Taihu site averages 0.181 µm over a range of AOD from 0.6 to 1.0, while it is 0.168 µm in northern China. In ASY retrieval, ASY is assumed to vary enough to match F_u_toa while ensuring the accuracy of all other inputs (e.g., AOD, SSA). This assumption can deviate from the reality if there are obvious differences between real and retrieval values of other inputs. This above condition can easily occur in the process of ASY retrieval, when ASY cannot be retrieved (ASY = NaN). Even if ASY can be obtained, ASY can be inaccurate when other inputs have large biases. The uncertainty of ASY caused by the other inputs (AOD, SSA, albedo, CERES F_u_toa) will be quantified in the following uncertainty analysis (Sect. 4.3).

Figure 5 The scatter plots between observed F_d_sur by pyranometers and simulated F_d_sur by SBDART in Baoshan, Fuzhou, and Yong'an. The blue line is the fitting curve and the dashed line represents y=x. The red dots denote the specific case in which the pyranometer captures the fluctuation of F_d_sur by clouds during 1 h. The olive green dots denote the specific case in which the site is completely covered by clouds, deduced from a MODIS true color map composed of channels 1, 4, and 3. The blue dots represent the other ordinary case.

After aerosol optical properties were obtained, these parameters from multiplatform datasets can be input into the SBDART model to simulate surface radiative fluxes and ADRF in eastern China according to the methodology in Sect. 3.

4.2 Validation of the method

Before conducting ADRF simulation in each grid of eastern China during 2000–2016, this method was first applied in the three grids of selected sites to assess the performance of ADRF retrieval. Three radiation sites in Baoshan (121.45^∘ E, 31.4^∘ N), Fuzhou (119.29^∘ E, 26.08^∘ N), and Yong'an (117.37^∘ E, 25.98^∘ N) were chosen to make the comparisons between calculated F_d_sur and surface observation by the pyranometers (FS-S6, China) during 2014–2016. Red circles in Fig. 1 denote the specific locations of pyranometers. Baoshan and Fuzhou are urban and coastal sites while Yong'an represents suburb and inland sites. The different aerosol concentration levels and abundant aerosol types in these sites can represent the most of aerosol properties in eastern China. These pyranometers had regular maintenance and were calibrated annually through intercomparisons with the basic-reference station. Additionally, quality control has been performed at these sites according to Long and Shi (2008), including the removal of physical possible limits as determined by the Baseline Surface Radiation Network (BSRN) and use of configurable limits based on climatological analysis of measurement data. The uncertainty in the pyranometers is expected to be 5 % (Song, 2013). Simulated F_d_sur was averaged in the scope of a 40 km side length with the center at the pyranometer, and the measured F_d_sur was averaged within ±30 min of the satellite overpass (Ichoku, et al., 2002).

Figure 6 Yearly mean ADRF distributions during 2000–2016 over eastern China (unit: W m⁻²).

Figure 5 displays the comparison results between simulated F_d_sur and observed F_d_sur by pyranometers at three sites. The simulated F_d_sur is fairly consistent with the observations, with correlation coefficients of 0.87 in Baoshan (Fig. 5a) and Fuzhou (Fig. 5b) and 0.90 in Yong'an (Fig. 5c). Root-mean-squared error (RMSE) is a good indicator for measuring the discrepancy between observed and simulated F_d_sur data. The RMSE is 7.9 W m⁻² in Baoshan, 7.5 W m⁻² in Fuzhou, and 5.6 W m⁻² in Yong'an. This discrepancy only accounts for 3 %–5 % of ADRF, indicating that this retrieval method has a relatively higher accuracy than those in other studies (e.g., Thomas et al., 2013; Fu et al., 2017). Additionally, all slopes are less than 1, which implies that the method has systematic biases at these sites. A similar tendency was found in the comparison between MODIS AOD and sun photometers in eastern China by He et al. (2010); it is speculated that the main systematic error in ADRF simulation may come from the input, MODIS AOD. Nevertheless, satisfactory comparison results indicate the suitability and feasibility of ADRF retrieval in or near the sea and at urban–suburban sites of eastern China, although the type of underlying surface and aerosol properties is evidently different in these areas.

Figure 7 Averaged spatial distribution of (a) ADRF (unit: W m⁻²), (b) AOD, and (c) SSA during 2000–2016 in eastern China.

To further assess the discrepancy between simulated F_d_sur and the observations, the relative errors of each case at the three sites were calculated. The results suggest that underestimated cases (negative relative errors) account for 61 % of the total cases, and overestimated cases (positive relative errors) account for 39 %. According to the validation results, the sources of error in the simulation may be attributed to the following reasons.

Cloud contamination. An examination of cloudiness was carried out at the three sites. According to the empirical clear-sky detection method, 1 h radiation data of a pyranometer was used to discriminate clear-sky observations (Xia et al., 2007). The red dots in Fig. 5 represent the cloudiness case detected by the pyranometer. Meanwhile, from the MODIS true color map composed of channels 1, 4, and 3 (not shown), the olive green dots denote the specific case in which the site is completely covered by clouds. Taking one olive green case (Baoshan, 18 October 2014) as an example. As shown in Fig. S2, it is obvious that a large amount of cloud exists in the area of 29–31^∘ N and 120–122^∘ E, and the Baoshan site is at the edge of the cloud. In this case, MODIS AOD was overestimated compared with sun photometer AOD, because some cloud effects were not completely removed from the MODIS AOD calculation. Therefore, a large discrepancy can occur in these cases between simulated F_d_sur and observation. The cloud effect, especially residual thin cirrus clouds, is difficult to completely remove from MODIS AOD (Kaufman et al., 2005). Moreover, the cloud mask algorithm in the MODIS aerosol inversion sometimes fails to distinguish fog or haze in high-humidity conditions. Many more fog days can be observed in Fuzhou than the other two sites, and fogginess can significantly reduce the accuracy of the simulation (Ye et al., 2010). In addition, the error source of MODIS AOD is also from errors in the aerosol model assumption and surface reflectivity (Xie et al., 2011).

Different spatial and temporal representativeness. In the validation, the area measurement (satellite and reanalysis data) was compared to point measurements (pyranometer). For temporal matching, the pyranometer can capture the process of perturbation induced by air mass movement within 1 h, whereas satellites can only provide the instantaneous condition. Hence, this comparison method inevitably introduces some degree of uncertainty.

Instrument and radiative transfer errors. One error source in pyranometers is the thermal offset effect. This spurious signal is due to the difference in temperature between the inner dome and the detector of a pyranometer and can lead to additional errors in the irradiance measurements, especially diffuse irradiance (Sanchez et al., 2015). To reduce this effect, a pyranometer should be installed in a transparent ventilation hood. Alternatively, several statistical methods have also been proposed to suppress the thermal offset effect (e.g., Song, 2013; Cheng et al., 2014). In this study, the correction of the thermal offset was not performed because of the lack of additional observation data. Aside from the instrument error, the model simulation discrepancy also depends on the radiative transfer models. They are based on some simplifications, including the sphericity of aerosol particles and the directional reflectance of the surface. Derimian et al. (2016) found that neglecting aerosol particle nonsphericity can overestimate the aerosol cooling effect. Furthermore, simulation results vary slightly among different models due to their different assumptions in radiative transfer. For instance, Yu et al. (2007) compared three models (second simulation of the satellite signal in the solar spectrum (6S), MODerate resolution atmospheric TRANsmission (MODTRAN), and SBDART) at Xianghe station and showed that approximately 80 % of the cases simulated by SBDART were lower than the surface observations, while the 6S simulation results were higher.

Figure 8 Time series of monthly mean ADRF (blue) and AOD (red) in eastern China from 2000 to 2016. These data are deseasonalized. Dashed lines represent the Mann–Kendall (MK) fitting trend of ADRF and AOD.

4.3 Sensitivity test and uncertainty analysis

To determine the uncertainty of the method for ADRF simulation caused by each input parameter, a sensitivity test for input parameters was carried out. A specific case in Shanghai on 11 October 2015, was used with the following values: AOD = 0.62, SSA = 0.85, ASY = 0.69, surface albedo = 0.13, total column water vapor = 0.69 g cm⁻², and total column ozone = 0.28 atm cm⁻¹. Figure S3 portrays the responses of F_d_sur, F_u_sur, and ADRF to changes in one parameter while holding the other parameters constant. To remove the impact of units, all the parameters are dimensionless; that is, the ratio of the input to the actual value is used as the x-axis value. The absolute value of every slope describes the impact of every parameter on the dependent variables (F_d_sur, F_u_sur, and ADRF). Figure S3 presents the actual condition of this case when the value of the x axis equals 1, in which F_d_sur is 629.15 W m⁻², F_u_sur is 83.52 W m⁻², and ADRF is −149.39 W m⁻². This situation denotes a strong cooling effect of aerosols at the surface. Apparently, different parameters impose diverse influences on the radiative values (F_d_sur, F_u_sur, and ADRF). As depicted in Fig. S3, AOD, SSA, and ASY are three crucial parameters that greatly influence F_d_sur. P. Wang et al. (2009) conducted the radiative closure experiment in the Netherlands and further found that AOD can affect the changes of direct and diffuse irradiation, while SSA and ASY only affect the diffuse irradiance. For F_u_sur, albedo, AOD, and SSA are more important parameters. The impact of surface albedo is much larger than the others because albedo actually determines how much of the irradiance is reflected by the surface. For ADRF, SSA, AOD, and ASY are major factors in determining ADRF. Additionally, AOD enhancement can increase the aerosol cooling effect at the surface, whereas increases in SSA and ASY can result in decreases in the aerosol cooling effect. In general, a sensitivity test shows that ADRF depends highly on AOD, SSA, ASY, and albedo. Two parameters (atmospheric profile and aerosol vertical profile) are not discussed because these parameters have little impact on clear-sky ADRF in the above case. The atmospheric profile has a minor effect on the perturbations of ADRF compared with the total columns of the atmospheric component (water vapor and ozone). This result has also been proven by Yu et al. (2007) and Li et al. (2016). As for the aerosol profile, two typical shapes were input to SBDART for the sensitivity test. The first type (type I) has an elevated aerosol layer, and the second type (type II) is the two-layer aerosol model as mentioned above (Fig. S1). The changes of the elevated layer height (type I) or PBL–ALH (type II) have very little impact on ADRF, and the according maximum value of the ADRF difference can only reach 0.5 W m⁻². This conclusion is consistent with Guan et al. (2009). However, this impact becomes much stronger in the presence of absorbing aerosols, especially in some extreme cases such as dust storms and biomass burning (Wang and Christopher, 2006). Reddy et al. (2013) also demonstrated that surface aerosol radiative forcing can be enhanced by 25 % due to the insertion of the extinction profile of absorbing aerosols to replace the default profile.

Figure 9 The spatial distribution of the ADRF trend in eastern China during 2000–2016 (unit: W m⁻² month⁻¹). Hatched regions represent those exceeding the 90 % significance level.

Table 3 Errors induced by different input parameters in ASY, radiative flux (F_d_sur, F_u_sur), and ADRF. Here, the uncertainties of input parameters (AOD, albedo, CERES F_u_toa) are from literatures and the uncertainty of SSA is from validation in Sect. 4.

^a He et al. (2010). ^b Cescatti et al. (2012). ^c Su et al. (2015).

On the basis of these four high-sensitivity factors, the uncertainties in ASY and ADRF due to these parameters were quantitatively assessed. According to data uncertainty mentioned in Sect. 2 and the SSA validation, the relative errors of AOD, SSA, albedo, and CERES F_u_toa are 20 %, 10 %, 5 %, and 1.6 %, respectively. This lower or upper limit of parameter errors was input to the ADRF calculation, and the associated uncertainty was calculated by the difference between the simulated radiative flux with parameter errors and without errors. Notably, the uncertainty analysis is based on extreme conditions, and the associated errors are much larger than the actual values. As displayed in Table 3, the uncertainty in ASY induced by SSA can reach up to 23 %, indicating that SSA is a decisive factor in ASY retrieval when using the CERES F_u_toa constraint. SSA also has the largest effect in regulating aerosol radiative forcing, which is consistent with the research on dust aerosols by Huang et al. (2009). AOD contributes uncertainties of 3.7 % in ASY and 15.4 % in ADRF. Albedo introduces −3.7 %–1.7 % uncertainty in ASY and approximately 3 % in ADRF. The error of the CERES product produces approximately 1.7 % uncertainty in ASY and 1.5 % in ADRF. The results of uncertainty analysis agree well with those of previous studies. For example, Xia et al. (2016) revealed that AOD and SSA together can account for 94 % of surface ADRF. Zhuang et al. (2018) further noted that the error sources from the absorbing component of AOD and coarse-aerosol SSA contributed to the greater uncertainty in the ADRF. Therefore, improving the precision of the input parameter is helpful for obtaining reliable ADRF estimation. As Michalsky et al. (2006) demonstrated, when using high-quality measurements as inputs to the model, the biases between modeled and measured irradiance can decrease to 1.9 %. In addition to these factors, Wang and Martin (2007) also revealed the effects of aerosol hygroscopicity on the aerosol phase function and the increase in SSA with RH enhancement, suggesting that relative humidity (RH) is also closely related to ADRF.

4.4 Long-term ADRF retrieval in eastern China

The above evaluations show the method for ADRF simulation is feasible and has high accuracy in eastern China; thus this method was further applied in each grid cell of eastern China to obtain a full coverage of ADRF during 2000–2016. Figure 6 outlines an overall picture of annual mean ADRF at the surface over eastern China during the past 17 years. It provides valuable information about aerosol radiative effect not only in the urban areas with intensive human activities, but also in the suburb with unavailable observational data. ADRFs in all grids are negative, ranging from −220 to −20 W m⁻², implying that aerosols have a cooling effect at the surface over eastern China. The yearly mean ADRF is −100.21 W m⁻². The magnitude of ADRF is higher than most cities in the world, such as Spain (Esteve et al., 2014), Gasan (Kim et al., 2006) and Karachi (Alam et al., 2011). The main reason is that AOD in eastern China is much larger than these cities, since eastern China has experienced rapid urbanization and economic development in the past 17 years and a robust increase can be found in anthropogenic emissions. For example, mean AOD in eastern China is 0.62 in this study during 2003–2011 while AOD is 0.19 in Spain during 2003–2011 (Esteve et al., 2014). The red area denotes the high absolute value of ADRF (Fig. 6), which is found in the densely populated and industrialized areas, including western Shandong Province, YRD, and Poyang Lake Plain. A low value (blue area) is observed in the southern part, such as Fujian and southern Zhejiang Province. An obvious difference of ADRF distributions is found between the northern and southern parts of eastern China, and the magnitude of ADRF increases from south to north. This pattern is consistent with site observations in Che et al. (2018), in which surface ADRF ranges from −150 to −100 W m⁻² at the northern sites of eastern China (Huainan and Hefei in Anhui Province) while ADRF ranges from −100 to −50 W m⁻² at the southern sites of eastern China (Jiande, Chunan, and Tonglu in Zhejiang Province). To further explore this difference, eastern China was divided into two parts: the north and south, with the boundary of 30^∘ N. The occurrence frequencies of annual ADRF for each grid cell in the north and south were calculated in Fig. S4. The occurrence frequency shows a broad range from −300 W m⁻² to 0 and the interval is 20 W m⁻². In the north, the largest proportion of ADRF, with the value of 76.47 %, falls in the range of −100 to −80 W m⁻², while the largest proportion (64.71 %) of ADRF falls in the range of −60 to −40 W m⁻² in the south. The extreme value over −250 W m⁻² may result from severe haze in the winter. Aerosol cooling radiative effect can sharply increase with large aerosol loadings. According to Yu et al. (2016b), surface ADRF can reach up to −263 W m⁻² on the haze days, while on the non-haze days, it can decrease to −45 W m⁻² in Beijing on January 2013. Usually in the heavy haze, the enhanced surface cooling, combined with atmosphere heating, can result in a more stable environment. It is unfavorable for the diffusion and dispersion of the aerosols and can further cause air accumulation and enhance aerosol ADRF (Wu et al., 2016). Meanwhile, positive ADRF was also found in a few grid cells, although it is not shown in Fig. S4. This condition occurs over bright surface in eastern China, especially with the abundance of absorbing aerosols (Sundström et al., 2015). According to the uncertainty analysis, ADRF is closely associated with the inputs (SSA and AOD). Based on this, comparison was conducted among the mean spatial distribution of ADRF, AOD, and SSA during 2000–2016 (Fig. 7). It is clear to see that the ADRF pattern is very similar to the negative phase of the AOD pattern; that is, the areas of high AOD have low ADRF. As for SSA, the higher value can be found in the south than the north, which indicates the aerosols in the south generally scatter more than in the north. Therefore, the large difference between north and south can be mainly attributed to the difference in AOD. The industry locations and topography between the north and south are obviously different. With the development of economy and urbanization, large amounts of anthropogenic aerosols in the north have imposed a strong cooling radiative effect in the past 2 decades. It is worth noting that, although western Shandong has lower urbanization compared with YRD, the aerosol cooling effect in western Shandong is even larger than in YRD. This is because Yimeng mountain (all places mentioned are shown in Fig. 1) located in the middle of Shandong blocks the west flow, leading to the enhancement of the aerosol accumulations and high AOD near its western border (He et al., 2012b). Meanwhile, Shandong is also easily impacted by air pollution transported from northern China. In addition, the high absolute value of ADRF is also found in Poyang Lake in Jiangxi with an abundance of anthropogenic aerosols, and these areas are surrounded by the mountains, where the poor ventilation conditions enhance aerosols. Compared with the north, the south is characterized by more extensive vegetation coverage and less human activities, and AOD is relatively lower in the south (Fig. 7b) and aerosols have a weaker cooling effect.

Apart from spatial changes, temporal changes of ADRF during 2000–2016 were also analyzed. Figure 8 displays the time series of monthly mean ADRF and AOD. For comparison, the blue line represents ADRF and the red line denotes AOD. They both show a fluctuation pattern, and they have an obvious negative phase with the correlation coefficient of −0.72. This indicates that the temporal change of ADRF is mainly attributed to the change of AOD. MK trends of ADRF and AOD are both positive but insignificant at the 90 % confidence level; in particular for the AOD trend, the value is close to zero. It shows AOD and ADRF did not change significantly during 2000–2016 in eastern China. Paulot et al. (2018) also proved this insignificant trend of ADRF in China based on chemical-climate models. Concerning AOD, Zhang et al. (2017) found that the AOD trend increased in 2000–2007 and then decreased in eastern China based on satellite observations. It is well known that the changes of AOD are closely linked with the change of anthropogenic emissions, especially in developing countries. Che et al. (2019) calculated that SO₂ has been the dominant anthropogenic emissions factor for AOD in China during the past few decades. Further, model simulations also indicate the changes of sulfate aerosols are the largest contributor to AOD and aerosol effect in China (Paulot et al., 2018). MK trends of monthly mean ADRF in each grid cell during 2000–2016 were also calculated (Fig. 9). Hatched regions indicate those exceeding the 90 % significance level. A high positive trend can be found in Anhui and Jiangxi, indicating the aerosol cooling effect is weaker in this region during 2000–2016. However, a few regions experience this cooling effect more strongly, especially in the northeast and southern areas of Yimeng mountain in Shandong. In general, the changes of ADRF during the past 17 years are mainly due to the anthropogenic emissions in eastern China. In addition, Paulot et al. (2018) further pointed that there is a nonlinear relationship between anthropogenic emissions and AOD–ADRF when considering the mix and oxidation of different emissions.

5 Conclusions

In this study, based on multiplatform datasets, high-accuracy ADRF distributions over eastern China during 2000–2016 were portrayed. MERRA-2 SSA data were first compared with sun photometer data (Taihu, Xuzhou, Pudong), and the validation results show that the relative error of the MERRA-2 SSA is ±10 % over eastern China. Then, ASY in each grid was retrieved by matching the simulated F_u_toa by SBDART with satellite observations. Then, aerosol optical properties (AOD from MODIS, SSA from MERRA-2, and ASY from the retrieval), surface albedo (from MODIS), aerosol vertical profile (from NCEP), atmospheric profiles (from ECMWF), total column ozone, and water vapor (from ECMWF) served as input parameters for SBDART to simulate ADRF in each grid cell of eastern China during 2000–2016. The validation result of this method at three sites (Baoshan, Fuzhou, and Yong'an) reveals that simulated F_d_sur is highly correlated with the pyranometer data during 2014–2016, with correlation coefficients of 0.87 in Baoshan and Fuzhou and 0.90 in Yong'an. The RMSEs are 7.9 W m⁻² in Baoshan, 7.5 W m⁻² in Fuzhou, and 5.6 W m⁻² in Yong'an. It shows that ADRF retrieval is feasible and has high accuracy over eastern China. In addition, associated factors, including cloud contamination, instrument and radiative transfer errors, and different spatial and temporal representativeness, were confirmed to produce additional uncertainty in ADRF simulations. A sensitivity test shows that ADRF highly depends on AOD, SSA, ASY, and albedo. Uncertainty analysis shows the uncertainty in ADRF retrieval induced by SSA is calculated to be 24 % and that by AOD is 15.4 %. Finally, ADRF simulation was conducted in each grid of eastern China during 2000–2016. Long-term ADRF distributions over eastern China were presented for the first time. ADRFs in all grids are negative and the range of ADRF is between −220 and −20 W m⁻², implying that aerosols have a cooling effect on the surface over eastern China. Aerosols are found to have a stronger cooling effect in the north compared with the south. The ADRF spatial pattern is consistent with the negative phase of the AOD pattern, and the temporal changes of ADRF also have a close relationship with AOD. They indicate that the changes of ADRF in eastern China can mainly be attributed to the changes of AOD. Furthermore, the spatiotemporal changes of AOD and ADRF are controlled by anthropogenic emissions, especially sulfate emissions in eastern China during the past 17 years.

In summary, this study suggests that the method for ADRF retrieval is feasible in eastern China. Especially in suburbs with no monitoring resources, our study offers valuable information on the direct radiative impact of aerosols. It is noted that, in our study, ADRF was calculated during the time that a satellite passes over rather than the whole day. More additional observation data from the sites are needed to further verify the performance of the ADRF retrieval and constrain these multiplatform datasets to improve the ADRF accuracy. In addition, it is necessary to improve the satellite instruments and the retrieval algorithm of aerosol properties; more novel methods, such as machine learning, can be involved in the ADRF estimates (Yin, 2010; Yu and Song, 2013). In the future work, aerosol-induced changes in the surface radiation under climate change and agricultural economic impact will be studied. This work can provide a deep understanding of aerosol radiative effects and is also helpful for aerosol modeling over eastern China.