Interannual Variations of TOA Albedo over the Arctic, Antarctic and Tibetan Plateau in 2000–2019

: Recent changes in Earth’s climate system have signiﬁcantly a ﬀ ected the radiation budget and its year-to-year variations at top of the atmosphere (TOA). Observing high-latitude TOA ﬂuxes is still challenging from space, because spatial inhomogeneity of surface / atmospheric radiative processes and spectral variability can reﬂect sunlight very di ﬀ erently. In this study we analyze the 20-year TOA ﬂux and albedo data from CERES and MISR over the Arctic, the Antarctic, and Tibetan Plateau (TP), and found overall great consistency in the TOA albedo trend and interannual variations. The observations reveal a lagged correlation between the Arctic and subarctic albedo ﬂuctuations. The observed year-to-year variations are further used to evaluate the reanalysis data, which exhibit substantial shortcomings in representing the polar TOA ﬂux variability. The observed Arctic ﬂux variations are highly correlated with cloud fraction (CF), except in the regions where CF > 90% or where the surface is covered by ice. An empirical orthogonal function (EOF) analysis shows that the ﬁrst ﬁve EOFs can account for ~50% of the Arctic TOA variance, whereas the correlation with climate indices suggests that Sea Ice Extent (SIE), North Atlantic Oscillation (NAO) and 55 ◦ N–65 ◦ N cloudiness are the most inﬂuential processes in driving the TOA ﬂux variabilities. Conceptualization, validation, Y.-K.L.; K.-M.K.; writing—original draft preparation, D.L.W.; writing—review and D.L.W., J.N.L., K.-M.K., Y.-K.L.; visualization, D.L.W.,


Introduction
The sea ice extent (SIE) has been decreasing substantially as a result of the recent polar warming and is subject to large year-to-year variability in Earth's climate system [1]. The sea ice loss, glacier retreat, and snow cover have become more sensitive to climate forcing and interannual temperature fluctuations than several decades ago, as sea ice becomes thinner [2] and snowfall becomes more intense and intermittent [3]. These changes have led to a significant impact on Earth's radiation budget. The estimated annual-mean absorbed solar radiation over 75 • N-90 • N would increase by 2.5 W/m 2 for every 10 6 km 2 SIE decrease in September [4]. The decreasing Arctic albedo in the recent years appears to be consistent with the SIE retreat, despite increases in cloudiness over open water [5,6]. Yet, these disturbances are not confined to the polar regions as Earth's climate system is highly coupled, including poleward heat transport [7], remote oceanic influences [8], and potential Arctic impacts on mid-latitude weather and climate [9][10][11]. Because the radiation balance at top of the atmosphere (TOA) is sensitive to the surface/cloud changes induced by the polar warming, the oscillatory response and feedback to these changes are important coupling processes in the Earth's climate system, and they become increasingly detectable in the rapid warming period. The climate community are now facing several critical challenges. On one hand, most of the climate models struggle to reproduce the observed TOA Figure 1 also reveal a lagged correlation between the mid-and high-latitude flux oscillations with the Arctic fluctuations leading the sub-Arctic (e.g., 60°N) ones by 1~2 years. Since 2009, this lagged relationship has become more evident due to the increased amplitude of Arctic fluctuations. It remains debatable, however, whether the changing Arctic or its increased variabilities can significantly impact the mid-latitude climate or weather systems [10], and vice versa [8]. Although TOA fluxes are not the only factors that drive the near-surface variabilities, these precise measurements are indicative of some key climate-system changes that have caused an additional imbalance in the Arctic at a decadal scale.
The CERES EBAF data also contain monthly cloud fraction (CF) on a 1° × 1° longitude-latitude grid, which is derived essentially from MODIS (Moderate Resolution Imaging Spectroradiometer). As an important variable to interpret the observed variations in the TOA SW flux, CF can have a major influence on how much sunlight is reflected back to space. Although clouds, unlike sea ice and snow cover, do not always occur at the same location, their impacts on the SW flux are realized through the frequency of occurrence, which is effectively integrated in the monthly CF data. The distribution of monthly cloud amount, as well as cloud occurrence frequency, is useful diagnosis variables of the TOA SW variability.

TOA Albedo Data
Although the CERES EBAF Ed4.1 data are significantly improved over the earlier versions, producing accurate all-sky TOA radiative fluxes requires an algorithm to take into account correction of geostationary-Earth-orbit (GEO) satellite calibration artifacts, CERES-MODIS narrowband-tobroadband (N2B) radiance conversion, and cloud mask and scene classification [14]. In the polar region, the TOA SW flux measurement may subject to greater uncertainty than the middle and lower latitudes Figure 1 also reveal a lagged correlation between the mid-and high-latitude flux oscillations with the Arctic fluctuations leading the sub-Arctic (e.g., 60 • N) ones by 1~2 years. Since 2009, this lagged relationship has become more evident due to the increased amplitude of Arctic fluctuations. It remains debatable, however, whether the changing Arctic or its increased variabilities can significantly impact the mid-latitude climate or weather systems [10], and vice versa [8]. Although TOA fluxes are not the only factors that drive the near-surface variabilities, these precise measurements are indicative of some key climate-system changes that have caused an additional imbalance in the Arctic at a decadal scale.
The CERES EBAF data also contain monthly cloud fraction (CF) on a 1 • × 1 • longitude-latitude grid, which is derived essentially from MODIS (Moderate Resolution Imaging Spectroradiometer). As an important variable to interpret the observed variations in the TOA SW flux, CF can have a major influence on how much sunlight is reflected back to space. Although clouds, unlike sea ice and snow cover, do not always occur at the same location, their impacts on the SW flux are realized through the frequency of occurrence, which is effectively integrated in the monthly CF data. The distribution of monthly cloud amount, as well as cloud occurrence frequency, is useful diagnosis variables of the TOA SW variability.

TOA Albedo Data
Although the CERES EBAF Ed4.1 data are significantly improved over the earlier versions, producing accurate all-sky TOA radiative fluxes requires an algorithm to take into account correction of geostationary-Earth-orbit (GEO) satellite calibration artifacts, CERES-MODIS narrowband-to-broadband (N2B) radiance conversion, and cloud mask and scene classification [14]. In the polar region, the TOA SW flux measurement may subject to greater uncertainty than the middle and lower latitudes due to scene classification and diurnal sampling correction (without GEO). To evaluate and verify the interannual variations in the CERES TOA SW flux, we compared the result with the observations obtained by Multiangle Imaging SprectroRadiometer (MISR) on the same platform but using a very different technique. The comparison is carried out in terms of the TOA albedo, which can be calculated by dividing the reflected TOA SW flux over the incident solar insolation.
Different from the CERES technique, which relies on empirical angular distribution models (ADMs) to convert the single-angle broadband radiance to the TOA flux [15,16], the MISR technique can measure instantaneous ADM at four narrow spectral bands with nine-angular radiance measurements in the along-track direction [17]. To obtain the TOA albedo, however, the MISR technique needs to use modeled N2B functions, which are scene-dependent, to convert its narrow-band radiances to the broadband SW radiances. Thus, the MISR algorithm also requires scene classification. Because the scene classification (into thousands of types) algorithms may differ depending on observing techniques, some differences in the TOA albedo are expected. For example, it is challenging to classify clouds over sea ice. Hence, the two techniques may yield different SW albedos, as the scene type classification can be different. Mountainous terrains with substantial 3D inhomogeneity can be particularly challenging for the techniques that rely on single-view radiance measurements for the SW flux determination.
In addition, CERES and MISR have quite different pixel resolutions (10 km vs. 275 m) and sampling swaths (3000 km vs. 400 km). The larger swath allows CERES to cover both North and South Pole during the summer months while MISR has a gap at latitudes from 84 • poleward. Observations from a narrow swath also have a larger sampling error due to fast processes like clouds. Thus, each technique has its own advantages and limitations in handling Arctic atmospheric/surface inhomogeneity and variability, and differences between the two TOA albedo products are expected [e.g., Zhan et al., 2018].
Here, an objective of MISR-CERES albedo intercomparisons is to identify robust and consistent features over the Arctic, Antarctic and Tibetan Plateau, given their sampling and conversion method differences. Because of differences in scene classification, viewing geometry, swath sampling and narrow-vs-broad radiances, verification of independent CERES and MISR measurements provides fidelity of the observed TOA SW features and their variabilities.
This study also compares the observed TOA albedo variations with those derived from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) and the fifth generation of European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyses (ERA5). MERRA-2 derives the TOA radiations fluxes by assimilating available satellite radiance and conventional observations. The MERRA-2 background atmospheric model is the Goddard Earth Observing System model (GEOS) [18]. MERRA-2 is produced on a 5/8 • longitude × 1/2 • latitude grid in 72 vertical levels. For MERRA-2, the surface representation has been substantially revised from MERRA to include snow hydrology processes, a prognostic albedo, and energy conductivity through snow and ice layers at a high vertical resolution [19]. Other improvements in MERRA-2 include modifications to the representation of turbulence and moisture-related processes. In this study we use the "tavgM_2d_rad_Nx: 2d, monthly mean, 1-hourly time averaged, assimilation, single-level diagnostics V5.12.4 [20], to calculate the TOA albedo. Differences between the TOA incoming solar flux values used by EBAF and the MERRA-2 reanalyses can cause a bias in the TOA SW fluxes. While a~1 W m −2 SW downward flux adjustment was made for MERRA-2 compared to that of MERRA, there still exists a difference of 6.1 W m −2 in TOA SW upward flux from that of EBAF [21]. In this study we focus primarily on comparisons of the interannual variations from these data sets.
In ERA5, the radiative flux variables have a spatial resolution of~31 km (0.28125 • × 0.28125 • ) and an hourly temporal resolution [22]. ERA5 uses the Integrated Forecasting System (IFS) with a 4-D variational analysis (4DVAR) assimilation system. We used the ERA5 Top net Solar Radiation (TSR) and Top net Thermal Radiation (TTR) [23].

Index Data for Arctic Analyses
Arctic sea ice extent (SIE) is one of the key variables that drives the summertime TOA albedo. We use the monthly SIE data archived at the National Snow and Ice Data Center (NSIDC) as an index to characterize the year-to-year variation of sea ice change. As discussed later in the paper, the TOA fluxes will be regressed to the SIE index, to understand its impacts not only on the Arctic but also on the surrounding sub-Arctic and mid-latitude regions.
For this reason, we employ another key regional index, the North Atlantic Oscillation (NAO), which is the difference of atmospheric sea level pressure between the Azores High and the Icelandic Low. Although NAO is derived from regional weather patterns, it has a profound impact on remote area including the Arctic basin. We also employ the Arctic Oscillation (AO) index as a representative of the Northern Annular Mode that describes the annular pattern of atmospheric circulation anomalies across the mid-latitudes and the Arctic. In this study we obtained the monthly NAO index from NOAA National Centers for Environmental Prediction (NCEP).
Cloud and snow cover are important factors that can affect the TOA albedo. To characterize northern hemisphere (NH) cloudiness, we define a cloudiness index (CI) as the gridbox mean CF averaged from latitudes of 55 • N-65 • N in the CERES EBAF data (see Section 4 for more details). This averaging is weighed by the gridbox area. For snow cover (SC), we use the total SC at latitudes of 40 • N poleward, which is derived from MODIS observations [24], as an index to characterize the land area covered by snow. As a secondary effect, ice/snow albedos may differ slightly, depending on their compositions (e.g., aerosol deposition, new/old ice) and mixtures (e.g., ice with leads).
To evaluate impacts of the El Niño-Southern Oscillation (ENSO) on the Arctic TOA fluxes, we employ the NOAA monthly Oceanic Nino Index (ONI) index, which is defined as the average sea surface temperature (SST) anomalies in the Niño 3.4 region (5 • S-5 • N; 170 • W-120 • W). Interannual variations at mid-to-high latitudes can be linked remotely to the tropical ocean change through teleconnection processes. Regression to the ONI index will help to quantify such a teleconnection impact. Similarly, a spring-summer-fall pattern in the surface temperatures over East Pacific-North Pacific (EPNP) provides a key index for intense storm activity over the mid-latitudes of the North Pacific [25].

Results
The TOA fluxes over the Arctic, Antarctic, and Tibetan Plateau (TP) are experiencing a significant change in recent years due to the vulnerability of their underlying sea ice and snow cover to the regional warming. During 2000-2019, the period with CERES and MISR observations, the global surface temperature has increased by~0.5 • C and the Arctic by~1.5 • C [26]. These warming trends as well as interannual variability of the warming play a fundamental role in the observed TOA SW flux.

Arctic
The observed mean JJA albedo distribution over the Arctic is shown Figure 2, along with the time series of zonal mean albedo perturbations as a function of latitude from MISR, CERES, ERA5 and MERRA2. Both MISR and CERES observations reveal a similar TOA albedo distribution with the brightest over Greenland and Arctic sea ice. The CERES albedo over the Pacific and Atlantic is mostly from clouds, which is slightly brighter than MISR measurements. On the other hand, MISR albedo perhaps has slightly more variability over landmasses. A decreasing TOA albedo trend is clearly evident in both CERES and MISR JJA data at latitudes greater than 60 • N where the Arctic change is dominated by the summer sea ice loss. The lowest Arctic albedo for JJA occurred in 2011 during the two-decade observation period. In addition, there is a strong interannual oscillation on the top of the trend, showing low albedos in 2008, 2011 and 2015 at these latitudes. It is noteworthy that these low-albedo years are one year ahead of the large mass loss anomalies over Greenland Ice Sheet (GrIS) in 2009, 2012 and 2016 [27]. At the 60 • N-80 • N latitude bin, GrIS is the dominant source of the TOA albedo. Despite their very different remote sensing techniques, CERES and MISR observations reveal an overall consistent pattern of the year-to-year TOA albedo variations at the latitudes (50 • N-82 • N). This gross consistency provides high fidelity on the TOA SW oscillations as observed by these sensors from 2000 to 2019, as well as the TOA albedo distribution at high latitudes.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 22 an overall consistent pattern of the year-to-year TOA albedo variations at the latitudes (50°N-82°N). This gross consistency provides high fidelity on the TOA SW oscillations as observed by these sensors from 2000 to 2019, as well as the TOA albedo distribution at high latitudes. Compared to the high-latitude (70°N-80°N) oscillations, the sub-Arctic (60°N-70°N) variations appear to have a lag of 1-2 years. The lagged correlation is more evident after 2011 when the interannual oscillations intensify. The oscillations at latitudes > 70°N show a period of ~4 years during 2006-2019, whereas the lagged oscillations in the sub-Arctic seem to preserve the periodicity during the period after 2011. As discussed above, this lagged albedo variation in Figure 2 is largely driven by the variabilities of GrIS at 65°N and sea ice at 80°N.
The interannual oscillations of TOA albedo in 2000-2019 are captured, to some extent, in the MERRA2 and ERA5 reanalysis data. Similar perturbations in the MERRA2 and ERA5 albedos are found at latitudes < 70°N. However, for latitudes > 70°N the reanalysis data differ from each other substantially. Their differences are much larger than those found between CERES and MISR, indicating problems in the reanalysis at high latitudes. For example, at latitudes > 80°N, the oscillations from MERRA2 and ERA5 even disagree in sign, with MERRA2 more in agreement with the CERES observation. The CERES observations, as well as in the MISR observations up to 82°N, show that the interannual oscillations have a coherent extension to the pole, whereas the ERA5 and MERRA2 reanalysis data exhibit a break at 80°N. Since the albedo at these latitudes come mostly from the sea-ice-covered Arctic Ocean, it raises a question whether some important processes (e.g., melt pond on sea ice) are missing in the underlying model of the reanalysis. The melt pond fraction is one of the key factors needed to determine the heat budget of ice-covered area [28], but it is difficult to measure from space directly. The interannual oscillations of TOA albedo in 2000-2019 are captured, to some extent, in the MERRA2 and ERA5 reanalysis data. Similar perturbations in the MERRA2 and ERA5 albedos are found at latitudes < 70 • N. However, for latitudes > 70 • N the reanalysis data differ from each other substantially. Their differences are much larger than those found between CERES and MISR, indicating problems in the reanalysis at high latitudes. For example, at latitudes > 80 • N, the oscillations from MERRA2 and ERA5 even disagree in sign, with MERRA2 more in agreement with the CERES observation. The CERES observations, as well as in the MISR observations up to 82 • N, show that the interannual oscillations have a coherent extension to the pole, whereas the ERA5 and MERRA2 reanalysis data exhibit a break at 80 • N. Since the albedo at these latitudes come mostly from the sea-ice-covered Arctic Ocean, it raises a question whether some important processes (e.g., melt Remote Sens. 2020, 12, 1460 7 of 21 pond on sea ice) are missing in the underlying model of the reanalysis. The melt pond fraction is one of the key factors needed to determine the heat budget of ice-covered area [28], but it is difficult to measure from space directly.

Antarctic
In the Antarctic the CERES and MISR TOA DJF albedo data show similar interannual variations ( Figure 3). Both data sets reveal an albedo increase in 2013-2015 at 65 • S, followed by a sharp reversal in 2016. While the 2013-2015 albedo increase is confined near 65 • S, the reversal occurs over a wide range of latitudes between 60 • -80 • S. The 65 • S latitude belt is where most of the Antarctic sea ice variability is seen. The dramatic decrease in Antarctic sea ice in 2016 has been intensively studied and attributed potentially to a number of factors including Southern Annular Mode (SAO) and Interdecadal Pacific Oscillation (IPO) [29], as well as Pacific South American (PSA) teleconnection [30].

Antarctic
In the Antarctic the CERES and MISR TOA DJF albedo data show similar interannual variations ( Figure 3). Both data sets reveal an albedo increase in 2013-2015 at 65°S, followed by a sharp reversal in 2016. While the 2013-2015 albedo increase is confined near 65°S, the reversal occurs over a wide range of latitudes between 60°-80°S. The 65°S latitude belt is where most of the Antarctic sea ice variability is seen. The dramatic decrease in Antarctic sea ice in 2016 has been intensively studied and attributed potentially to a number of factors including Southern Annular Mode (SAO) and Interdecadal Pacific Oscillation (IPO) [29], as well as Pacific South American (PSA) teleconnection [30]. The TOA albedo oscillations induced by the Antarctic sea ice are perhaps represented better by ERA5 than by MERRA2. The 2016 albedo reversal in ERA5 is more consistent with the CERES and MISR observations, but it appears to be delayed by a year. On the other hand, MERRA2 lacks variability in the TOA albedo during 2000-2019. In addition, MERRA2 shows a much larger trend in the earlier period  at 60°S-70°S latitudes, compared to the ERA5 data, indicating a great challenge in these reanalyses to represent the energy balance in the Antarctic.

Tibetan Plateau (TP) and High Mountain Asia (HMA)
Unlike the Arctic and Antarctic, the TP and HMA (e.g., Afghanistan, Pakistan, and Tajikistan) have highly variable terrains where the 3D effects can induce a larger uncertainty than the flat surface cases. Thus, the TP and HMA regions provide a good testbed to verify the TOA albedo observed by CERES and MISR. As shown in Figure 4, the CERES and MISR albedo measurements reveal a similar distribution and interannual variations. Although the CERES albedo values may be slightly higher than MISR, the distribution patterns over the TP and HMA regions are quite consistent between the two data sets. The TOA albedo oscillations induced by the Antarctic sea ice are perhaps represented better by ERA5 than by MERRA2. The 2016 albedo reversal in ERA5 is more consistent with the CERES and MISR observations, but it appears to be delayed by a year. On the other hand, MERRA2 lacks variability in the TOA albedo during 2000-2019. In addition, MERRA2 shows a much larger trend in the earlier period  at 60 • S-70 • S latitudes, compared to the ERA5 data, indicating a great challenge in these reanalyses to represent the energy balance in the Antarctic.

Tibetan Plateau (TP) and High Mountain Asia (HMA)
Unlike the Arctic and Antarctic, the TP and HMA (e.g., Afghanistan, Pakistan, and Tajikistan) have highly variable terrains where the 3D effects can induce a larger uncertainty than the flat surface cases. Thus, the TP and HMA regions provide a good testbed to verify the TOA albedo observed by Remote Sens. 2020, 12, 1460 8 of 21 CERES and MISR. As shown in Figure 4, the CERES and MISR albedo measurements reveal a similar distribution and interannual variations. Although the CERES albedo values may be slightly higher than MISR, the distribution patterns over the TP and HMA regions are quite consistent between the two data sets.
as these observations have to overcome 3D effects of the complex terrain with very different techniques. December is the month when the largest albedo decreasing trend as well as year-to-year variations are observed by MISR and CERES during 2000-2019. The December oscillations reveal a low albedo in 2005, 2010, and 2016, with a period of ~5 years. The 5-year fluctuation is mixed with a quasi-biennial oscillation on the top of a decreasing trend, driving most of the December variabilities. HMA has the largest glacier and perennial snow cover concentration outside the polar regions. The TOA albedo over the TP and HMA is largely determined by the glacier and snow cover in the region. The MODIS snow cover variation (not shown) is found to be highly correlate with the TOA albedo oscillation seen in Figure 4. Recent studies also show that the snowline altitude at the end of melt season over HMA has been rising during 2001-2016 [31], suggesting a reduction in regional snow cover and surface albedo.

Trends in 2000-2019
To characterize regional variations of the TOA albedo in further detail, we analyze the observed trends in 2000-2019 on a monthly basis for the Arctic, Antarctic and TP/HMA regions. By breaking down the trend distribution by month, it helps to identify the underlying processes that may evolve with time, in particular those from the surface processes that can vary substantially from month to month and region to region.
As shown in Figure 5, decreasing trends of the TOA albedo dominate everywhere in the Arctic in every month. During June-October the obvious decreasing trend is near the sea ice margin zone where sea ice retreat and early start of the melt season have contributed to the large trend seen in 2000-2019. This trend is narrowly following the margin zone in August-September but widen in October. In May-June a significant trend can also be found in the northern Canada, Victoria and Queen Elizabeth Islands, Beaufort Sea, Laptev Sea, and the eastern Europe, where early snow melt The CERES-MISR consistency is perhaps more striking in the observed interannual variations, as these observations have to overcome 3D effects of the complex terrain with very different techniques. December is the month when the largest albedo decreasing trend as well as year-to-year variations are observed by MISR and CERES during 2000-2019. The December oscillations reveal a low albedo in 2005, 2010, and 2016, with a period of~5 years. The 5-year fluctuation is mixed with a quasi-biennial oscillation on the top of a decreasing trend, driving most of the December variabilities. HMA has the largest glacier and perennial snow cover concentration outside the polar regions. The TOA albedo over the TP and HMA is largely determined by the glacier and snow cover in the region. The MODIS snow cover variation (not shown) is found to be highly correlate with the TOA albedo oscillation seen in Figure 4. Recent studies also show that the snowline altitude at the end of melt season over HMA has been rising during 2001-2016 [31], suggesting a reduction in regional snow cover and surface albedo.

Trends in 2000-2019
To characterize regional variations of the TOA albedo in further detail, we analyze the observed trends in 2000-2019 on a monthly basis for the Arctic, Antarctic and TP/HMA regions. By breaking down the trend distribution by month, it helps to identify the underlying processes that may evolve with time, in particular those from the surface processes that can vary substantially from month to month and region to region. As shown in Figure 5, decreasing trends of the TOA albedo dominate everywhere in the Arctic in every month. During June-October the obvious decreasing trend is near the sea ice margin zone where sea ice retreat and early start of the melt season have contributed to the large trend seen in 2000-2019. This trend is narrowly following the margin zone in August-September but widen in October. In May-June a significant trend can also be found in the northern Canada, Victoria and Queen Elizabeth Islands, Beaufort Sea, Laptev Sea, and the eastern Europe, where early snow melt may play a role in reducing the surface albedo. In December a decreasing trend in Europe, particularly over the Alps, reflects a lack of snowfall or snow cover in the recent years.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22 may play a role in reducing the surface albedo. In December a decreasing trend in Europe, particularly over the Alps, reflects a lack of snowfall or snow cover in the recent years. In Figure 6, the Antarctic TOA albedo is dominated by a decreasing trend in all months, which is largely driven by the sea ice trend [32]. During the summer months (e.g., February) the sea ice extent retreats near the coastal line (~67°S), whereas in the winter months (e.g., September) it expands to the latitudes of ~60°S. The decreasing trend over West Antarctica is quite persistent during November-April, which is consistent with the report from CryoSat-2 that the West Antarctic Ice Sheet (WAIS) is disappearing [33]. Over East Antarctica, the decreasing trends appear to be associated with the sea ice margin zones, moving from the Southern and Indian Ocean in August-November to Weddell Sea in December-January. There is a significant decreasing trend over the Atlantic Ocean in September and December, which is likely due to cloudiness changes. Compared to the decreases in sea-ice-induced TOA albedo, the cloud effects are more scattered and diffused in a wider area. Finally, a significant trend is observed over the Chilean Patagonia during the winter-spring months (June-October), which seems to be related to the decreases in snow cover from high-altitude warming [34]. In Figure 6, the Antarctic TOA albedo is dominated by a decreasing trend in all months, which is largely driven by the sea ice trend [32]. During the summer months (e.g., February) the sea ice extent retreats near the coastal line (~67 • S), whereas in the winter months (e.g., September) it expands to the latitudes of~60 • S. The decreasing trend over West Antarctica is quite persistent during November-April, which is consistent with the report from CryoSat-2 that the West Antarctic Ice Sheet (WAIS) is disappearing [33]. Over East Antarctica, the decreasing trends appear to be associated with the sea ice margin zones, moving from the Southern and Indian Ocean in August-November to Weddell Sea in December-January. There is a significant decreasing trend over the Atlantic Ocean in September and December, which is likely due to cloudiness changes. Compared to the decreases in sea-ice-induced TOA albedo, the cloud effects are more scattered and diffused in a wider area. Finally, a significant trend is observed over the Chilean Patagonia during the winter-spring months (June-October), which seems to be related to the decreases in snow cover from high-altitude warming [34]. Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 22 As demonstrated in Figure 4, the consistent TOA albedo observations from CERES and MISR over the TP and HMA allow more quantitative assessment of the regional trend over these mountainous areas. June-August (JJA) is the TP's wet season, followed by a dry season in September-November (SON) [35]. An enhanced warming in TP is strongly coupled with surface albedo changes [36]. HMA gets its most precipitation in December-February (DJF) whereas Taklamakan receives its mostly in March-May (MAM). During its growing season (May-September) decreasing trends in June and August are observed in the central and eastern TP (Figure 7). The decreasing albedo could be caused by a number of factors including the enhanced evapotranspiration cooling on vegetation [37], black carbon (BC) and dust aerosols on glaciers [38], and increases of atmospheric anthropogenic aerosols [39]. There are decreasing trends in November and January over the eastern TP region (e.g., Bhutan, Myanmar and Bangladesh), which is likely related to the observed decreases in snow cover [40] and glacier area [41].
The decreasing trend over HMA in December is highly correlated with the decreasing snow cover [42] and the shrinking glacier area in the region [42,43]. The glacier losses appear to accelerate in the 21st century with a higher rate at lower elevations [44]. The thinning glaciers [45] and reduced snowfall [46] in HMA region may have all contributed to the reduction of surface albedo, which results in the decreasing trend in the TOA albedo. As demonstrated in Figure 4, the consistent TOA albedo observations from CERES and MISR over the TP and HMA allow more quantitative assessment of the regional trend over these mountainous areas. June-August (JJA) is the TP's wet season, followed by a dry season in September-November (SON) [35]. An enhanced warming in TP is strongly coupled with surface albedo changes [36]. HMA gets its most precipitation in December-February (DJF) whereas Taklamakan receives its mostly in March-May (MAM). During its growing season (May-September) decreasing trends in June and August are observed in the central and eastern TP (Figure 7). The decreasing albedo could be caused by a number of factors including the enhanced evapotranspiration cooling on vegetation [37], black carbon (BC) and dust aerosols on glaciers [38], and increases of atmospheric anthropogenic aerosols [39]. There are decreasing trends in November and January over the eastern TP region (e.g., Bhutan, Myanmar and Bangladesh), which is likely related to the observed decreases in snow cover [40] and glacier area [41].
The decreasing trend over HMA in December is highly correlated with the decreasing snow cover [42] and the shrinking glacier area in the region [42,43]. The glacier losses appear to accelerate in the 21st century with a higher rate at lower elevations [44]. The thinning glaciers [45] and reduced snowfall [46] in HMA region may have all contributed to the reduction of surface albedo, which results in the decreasing trend in the TOA albedo. Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 22

Discussions
In this section we attempt to provide a further analysis on the attributions of the interannual variations of Arctic summertime TOA fluxes. We focus on components that contribute most to the total variances seen in 2000-2019 at latitudes of 40°N poleward, such as spatial climatological patterns and weather indices. Two methods are employed to break down the total variability: empirical orthogonal function (EOF) analysis and climate index regression. The EOF analysis will provide the orthogonal variation patterns that yield the most contribution to the total variance, but will not provide any attributions to the physical parameter/process responsible for the derived patterns. On the other hand, the index regression will help sort out the contributions and their correlation to each physical index/process, but may not capture most of the variances in Arctic TOA fluxes. In addition, indices may not be fully independent of each other, which can lead to a large unwanted covariance. Figure 8 shows three of the five leading EOF patterns of the CERES total, LW and SW TOA fluxes. As seen in Table 1, the five EOFs contribute about ~50% of the total variance from latitudes of 40°N poleward, whereas the combination of the first three EOFs gives ~32%. It would take 10 EOFs to capture ~73% and 18 EOFs for ~98% of the total variance. After EOF#5, the variance tends to fall exponentially as a function of , where N is the EOF number and b is an empirical coefficient that can be found in Table 1.

Leading EOF Patterns
These EOFs are derived with an area-weighted analysis, and the first EOF takes 16%, 13% and 15% of the variance respectively for total, LW and SW fluxes. Without weighted by area, the EOFs would be biased to the variabilities at high latitudes, for which the first EOF would capture as more variance as 24%, 19% and 23% respectively. Note that the amplitudes and patterns of the EOFs in the

Discussion
In this section we attempt to provide a further analysis on the attributions of the interannual variations of Arctic summertime TOA fluxes. We focus on components that contribute most to the total variances seen in 2000-2019 at latitudes of 40 • N poleward, such as spatial climatological patterns and weather indices. Two methods are employed to break down the total variability: empirical orthogonal function (EOF) analysis and climate index regression. The EOF analysis will provide the orthogonal variation patterns that yield the most contribution to the total variance, but will not provide any attributions to the physical parameter/process responsible for the derived patterns. On the other hand, the index regression will help sort out the contributions and their correlation to each physical index/process, but may not capture most of the variances in Arctic TOA fluxes. In addition, indices may not be fully independent of each other, which can lead to a large unwanted covariance. Figure 8 shows three of the five leading EOF patterns of the CERES total, LW and SW TOA fluxes. As seen in Table 1, the five EOFs contribute about~50% of the total variance from latitudes of 40 • N poleward, whereas the combination of the first three EOFs gives~32%. It would take 10 EOFs to capture~73% and 18 EOFs for~98% of the total variance. After EOF#5, the variance tends to fall exponentially as a function of e −bN , where N is the EOF number and b is an empirical coefficient that can be found in Table 1.

Leading EOF Patterns
changes across the Northeast Pacific/Alaska and Northern Canada (e.g., Archipelago region). The correlations between the PC time series and the EPNP index are 0.34 (total) and 0.32 (SW). Although the correlation values are not remarkably high, the impact of EPNP appears to be reflected in the EOF#2 pattern, as negative anomalies of total and SW flux along the North America northwest coast and positive anomalies over Northern Canada. Such an oscillating pattern corresponds to the positive phase of EPNP. While there is no significant trend in the time series of EOF#1 and #2 (Figure 9), the PC amplitudes of EOF#3 in 2000-2019 do exhibit a significant trend. The EOF#3 patterns in Figure 9 indicate the vulnerable regions of the future SW flux variation if the albedo continues to decrease. One of these regions is the western Greenland and northern Canada where the SW flux or albedo had  These EOFs are derived with an area-weighted analysis, and the first EOF takes 16%, 13% and 15% of the variance respectively for total, LW and SW fluxes. Without weighted by area, the EOFs would be biased to the variabilities at high latitudes, for which the first EOF would capture as more variance as 24%, 19% and 23% respectively. Note that the amplitudes and patterns of the EOFs in the total and SW fluxes are similar and comparable, indicating that the total TOA flux variability is dominated by the SW. The LW EOF patterns generally have an opposite sign to the SW ones, providing a damping effect to the SW forcing. As discussed later in this section, clouds play a major role in creating such an offset effect.
A significant portion of EOF#1 can be explained by the TOA flux anomaly distributions associated with the variation of NAO phase. The temporal correlation between the time series of EOF#1 PC and summertime NAO index is high: −0.79 (total), −0.51 (LW), and −0.76 (SW). The high NAO-EOF#1 correlation suggests a broad influence of the negative phase of the NAO on the total and SW flux anomalies over the Arctic, which extends to Greenland as seen in EOF#1 (left column of Figure 8). This characteristic feature will be discussed again in Section 4.3. EOF#2 may have a moderate relationship with the EPNP that generally accounts for the large-scale mass/circulation changes across the Northeast Pacific/Alaska and Northern Canada (e.g., Archipelago region). The correlations between the PC time series and the EPNP index are 0.34 (total) and 0.32 (SW). Although the correlation values are not remarkably high, the impact of EPNP appears to be reflected in the EOF#2 pattern, as negative anomalies of total and SW flux along the North America northwest coast and positive anomalies over Northern Canada. Such an oscillating pattern corresponds to the positive phase of EPNP.
While there is no significant trend in the time series of EOF#1 and #2 (Figure 9), the PC amplitudes of EOF#3 in 2000-2019 do exhibit a significant trend. The EOF#3 patterns in Figure 9 indicate the vulnerable regions of the future SW flux variation if the albedo continues to decrease. One of these regions is the western Greenland and northern Canada where the SW flux or albedo had been decreasing in 2000-2019. The other regions, which appear in a wavy pattern, include Europe (decreasing), Finland (increasing), Russia (decreasing), the northeastern China (increasing), and Kara and Beaufort Seas (decreasing). Despite the small (8%) contribution of EOF#3, its trend amplitude must not be neglected since it may reveal severe impacts of the underlying processes (e.g., sea ice loss) on these regions. . Despite the small (8%) contribution of EOF#3, its trend amplitude must not be neglected since it may reveal severe impacts of the underlying processes (e.g., sea ice loss) on these regions.

Impacts of Cloud Cover
Clouds play a significant role in the overall Arctic and subarctic TOA albedo and its variability. While cloud cover may vary from year to year in response and/or feedback to the changes from atmospheric, hydrological and cryospheric processes, their contributions to the TOA SW flux are more than the surface contributors such as those from Arctic sea ice and Greenland ice sheet (GrIS) changes [47].
The authors of Figure 10, characterize the importance of cloud cover by mapping out the correlation between the TOA SW flux and the monthly gridded cloud fraction (CF) from the CERES EBAF data. The climatological mean distributions of the SW flux and CF resemble each other over the regions without sea ice and outside Greenland. The year-to-year variabilities of the CF and SW flux are highly correlated in these regions with the exception of those where the mean CF becomes greater than 90%, such as the Pacific and Atlantic storm tracks. Because these high-CF regions are constantly covered by clouds, it is expected that cloud-saturated regions produce a smaller year-toyear variation. On the other hand, the regions with less-saturated CF, such as landmasses, had experienced a relatively large CF variability that is correlated with the SW flux fluctuations. As expected, the regions covered by sea ice and permanent ice sheets have a low correlation with CF.

Impacts of Cloud Cover
Clouds play a significant role in the overall Arctic and subarctic TOA albedo and its variability. While cloud cover may vary from year to year in response and/or feedback to the changes from atmospheric, hydrological and cryospheric processes, their contributions to the TOA SW flux are more than the surface contributors such as those from Arctic sea ice and Greenland ice sheet (GrIS) changes [47].
The authors of Figure 10, characterize the importance of cloud cover by mapping out the correlation between the TOA SW flux and the monthly gridded cloud fraction (CF) from the CERES EBAF data. The climatological mean distributions of the SW flux and CF resemble each other over the regions without sea ice and outside Greenland. The year-to-year variabilities of the CF and SW flux are highly correlated in these regions with the exception of those where the mean CF becomes greater than 90%, such as the Pacific and Atlantic storm tracks. Because these high-CF regions are constantly covered by clouds, it is expected that cloud-saturated regions produce a smaller year-to-year variation. On the other hand, the regions with less-saturated CF, such as landmasses, had experienced a relatively large CF variability that is correlated with the SW flux fluctuations. As expected, the regions covered by sea ice and permanent ice sheets have a low correlation with CF.
In the regions where CF and SW are highly correlated, the TOA SW and LW fluxes are anti-correlated, suggesting that cloudy-sky produces less outgoing LW radiation. In other words, the cloud top temperature in JJA is generally colder than the clear-sky surface temperature. This anti-correlation between the LW and SW fluxes is found mostly over land and near the coasts, where the cold cloud top and warm surface conditions are often valid. Again, the LW and SW correlation is poor in the regions where CF is greater than 90%.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 22 As seen in Figure 10, this CI index is designed to capture most of the cloud-induced SW variabilities over the subarctic landmasses, Figure 10. CERES mean TOA SW flux and cloud fraction (CF) distributions for JJA (top panels), and correlation maps between year-to-year variations in the LW and SW fluxes and between CF and SW (bottom panels) at latitudes of 40°N poleward. The red contour is the JJA mean sea ice extent at 80%. The regional mean correlation coefficients are indicated in the label of the bottom panels, showing dominated anti-correlation between LW and SW and positive correlation between CF and SW.

Impacts of SIE, SC, NAO, ONI, CI, EPNP, and AO
Since the EOF analysis cannot provide clear attribution to what physical processes are responsible for the interannual oscillations seen in the TOA flux data, we regress CERES TOA fluxes to several climate indices, namely, SIE, SC, NAO, ONI, CI, EPNP, and AO, to evaluate their connection to the Arctic flux variations. SIE, SC, an CI can directly affect the TOA albedo because of relatively high albedos from sea ice, snow and clouds. Other atmospheric (e.g., water vapor, aerosol) and surface (e.g., wetlands, vegetation, melt ponds) variables may directly or indirectly alter the TOA albedo. Collectively, these albedo contributors are subject to atmospheric, hydrological, and cryospheric processes, which can hardly be represented by a single dynamical index. Climatologically, the NAO, ONI, EPNP, and AO indices are found to have a significantly influence on the NH mid-to-high latitudes.
In Figure 11, the authors show the time series of these indices in 2000-2019 during which only SIE exhibits a significant trend. Two strong peaks in NAO since 2010 are separated by ~4 years, which are apparently cross-correlated with the oscillations in CI and SC. The cross correlation among indices implies the existence of a significant covariance and mutual dependence on similar Because CF is highly correlated with the SW flux in the subarctic regions, we create a cloudiness index (CI) to represent the gridbox mean CF at latitudes of 55 • N-65 • N using the CERES EBAF data. As seen in Figure 10, this CI index is designed to capture most of the cloud-induced SW variabilities over the subarctic landmasses,

Impacts of SIE, SC, NAO, ONI, CI, EPNP, and AO
Since the EOF analysis cannot provide clear attribution to what physical processes are responsible for the interannual oscillations seen in the TOA flux data, we regress CERES TOA fluxes to several climate indices, namely, SIE, SC, NAO, ONI, CI, EPNP, and AO, to evaluate their connection to the Arctic flux variations. SIE, SC, an CI can directly affect the TOA albedo because of relatively high albedos from sea ice, snow and clouds. Other atmospheric (e.g., water vapor, aerosol) and surface (e.g., wetlands, vegetation, melt ponds) variables may directly or indirectly alter the TOA albedo. Collectively, these albedo contributors are subject to atmospheric, hydrological, and cryospheric processes, which can hardly be represented by a single dynamical index. Climatologically, the NAO, ONI, EPNP, and AO indices are found to have a significantly influence on the NH mid-to-high latitudes.
In Figure 11, the authors show the time series of these indices in 2000-2019 during which only SIE exhibits a significant trend. Two strong peaks in NAO since 2010 are separated by~4 years, which are apparently cross-correlated with the oscillations in CI and SC. The cross correlation among indices implies the existence of a significant covariance and mutual dependence on similar underlying process(es). Another cross-correlation example is shown in Figure 12 where the zonal mean CF has latitude-dependent correlation with the AO and ONI indices. It suggests that CFs in the core of arctic (75 • N-90 • N) are highly correlated with AO while CFs are more correlated with ONI in the subarctic region (65 • N-75 • N).
Despite no significant dominance from each individual index, NAO, CI, and SIE are found to contribute most to the Arctic TOA flux variations, followed by ONI and EPNP ( Table 2). For the TOA SW flux, the designated CI outperforms NAO to make the most contribution. For the surface upward SW flux, which can be obtained from the CERES EBAF data, SIE provides the most contribution (24%) as expected for the dominant surface process associated with the Arctic sea ice loss. The combined contribution from the seven indices (SIE, NAO, CI, SC, ONI, EPNP, and AO) counts for~47% of TOA total and SW flux variance. Similarly, the surface SW flux covered by multiple indices can go up to 56%.
The correlation maps of the TOA fluxes with these indices reveal the regions where each index has most impacts ( Figure 13). For SIE, strong correlations with the total and SW fluxes are seen over the Arctic sea ice margin zone as well as over Greenland and Russia. The LW flux has strong anti-correlation with SIE over Greenland, Russia and Northern Canada, but not much over the sea ice margin zone.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 underlying process(es). Another cross-correlation example is shown in Figure 12 where the zonal mean CF has latitude-dependent correlation with the AO and ONI indices. It suggests that CFs in the core of arctic (75°N-90°N) are highly correlated with AO while CFs are more correlated with ONI in the subarctic region (65°N-75°N). Despite no significant dominance from each individual index, NAO, CI, and SIE are found to contribute most to the Arctic TOA flux variations, followed by ONI and EPNP (Table 2). For the TOA SW flux, the designated CI outperforms NAO to make the most contribution. For the surface upward SW flux, which can be obtained from the CERES EBAF data, SIE provides the most contribution (24%) as expected for the dominant surface process associated with the Arctic sea ice loss. The combined contribution from the seven indices (SIE, NAO, CI, SC, ONI, EPNP, and AO) counts for ~47% of TOA total and SW flux variance. Similarly, the surface SW flux covered by multiple indices can go up to 56%.
The correlation maps of the TOA fluxes with these indices reveal the regions where each index has most impacts ( Figure 13). For SIE, strong correlations with the total and SW fluxes are seen over the Arctic sea ice margin zone as well as over Greenland and Russia. The LW flux has strong anticorrelation with SIE over Greenland, Russia and Northern Canada, but not much over the sea ice margin zone.
NAO produces overall strong positive correlation with the total and SW fluxes over the Arctic Ocean, Greenland and Northern Canada, but anti-correlated with the LW flux in these regions. In the northern Europe, NAO is negatively correlated with the total and SW fluxes, with a slightly positive correlation with the LW flux.       NAO produces overall strong positive correlation with the total and SW fluxes over the Arctic Ocean, Greenland and Northern Canada, but anti-correlated with the LW flux in these regions. In the northern Europe, NAO is negatively correlated with the total and SW fluxes, with a slightly positive correlation with the LW flux.

Interactions between Mid-and High-Latitude Disturbances
The observed TOA radiation changes indicate a significant imbalance of the summer Arctic energy on an interannual-to-decadal time scale. The oscillatory patterns seen in Figure 1 between midand high-latitudes suggest some complex teleconnection or coupled processes at work. The lagged correlation between subarctic and Arctic TOA fluxes indicates potential influences of Arctic sea ice loss on subarctic climate on an interannual scale. To further characterize the coupling between low and high latitude variations, here we apply a two-dimensional (2D) Fast Fourier Transform (FFT) decomposition method to the oscillatory patterns in Figure 1, and derive the poleward and equatorward progressing components. As a typical spectral analysis, this FFT method splits the input 2D wave into standing, y+ (poleward) and y− (equatorward) components, where the x-axis is time. For the August TOA SW flux, the equatorward component has the higher contribution with 41% of the variance, compared to 29% and 28% of the standing and poleward components.
The connection between mid-and high-latitude oscillations becomes much clearer in the decomposed equatorward and poleward components with extended coherent patterns ( Figure 14  It is worth noting that Europe experienced some of the worst heat waves in the summer of 2003 and 2015, which in both cases came 2 years after a strong TOA SW anomaly at the mid-latitudes. How the heat wave patterns are related interannually between different regions warrants a further investigation. For the poleward component, the SW flux and CF have a relatively weak influence from mid to high latitudes, nor significant influence extended from the tropical ENSO.

Conclusions
In this study we compared the TOA albedo observed by MISR and CERES in 2000-2019 over the Arctic, the Antarctic and the HMA/TP (High Mountain Asia and Tibetan Plateau), and found generally consistent distributions and interannual variabilities, despite the fact that the observations were made from very different remote sensing techniques. It remains challenging to make reliable observations of the TOA fluxes from space, because the assumptions used in the retrievals for angular distribution functions (ADMs) and narrow-to-broadband (N2B) conversion can subject to various issues at high latitudes. Scene type classification, complex 3D surface topography, spatiotemporal sampling differences can all contribute to errors or differences of the albedo measurement made by MISR and CERES.
Despite these limitations and different observing techniques, both CERES and MISR data showed consistent TOA albedo distribution and interannual variations over high latitudes, as well as over HMA/TP, which provided high fidelity about the observed patterns. Several salient features emerge from these measurements: • A summertime TOA albedo decreasing trend is consistently seen in both MISR and CERES data at latitudes greater than 60°N where the Arctic change is dominated by

Conclusions
In this study we compared the TOA albedo observed by MISR and CERES in 2000-2019 over the Arctic, the Antarctic and the HMA/TP (High Mountain Asia and Tibetan Plateau), and found generally consistent distributions and interannual variabilities, despite the fact that the observations were made from very different remote sensing techniques. It remains challenging to make reliable observations of the TOA fluxes from space, because the assumptions used in the retrievals for angular distribution functions (ADMs) and narrow-to-broadband (N2B) conversion can subject to various issues at high latitudes. Scene type classification, complex 3D surface topography, spatiotemporal sampling differences can all contribute to errors or differences of the albedo measurement made by MISR and CERES.
Despite these limitations and different observing techniques, both CERES and MISR data showed consistent TOA albedo distribution and interannual variations over high latitudes, as well as over HMA/TP, which provided high fidelity about the observed patterns. Several salient features emerge from these measurements: • A summertime TOA albedo decreasing trend is consistently seen in both MISR and CERES data at latitudes greater than 60 • N where the Arctic change is dominated by the summer sea ice loss. The lowest Arctic albedo for JJA occurred in 2011 during the two-decade observation period. We further evaluated ERA5 and MERRA2 TOA SW fluxes with the MISR and CERES observations. It was found that at high latitudes there are large discrepancies between the CERES/MIRS TOA albedo and ERA5/MERRA2 reanalyses, suggesting a serious shortcoming in these advanced reanalyses in terms of representing polar energy budget. These reanalyses failed to capture the lagged relation between the Arctic and subarctic TOA SW fluxes, and in some cases they produced an out-of-phase variations compared to the observations. The observed fluctuations of Arctic TOA fluxes have a quasi-4-year period during 2000-2019 and appear to intensify after 2009. Sea ice and cloud variabilities dominate year-to-year variations of the high-latitude TOA albedo. In the Arctic the CI, NAO, and SIE contribute most to the interannual variation, for which the PC variation of the first EOF is highly correlated with the NAO index. In the Antarctic, MERRA2 lacks the TOA SW variability induced by sea ice changes during 2000-2019, whereas ERA5 generally performs better in representing this variability.
We applied an 2D-FFT analysis to decompose the interannual variations of the Arctic TOA fluxes in the time and latitude domain. The method helps to reveal a clearer polar influence on the subarctic climate, as well as the lagged relationship between the mid and high latitudes.