Planetary albedo (PA; shortwave broadband albedo) and its long-term variations, which are controlled in a complex way by various atmospheric and surface properties, play a key role in controlling the global and regional energy budget. This study investigates the contributions of different atmospheric and surface properties to the long-term variations of PA based on 13 years (2003–15) of albedo, cloud, and ice coverage datasets from the Clouds and the Earth’s Radiant Energy System (CERES) Single Scanner Footprint edition 4A product, vegetation product from Moderate Resolution Imaging Spectroradiometer (MODIS), and surface albedo product from the Cloud, Albedo, and Radiation dataset, version 2 (CLARA-A2). According to the temporal correlation analysis, statistical results indicate that variations in PA are closely related to the variations of cloud properties (e.g., cloud fraction, ice water path, and liquid water path) and surface parameters (e.g., ice/snow percent coverage and normalized difference vegetation index), but their temporal relationships vary among the different regions. Generally, the stepwise multiple linear regression models can capture the observed PA anomalies for most regions. Based on the contribution calculation, cloud fraction dominates the variability of PA in the mid- and low latitudes while ice/snow percent coverage (or surface albedo) dominates the variability in the mid- and high latitudes. Changes in cloud liquid water path and ice water path are the secondary dominant factor over most regions, whereas change in vegetation cover is the least important factor over land. These results verify the effects of atmospheric and surface factors on planetary albedo changes and thus may be of benefit for improving the parameterization of the PA and determining the climate feedbacks.
The planetary albedo of Earth, which is the ratio of shortwave radiation reflected by Earth to the incoming shortwave radiation at the top of the atmosphere, may be considered a key parameter in regulating the global climate system and its variability due to its substantial role in controlling the global energy budget and surface temperature (Stephens et al. 2015; Wielicki et al. 2005; Donohoe and Battisti 2011). For example, a small change of 0.01 in the planetary albedo corresponds to a change in the shortwave net flux of 3.4 W m−2, which would approximately compensate for the radiative forcing of double the amount of CO2 in the atmosphere (IPCC 2001; Wielicki et al. 2005; Bender et al. 2006). Moreover, Budyko (1969) found that small variations in the planetary albedo could be sufficient for the development of Quaternary glaciations. Although many studies have indicated that the current albedo of Earth maintains a relative stable value (approximately 0.29) and displays a high degree of hemispheric symmetry (e.g., IPCC 2001; Loeb et al. 2009; Voigt et al. 2013; Stephens et al. 2015), our understanding of albedo remains limited owing to an incomplete understanding of its underlying physical processes. Thus, existing models exhibit relatively large discrepancies among simulations of the regional planetary albedo (e.g., Halthore et al. 2005; Wild 2005). For example, when analyzing albedo simulations from phase 3 of the Coupled Model Intercomparison Project (CMIP3), Donohoe and Battisti (2011) found that the intermodel spread in albedo was predominantly due to the differences in the atmospheric albedo among the different models. Therefore, a reasonable evaluation of the impacts of different feedback processes (e.g., cloud and surface properties) on the long-term albedo variability would be very helpful for improving our understanding of regional climate change and predicting how the planetary albedo will respond to climate changes.
Generally, changes in the planetary albedo are complicatedly controlled by both atmospheric and surface properties, for example, the cloud fractional coverage, cloud liquid water path, cloud ice water path, water vapor, and atmospheric aerosol amounts (Klein and Hartmann 1993; Held and Soden 2000; Christopher and Zhang 2002; Bender et al. 2006; Guo et al. 2011; Engström et al. 2015; Zhao and Garrett 2015; Lin et al. 2013; Xie et al. 2017), in addition to vegetation cover, land use, desertification, and snow and ice cover (Myhre and Myhre 2003; Zeng and Yoon 2009; Barnes and Roy 2008; Kashiwase et al. 2017). Joint influences from these factors and their interactions can complicate the simulations and predictions of planetary albedo. Furthermore, uncertainty in scene type can also complicate the observed flux. Until now, many efforts have been made to lessen the intermodel spread in albedo model simulations, based on satellite observations or in situ measurements (Qu and Hall 2005; Wang et al. 2006; Bender et al. 2006; Kato et al. 2006; Engström et al. 2015; Stephens et al. 2015; Loeb et al. 2016; Bender et al. 2017). For example, Donohoe and Battisti (2011) used the Clouds and the Earth’s Radiant Energy System (CERES; Wielicki et al. 1996; Corbett and Loeb 2015) flux data to quantify the relative contributions of the surface and atmosphere to planetary albedo on a global scale. Their results indicated that most of the observed planetary albedo is caused by atmospheric reflection and that the atmosphere attenuates the surface contribution to the planetary albedo. By using the International Satellite Cloud Climatology Project (ISCCP) D-series flux dataset, Qu and Hall (2005) found that the interannual variability in the planetary albedo within ice- and snow-covered regions is mainly attributable to variations in the surface albedo, but the atmospheric processes can attenuate 90% of the surface albedo effect on changes in planetary albedo. Pistone et al. (2014) pointed out that changes in cloudiness appear to play a negligible role in the observed Arctic darkening and that cloud albedo feedback may not be effective in offsetting Arctic warming. In addition, several other studies have focused mainly on the impacts of different atmospheric and surface properties on the components of planetary albedo (e.g., cloud and surface albedo). For example, as considerable portions of the surface of Earth are becoming greener as a result of climate change, causing a rise in CO2 concentration and nitrogen deposition (Piao et al. 2015), Forzieri et al. (2017) used the leaf area index (LAI) to study the effects of vegetation changes on the local climate and found that an increase in the LAI contributed to a reduction in the surface albedo. In addition, some studies have also shown that the macrophysical properties of clouds (e.g., the cloud fractional coverage and cloud liquid water path) dominate the atmospheric albedo (Stephens 2005) and that changes in the cloud fraction dominate changes in planetary albedo (Loeb et al. 2007; Bender et al. 2017). However, most of these studies are limited to specific locations, short investigation periods, or specific contributory variables. Systematic studies concerning the statistical relationships between long-term variations in the planetary albedo and different contributory variables at the regional scale have received far less attention.
To better understand the long-term variations in the regional planetary albedo, two key questions must be addressed in this investigation. First, what are the factors that drive the temporal variability in the planetary albedo at the regional scale? Second, which one of these factors is more important? In the following study, we will use multiple satellite datasets to build a regression relationship between planetary albedo and various variables to further quantify the relative contributions of different factors to the observed variability in planetary albedo. This paper is organized as follows. A brief introduction to all of the datasets and methods used in this study is given in section 2. Section 3a describes the global characteristics of planetary albedo and the difference between Aqua and Terra. Further analysis regarding the impacts of atmospheric and surface parameters on long-term variations in planetary albedo and the contribution evaluations are provided in section 3b. Finally, the conclusions and discussion are presented in section 4.
2. Datasets and methods
In the following study, 13 years (from 2003 to 2015) of data from multiple satellite datasets are collected to analyze the impacts of different factors on the long-term variability in regional-scale planetary albedo (PA).
a. Terra and Aqua
Terra was launched on 18 December 1999 and placed into a near-polar, sun-synchronous orbit at an altitude of 705 km with a 1030 local time (LT) descending node. Complementary to Terra, Aqua was launched on 4 May 2002, with a 1330 LT ascending node. Here, we use the products from two instruments (CERES and MODIS) carried by the Terra and Aqua satellites to provide the monthly mean radiative flux, cloud properties, and surface vegetation index.
The CERES instrument can accurately measure the top-of-atmosphere (TOA) radiances and convert radiances into fluxes via the use of angular dependence models (Su et al. 2015a,b). The instantaneous fluxes are converted to daily mean fluxes using empirical diurnal albedo models based on the cloud, atmosphere, and surface conditions at the time of the observation (Doelling et al. 2013; Loeb et al. 2018). In the following study, the TOA outgoing SW flux (all sky) in the CERES Single Scanner Footprint 1.0° (SSF1deg) monthly edition 4A (Ed4A) dataset is used to calculate the all-sky albedo. For clear-sky albedo, we use the clear-sky TOA SW fluxes in the CERES Energy Balanced and Filled (EBAF) monthly Ed4A dataset instead of CERES SSF monthly product. This is because the EBAF product combines both Terra and Aqua clear-sky measurements, and it also includes CERES subfootprint flux clear-sky measurements and thus greatly increases the sampling over those persistent cloud domains (e.g., the Southern Ocean). However, the TOA outgoing SW flux (clear sky) from CERES SSF1deg product is only used to study the Terra and Aqua cloud albedo forcing differences in the section 2c. The CERES SSF1deg product assumes that there is a solar zenith angle dependency of the clear-sky albedos, and it also assumes that the solar zenith angle dependency is symmetric about noon. In addition, note that the SSF product is a combination of CERES radiation data and coincident cloud properties from MODIS measurements (Sun et al. 2006; Zhan and Davies 2016). Thus, this product includes some cloud parameters [e.g., the daytime cloud area fraction (CF), ice water path (IWP), and liquid water path (LWP)], which can be used as cloud variables to assess the impacts of different cloud properties on the planetary albedo. It is worth noting that the time-averaged cloud parameters in the Ed4 are weighted by the cloud fraction. For example, the liquid and ice cloud properties were temporally weighted by the corresponding liquid or ice fractions to determine the daily or monthly mean (optical depth is the log of the optical depth and cloud fraction weighted). This ensures that the CF and LWP (or IWP) are the independent predictors of each other.
In addition, the CERES SSF dataset also provides the surface ice/snow percent coverage (I/SPC) data from the National Snow and Ice Data Center (NSIDC) (Nolin et al. 1998). This parameter is particular important at high latitudes and over the Tibetan Plateau, where snow and ice significantly affect the surface albedo (Pistone et al. 2014).
In addition to the snow/ice coverage, the local vegetation cover is also one of the most important factors in evaluating the variations in the surface albedo (Betts 2000; Sandholt et al. 2002). Here, we utilize the MODIS normalized difference vegetation index (NDVI) monthly product (MYD13C2.006) to describe the surface vegetation cover with a spatial resolution of 0.05° × 0.05°. The NDVI has been proven to be a robust indicator of the terrestrial vegetation productivity; it exhibits a sensitive response to vegetation dynamics and thus is considered as a useful tool for effectively reflecting the vegetation cover (Tucker et al. 2005; Beck et al. 2006). The quality of MODIS NDVI data is greatly improved as a consequence of the narrow MODIS red and NIR bandwidths (Huete et al. 2002). The definition of the NDVI is as follows:
where RNIR and Rred are the reflectances in the near-infrared (NIR) and red bands, respectively. Higher index values typically indicate a wider vegetation cover in a pixel.
b. CLARA-A2 dataset
In this investigation, we also use the monthly surface albedo information with a spatial resolution of 0.25° × 0.25° from the Cloud, Albedo, and Radiation dataset, version 2 (CLARA-A2) satellite product. The CLARA-A2 product is retrieved by the Advanced Very High Resolution Radiometer (AVHRR) operated onboard polar-orbiting NOAA satellites as well as by the MetOp polar-orbiting meteorological satellites operated by EUMETSAT. Compared with the CLARA-A1 (version 1) product, the CLARA-A2 product significantly enhances the quality of the surface albedo estimates through several important improvements (e.g., dynamic aerosol optical depths are used instead of a constant AOD, and wind speed data are used to describe the sea surface roughness and to retrieve the sea surface albedo) (Karlsson et al. 2017). The improved surface albedo product has been validated against in situ albedo observations and compared with a MODIS product (MCD43C3). This comparison showed that the CLARA-A2 surface albedo product is in good agreement with the surface albedo product from MODIS (MCD43C3) at the global scale and that the differences in the albedo estimates are less than 5% (Karlsson et al. 2017). Compared with the MCD43C3, an obvious advantage of the CLARA-A2 product is that albedo information over water bodies can also be derived. In addition, by checking time consistency (e.g., exclusively choosing afternoon satellite), the CLARA-A2 product achieves relatively homogeneous observation conditions to reduce the influence of orbital drift. And a major effort has been made to correct and homogenize the basic AVHRR level-1 radiance record (Karlsson et al. 2017).
In this study, the cloud albedo forcing αcloud is roughly estimated as follows:
where αall-sky and αclear-sky are the planetary albedo values under all-sky and clear-sky conditions, respectively. Here, it is worth noting that the cloud albedo forcing αcloud also includes some noncloud effects (e.g., aerosol direct radiative forcing) (Erlick and Ramaswamy 2003). In addition, the accurate estimation of the clear-sky two-way atmospheric transmittance T2 is rather difficult to obtain based on current datasets used in our study. To do this, we define the clear-sky reflected shortwave as
where S is the incoming shortwave at top of atmosphere and αsurface is the surface albedo. If αclear-sky is considered as
then T2 can be roughly approximated as
This parameter represents the atmospheric extinction ability to solar radiation when radiation from the TOA arrives at the surface and is then reflected back to the TOA. The impacts of water vapor and aerosols on the radiation are already included in the T2.
The planetary albedo is known to be affected by the variations in surface properties, which are determined by the land use (e.g., urbanization) as well as changes in the snow/ice and vegetation covers (e.g., desertification, forest cover changes) (National Research Council 2005; Wang et al. 2006). In this study, the I/SPC and NDVI (i.e., only for land areas, except for in polar regions) are used as proxies for representing the surface properties. Note that the surface albedo (SA) is considered a surface parameter only if both the NDVI and I/SPC fail to pass the t test. The atmospheric parameters include cloud properties and the atmospheric two-way transmittance T2. Recent studies have shown that cloud variability, especially the macrophysical cloud properties (i.e., the LWP and CF), dominate the variability of PA (Stephens et al. 2015; Seinfeld et al. 2016). As a result, the LWP, IWP, and CF are considered as cloud parameters in the following analysis. In addition, by performing collinearity diagnostics, we assess the strength and sources of collinearity among these variables at each grid and find that the predictors used in our study are almost noncollinear over all the regions (see Fig. S1 in the online supplemental material).
Because the stepwise regression method can effectively filter the predictors with collinearity and remove insignificant variables, we use this method to perform a multilinear regression analysis in each grid to construct a stable relationship between the planetary albedo anomalies and predictors and for further assessing the contributions of different variables to the planetary albedo anomalies. Figure 1 provides the valid sample number of the monthly parameters during 2003–15. (Note that Fig. 1 is based on the same regional monthly sampling as those of all Figs. 3–7, except Figs. 3c and 3d.) It is clear that the number of months sampled exceeds 144 months at mid- to low latitudes and tends to decrease with the latitude. At high latitudes, approximately half of the data (72 months) are usable. To avoid the bias of seasonally averaged values caused by the missing month in a season, we use the monthly means (i.e., 1-month increment) of different parameters to perform the following analysis.
Prior to the regression analysis, all of the variables are interpolated into a 1° × 1° grid to match the CERES grid size. Next, we select the predictor variables for each grid by calculating the temporal correlations between planetary albedo anomalies and different parameters’ anomalies. If the confidence level of the temporal correlation between two variables is less than 90% (i.e., p > 0.1), the variable is excluded. Note that the monthly anomalies of different predictors are already deseasonalized and that their long-term trends are removed and normalized. For details regarding the data-processing workflow, please see Fig. 2.
where c is a constant term; X1, X2, …, Xi are the predictor variables; i is the number of predictor variables; and is the planetary albedo anomaly. We use the stepwise method to remove the insignificant terms from multilinear models based on their statistical significance in the regression process (confidence level > 90%). Furthermore, the relative contribution of each variable to the regional planetary albedo anomaly can be derived from the following formula (Huang and Yi 1991):
where m is the length of the data series, a is the number of independent variables, Aij = bjxij, bj denotes the regression coefficients of each term, xij represents the predictor variables, and j is the number of predictor variables. Finally, we also calculate the coefficient of determination R2 and root-mean-square error (RMSE) between the observed and predicted planetary albedo anomalies in order to evaluate the performance of the regression models.
a. Global characteristics
Figures 3a and 3b show the global distributions of the averaged planetary albedo and cloud albedo forcings, respectively, during 2003 to 2015. Loeb et al. (2007) found that the main source of tropical albedo variability is attributable to cloudiness variations associated with the El Niño–Southern Oscillation (ENSO) phenomenon. However, the strongest El Niño since 1998 began in late 2015. It means that the variation of planetary albedo during the time period between 2003 and 2015 may be independent of the large climate forcing (e.g., El Niño). Generally, the global distribution of PA is in agreement with the findings of Donohoe and Battisti (2011). The mean PA ranges from 0.1 to 0.75 and tends to increase with the latitude. The geographical distribution of the PA is associated with the land–sea, cloud, and snow-/ice-cover distributions (Wielicki et al. 2005; Bender et al. 2006). For example, the PA is relatively large (>0.5) at high latitudes (poleward of 60°) in two hemispheres owing to the presence of wide ice sheets, which enhance the surface albedo. In particular, the PA exceeds 0.65 over Antarctica and Greenland, and the highest planetary albedo (>0.7) is predominantly located over northwestern Antarctica. The variation of PA in the Southern Hemisphere (SH) is much less than that in the Northern Hemisphere (NH) due to the differences in the land–sea distribution. The PA at low latitudes (between 30°N and 30°S) ranges from approximately 0.1 to 0.4. Compared with land, the values of the planetary albedo over the ocean are smaller and range from approximately 0.1 to 0.25. Lower PA values (<0.2) are predominantly located over subtropical oceans, where low-altitude cumulus clouds are frequently observed. Based on Fig. 3b, we can see that the cloud albedo forcing value over the ocean is larger than those values over the high-PA regions (e.g., ice-covered or semiarid/arid zones) since there is greater contrast between ocean and cloud reflectance than land and cloud reflectance because oceans have the darkest surface reflectance. Here, we also provide the global distributions of the PA and cloud albedo forcing differences between Aqua and Terra (Figs. 3c,d), respectively. Note that the Terra and Aqua cloud albedo forcing differences are calculated based on the CERES SSF product. Figure 3c shows that those obvious negative differences in PA were mainly located in typical marine stratocumulus regions (e.g., the Californian, Canarian, Namibian) owing to the decrease in the amounts stratocumulus clouds observed from dawn to afternoon (Garreaud and Muñoz 2004). Over the Tibetan Plateau, South Africa, the northern part of South America, and the western part of North America, the PA during the afternoon is noticeably larger than those values measured during the morning. This difference may be linked to the land afternoon convection (Yang et al. 2004). Indeed, Fig. S2 indicates that the CF and IWP during the afternoon are higher than those results found during the morning. By checking the cloud albedo forcing, we find that the global pattern of cloud albedo forcing difference is very similar to the result of the PA difference, except over the higher Southern Ocean latitudes. This means that the small differences of PA over the subtropical oceans are mainly caused by the weak contrast of cloud albedo forcing, which is because the differences of the macrophysical cloud properties (e.g., CF, LWP, and IWP) between Aqua and Terra are small between two instantaneous observation times (see Fig. S2). Over the Southern Ocean, the cloud albedo forcing during the afternoon is obviously larger than those values observed during morning, but PA exhibits a weak difference. The obvious positive difference of cloud albedo forcings over this region is mainly due to the lower planetary albedo under clear-sky condition during the afternoon. In addition, the differences of PA are small in the mid–high latitudes because the amount of convection is limited owing to the small amount of solar incident radiation. In the following study, it is noted that the Aqua and Terra observation datasets show similar results with regard to the regression analysis and contribution calculations; thus, only the statistical results from Aqua are presented in the current study.
The question then arises as to which one of the factors dominates the long-term variations in the PA over different regions. To address this issue, we first analyze the temporal correlations between the anomalies of PA and the surface albedo (and cloud albedo forcing) anomalies at a global scale (see Fig. 4). Note that the anomalies of each grid are already deseasonalized and detrended.
Statistical results indicate that the temporal correlations between the PA and surface albedo anomalies are almost positive all over the world, except over the central and western Pacific, which may be caused by the small magnitude of the surface albedo oscillations (see Fig. S10). The correlation coefficients increase with increases in latitude. At high latitudes, especially over the Arctic and Southern Ocean, these correlations are equal to or greater than 0.6, indicating that long-term changes in the PA are remarkably consistent with the changes in surface albedo. For the cloud albedo forcing anomalies (Fig. 4b), it is clear that the obvious positive temporal correlations are mainly located at those regions between 60°S and 60°N, especially over the oceans, where these correlations are almost above 0.9. Meanwhile, negative temporal correlations between PA and cloud albedo forcing anomalies can be found in regions where the surface albedo anomalies are significantly positively correlated with the planetary albedo anomalies (e.g., the higher Southern Ocean latitudes). Over these regions, cloud fraction anomalies are out of phase with snow and ice albedo anomalies (see Fig. S11), indicating that surface albedo feedback diminishes (Qu and Hall 2005). By using the CERES SSF product, Kato et al. (2006) found that the effect of decreased Arctic sea ice on albedo might be compensated for by an increase in cloud cover due to enhanced evaporation from the sea surface. Thus, they concluded that any ice-albedo feedback could be dampened due to an increased cloud cover. Hence, whether climate simulations capture this feature is critical for high-latitude regions.
b. Regression coefficients and contribution calculation
Figures 5a–g show the global distributions of the regression coefficient of each variable based on stepwise multiple linear regression models acquired from Eq. (6). Those regions without values indicate that the corresponding variable is not considered a predictor variable in the regression model (i.e., it fails to pass the t test). In addition, if both of the surface parameters (i.e., the NDVI and I/SPC) fail to pass the confidence test, the surface albedo anomaly is considered as a predictor variable when its confidence level exceeds 90%. From Figs. 5a–g, we can see that the coefficients from the regression model vary among the different variables and regions. The cloud properties (e.g., CF, LWP, and IWP) show positive coefficients, which means that larger CF, LWP, and IWP will enhance PA. High coefficients (e.g., >0.6) for CF are mainly located at tropical and subtropical regions, while the maximum values of the regression coefficients for LWP and IWP are concentrated in the Southern Ocean and Pacific warm pool, respectively. For SA, NDVI, and I/SPC, the coefficients only focus on some special regions. For example, a negative coefficient for NDVI can be found for a landmass at high latitude in the NH, where the increased vegetation trends reduce the surface albedo (Bonan 2008; Li et al. 2018). The positive correlation between I/SPC anomalies and PA anomalies is extremely apparent over the higher Southern Ocean latitudes and the Arctic. This means that wide I/SPC covers enhance the PA. However, for Antarctica the variation of I/SPC over this region is small enough that the correlation between the I/SPC and PA anomalies is insignificant. However, Fig. 5g indicates that the surface albedo anomalies have positive correlations with the anomalies of PA over Antarctica. Recent studies have indicated that Antarctica experienced a positive phase of the Antarctic Oscillation index (AAO) trend during the period of 1983 to 2009, which means that the whole of the South Pole exhibits a cooling trend and an increasing snowfall and ice mass when there is a strong polar vortex (Seo et al. 2016). Because of the fact that the albedo of snow varies with the snow condition during and after snowfall (Dang et al. 2016), we therefore speculate that the surface albedo anomalies over Antarctica are possibly caused by the changes in snow density.
Figure 5f shows that reduced T2 is associated with increased PA over the land of NH. It may be linked with the effect of aerosol on the radiation. Higher aerosol loading may decrease the transparency of the atmosphere to solar radiation and enhance the reflected solar radiation under clear-sky conditions. Finally, it shows a weak negative regression coefficient over the land of the NH. In addition, we also provide the confidence interval range (α = 0.1) for each of the regression coefficients (see Figs. S3–S9). Furthermore, the R2 and the RMSE values between the observed and regressed PA anomalies are given in Figs. 5h and 5i. Again, note that our analysis is based on deseasonalized and normalized anomalies by removing the long-term trends in every grid. Statistical results indicate that the regression model can capture changes in the observed long-term PA anomalies in most regions. Almost all of the R2 values are greater than 0.4, particularly at the low latitudes, where the R2 is greater than 0.7.
Finally, we quantify the contribution rates of different variable to the PA anomalies based on Eq. (7). Figure 6 shows the global distributions of the contributions for each predictor variable. Furthermore, the global distributions of the dominant factor and its contribution rate are shown in Fig. 7. Figures 6a and 7a clearly show that the CF anomalies dominate the PA anomalies, especially over the oceans at middle and low latitudes, where its contributions even exceed 70%. This means that climate models need to reasonably simulate the total cloud fraction for improving the prediction of PA in a warmer climate. However, many studies have shown that the reliable simulation of the total cloud fraction in the climate models is strongly dependent on the representation of cloud overlap properties, which is obviously dependent on the observation techniques used (Huang et al. 2005; Li et al. 2011, 2015). At high latitudes, however, it is clear that both the anomalies in CF and I/SPC (or SA) are important to the change of PA. Compared with the CF, the contributions from LWP and IWP anomalies are both secondary (see Figs. 6b,c and Figs. 7a,b). For most of the regions at middle and low latitudes, the contributions of the LWP and IWP anomalies are less than 30%, except for those in the stratus/stratocumulus, lower Southern Ocean latitudes, and the western Pacific warm pool, where the contributions from the LWP and IWP may reach 40%. Previous studies have verified that the obvious differences in the albedo feedbacks over the Southern Ocean found using different models are mainly caused by the inconsistent poleward redistribution of the cloud liquid water content (Tsushima et al. 2006; Hu et al. 2010) and, therefore, result in climate model errors in the predicted TOA fluxes over the Southern Ocean are the largest (Trenberth and Fasullo 2010). Based on CloudSat and CALIPSO cloud observations, Mason et al. (2014) indicated that mid-topped clouds are responsible for biasing the absorbed shortwave radiation in climate models. Indeed, our results indicate that the effect of LWP anomalies on the PA anomalies is nonnegligible over the lower Southern Ocean latitudes. This means that models need to discriminate the cloud phase over this region more reliably, especially with regard to supercooled water clouds (Hu et al. 2010).
The contributions from T2 anomalies are more obvious over land than over ocean. As stated in the above analysis, the contributions of T2 may be related to water vapor and aerosol loading. For the NDVI anomalies (Fig. 6d), we find that their contributions to the PA anomalies over most of the land regions are below 20%. Park et al. (2015) suggested that vegetation greening is observed in response to regional warming in the NH, indicating that vegetation feedback processes (e.g., surface albedo) would be strengthened. A recent finding by Li et al. (2018) verifies that the greening trend at high latitudes made a greater contribution to the decline in surface albedo in the NH (i.e., Siberia), which has experienced a pronounced decrease in surface albedo due to the increase in the extent of evergreen conifer forest cover resulting from the warming over the past several decades (Kharuk et al. 2005, 2007; He et al. 2017). However, our results indicate that the contribution (<10%) from the changes in vegetation cover plays a limited role in the long-term variations of the PA over landmass. Even so, the response of the PA to the vegetation changes over the arid and semiarid regions (e.g., the western United States) is still noteworthy. On the one hand, Forzieri et al. (2017) showed that an increasing vegetation trend contributed to a reduction in surface albedo and to evaporation-driven cooling in arid regions. An enhanced evaporation could additionally decrease the soil moisture and, therefore, lead to desertification (Huang et al. 2017). On the other hand, Huang et al. (2017) noted that dry soils and smaller leaf areas contributed to an increase in the surface albedo and a reduction in transpiration, thereby further intensifying drought.
The relative contributions of I/SPC and SA are shown in Figs. 6e and 6g, respectively. Statistical results clearly show that the contributions from I/SPC and SA anomalies are dominant over the higher Southern Ocean latitudes, Arctic, and Antarctic, respectively. Over the higher Southern Ocean latitudes, the contribution from I/SPC exceeds 40%, whereas the contributions are smaller than 30% over other land in the mid–high latitudes. Pistone et al. (2014) concluded that cloud albedo feedback is playing an insignificant role in the observed Arctic warming. However, our statistical results show that the both CF and I/SPC anomalies are very important to the long-term change of PA over the Arctic. This result is consistent with that of the study of Qu and Hall (2005), who found that the interannual variability in the PA over the cryosphere is dominated by surface albedo changes, but the atmospheric damping effect due to cloud fluctuations can significantly attenuate these surface changes by as much as 90%.
4. Conclusions and discussion
The planetary albedo, which is controlled in a complex way by various atmospheric and surface properties, plays a key role in regulating the global and regional energy budgets. However, our incomplete understanding of the related physical processes (e.g., cloud process) means that the reliable simulation and reproduction of PA in the climate models remain challenging (Bender et al. 2006). To improve the simulation of PA and predict how the planetary albedo will respond to climate change, one of the remaining issues is determining which factor dominates the temporal variability of the planetary albedo based on observations, especially at a regional scale. As a result, this study utilizes 13 years (2003–15) of data from multiple satellite datasets to evaluate the contributions of atmospheric and surface parameters to the long-term variations of the gridded planetary albedo.
By performing a temporal correlation analysis between the planetary albedo anomalies and surface albedo (and cloud albedo forcing) anomalies, we find that the variations in planetary albedo show an obvious positive correlation with the cloud albedo forcing anomalies over the regions between 60°S and 60°N, whereas the negative correlations are mainly located over the Arctic and high-latitude Southern Ocean. Meanwhile, surface albedo anomalies also exhibit an apparent positive correlation with planetary albedo anomalies at high latitudes. Based on a stepwise multiple linear regression analysis, our statistical results indicate that the variations in planetary albedo are closely related to the variations in cloud properties (e.g., cloud fraction, ice water path, and liquid water path), ice/snow percent coverage, and NDVI; however, their temporal relationships vary among the different regions. Generally, the regression model is able to capture the observed planetary albedo anomalies for most regions. The contribution calculation shows that the variations of cloud properties, especially in CF, dominate the long-term variations in the PA over the most of global areas. This conclusion is consistent with the results of previous studies (e.g., Seinfeld et al. 2016). Aside from those of CF, the contributions of I/SPC and SA anomalies are dominant over the high latitudes of the Southern Ocean and part of the Antarctic regions, respectively. In addition, the effects of LWP and IWP are nonnegligible, and their contributions are the second-most important factor globally, except for in the Arctic and Antarctic regions.
At present, the model simulation of PA still suffers from some uncertainties. For example, by using the multimodel dataset from phase 3 of the Coupled Model Intercomparison Project (CMIP3) and two satellite datasets [the Earth Radiation Budget Experiment (ERBE) and CERES], Bender et al. (2006) found that seasonal variations in the PA are captured to some extent by models that span the storm tracks, whereas the albedo is not well reproduced by models spanning the entire oceans. They speculated that this bias may be due to the poor simulation of the solar angle by the models. Meanwhile, they also noted that the poor simulation of cloud in models is responsible for the failure to reproduce the seasonal PA variations over the subtropical arid regions. In addition, Stephens et al. (2015) showed that the reflected energy simulated by the CMIP5 models failed to reproduce the observed hemispheric symmetry, and the modeled global albedo values were systematically higher (by almost 10%) than those from CERES observed during the boreal summer months. This bias has persisted from CMIP3 (Bender et al. 2006). These studies clearly showed that the uncertainties in the simulation of PA are closely related to the cloud and surface parameters, especially CF. In future work, the inclusion of the cloud overlap assumption, which considers the effects of dynamic factors, in the models should improve the simulation of the total cloud fraction and reduce the bias of PA caused by CF (Di Giuseppe and Tompkins 2015). However, over the semiarid (or high pollution) regions, vegetation cover and aerosol loading need to be further considered for improving the parameterization of the planetary albedo and determining climate feedbacks over these regions (Huang et al. 2014).
This research was jointly supported by the National Key R&D Program of China (2016YFC0401003), the Foundation for Innovative Research Groups of the National Science Foundation of China (Grant 41521004), National Science Foundation of China (Grant 41575015), and the Fundamental Research Funds for the Central Universities (Lzujbky-2017-64). We thank the CERES, MODIS, and CM SAF science teams for providing excellent and accessible data products that made this study possible. We are also grateful for the constructive comments from the editor and two anonymous reviewers to improve the quality of this paper.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-17-0848.s1.