Earth’s Albedo and Its Symmetry

The planetary albedo  , an intrinsic property of a planet, measures the fraction of incident radiant energy (or insolation) that it reflects back to space. Its complement, the co-albedo   (1 ) thus determines what fraction of that insolation, I , remains to heat the planet. Several components of the planet and interactions among them go into deciding the value of  . For Earth, the atmosphere and clouds are major contributors to albedo (Ramanathan, 1987), as are its surface properties (the land fraction, ice cover (Budyko, 1969) and even the biosphere (Betts, 2000)). Whereas the contributions of the constituent parts of the albedo have been studied in great detail, little attention has been devoted to understanding the properties of the albedo as a whole, as seen from space, and as one might do for another planet. Abstract The properties of Earth's albedo and its symmetries are analyzed using twenty years of space-based Energy Balanced And Filled product of Clouds and the Earth's Radiant Energy System measurements. Despite surface asymmetries, top of the atmosphere temporally & hemispherically averaged reflected solar irradiance R appears symmetric over Northern/Southern hemispheres. This is confirmed with the use of surrogate time-series, which provides margins of   2 0.1 0.28Wm for possible hemispheric differences supported by Clouds and Earth's Radiant System data. R time-series are further analyzed by decomposition into a seasonal (yearly and half yearly) cycle and residuals. Variability in the reflected solar irradiance is almost entirely (99%) due to the seasonal variations, mostly due to seasonal variations in insolation. The residuals of hemispherically averaged R are not only small, but also indistinguishable from noise, and thus not correlated across hemispheres. This makes yearly and subyearly timescales unlikely as the basis for a symmetry-establishing mechanism. The residuals however contain a global trend that is large, as compared to expected albedo feedbacks. It is also hemispherically symmetric, and thus indicates the possibility of a symmetry enforcing mechanism at longer timescales. To pinpoint precisely which parts of the Earth system establish the hemispheric symmetry, we create an energetically consistent cloud-albedo field from the data. We show that the surface albedo asymmetry is compensated by asymmetries between clouds over extra-tropical oceans, with southern hemispheric storm-tracks being 11% cloudier than their northern hemisphere counterparts. This again indicates that, assuming the albedo symmetry is not a result of chance, its mechanism likely operates on large temporal and spatial scales.

Clouds and Earth's Radiant System (CERES) (Kato et al., 2018;Loeb et al., 2018) datasets provide precise measurements of Earth's radiant energy balance as seen from space. From these measurements it is possible to deduce the magnitude of Earth's planetary albedo,   0.29, which varies surprisingly little across years (Stevens & Schwartz, 2012). Remarkably, the measurements show that on long-time averages the two hemispheres have the same albedo, which we refer to as the hemispheric albedo symmetry (Bender et al., 2017;Haywood et al., 2016;Stephens et al., 2015Stephens et al., , 2016Stevens & Schwartz, 2012;Voigt et al., 2013Voigt et al., , 2014. In the simplest approximation this arises from the asymmetric distribution of clouds that counterbalances surface albedo asymmetries (Voigt et al., 2013). As pointed out by Stevens and Schwartz (2012) (see also the review by Stephens et al. (2015)), these properties of Earth's albedo lack even the outline of a theoretical explanation, yet are fundamental to an understanding of Earth's climate. Consider the possibilities: (a) the albedo symmetry is by chance and no mechanism enforces it. Then the two hemispheres are two different "realizations" of a planet with the same shortwave energy balance, and can be used to test how models capture circulation and cloudiness at the largest scales. This also means that, due to conservation of energy, cross equatorial heat transport must be balanced by the hemispheric asymmetry in outgoing longwave radiation, which would constrain inter-hemispheric temperature gradients, but how? Now if (b), the symmetry is not by chance, then we need to know the mechanisms behind it.
There is a surprising scarcity of literature on the topic. Early theoretical studies of Earth's albedo mostly focused on surface contributions to the planetary albedo, specifically the ice-albedo feedback first postulated by Arrhenius (1896). Work on this question flourished in the early 1970s, see Budyko (1969); Sellers (1969); North (1975); Chýlek and Coakley (1975); Stone (1978) among others. These studies sought to understand how the surface albedo depends on, and in turn influences, the planetary temperature, but they did so by taking clouds for granted. In effect, they assumed that, despite being by far the dominant and most dynamic component of the planetary albedo (as later shown in pioneering work by Ramanathan, 1987), secular changes in cloudiness are negligible. The early work also assumed, albeit implicitly, hemispheric symmetry. Such an assumption seems natural, until one asks why clouds should vary between the hemispheres in a way that counterbalances a large, 6 Wm −2 , surface radiation asymmetry (see below), yet somehow be independent of a temporally changing climate state.
Our work is motivated by the intellectual tension that arises from the assumed constancy of clouds over long timescales on one hand, and the fact that they compensate the surface albedo asymmetry to result in overall symmetric albedo on the other hand. We build upon earlier observational studies by Donohoe and Battisti (2011); Loeb et al. (2019), who decomposed the top of the atmosphere (TOA) albedo value and time-anomalies into a surface and atmosphere contribution, and by Voigt et al. (2013), who showed that the hemispheric albedo symmetry is neither a trivial property of the Earth system, nor reproduced by comprehensive Earth system models. By using the new cloud-information provided in the latest release of CERES, we are able to create a new measure of cloudiness that allows us to better understand how, and to what extent, cloud variations influence albedo variations. Our analysis is further aided by more sophisticated methods of time-series analysis, and a near doubling of the length of the observational record as compared to what was available to Voigt et al. (2013). This allows a more rigorous quantification of properties of Earth's albedo that models or theories should explain, and provides context for observations from an increasing number of studies of the albedo of other planets (Cowan & Agol, 2011;Mansfield et al., 2019;Shields et al., 2013), some of which also identify a role for clouds (Kreidberg et al., 2014).
Our contributions are as follows. We first establish the margins of possible observable hemispheric differences as allowed by CERES data, which are 0.1  0.28 Wm −2 . The margins include 0, confirming the conjecture of "symmetry" stated in existing literature. We further analyze radiation time-series in search of indicators of dynamical communication mechanisms that establish the symmetry between the two hemispheres. The overwhelming majority of temporal variability is associated with the seasonal cycle. Residuals of the seasonal cycle are found to be indistinguishable from noise and as such provide no sign of a symmetry-establishing mechanism operating on short (monthly) timescales. However, a long-term trend in the albedo is shown to be hemispherically symmetric, which we would not expect in the absence of a symmetry-establishing process, hinting the existence of one at longer timescales. We then construct a "cloudiness" field, representing physical cloud albedo, which is a better representation of clouds' impact on the shortwave part of the radiation balance than cloud fraction. Using this "cloudiness" we demonstrate DATSERIS AND STEVENS 10.1029/2021AV000440 2 of 13 that, as expected, the hemispheric albedo symmetry is a result of a hemispheric cloud asymmetry (see below). The cloudiness allows us to further demonstrate the extent to which average, versus spatiotemporal, cloud properties are responsible for observed cloud asymmetries. Unexpectedly we find that the major source of asymmetry is in the last place one would expect: storms over the southern ocean are on average cloudier than their counterparts in the north. In the conclusion section we put our findings together and discuss which ones favor and which disfavor the existence of a symmetry-establishing mechanism, and what properties it should have if it exists.

Terminology
Notation and symbols used in the text are summarized in Table 1. The reference to clear-sky adopts values defined by CERES and is estimated from scenes identified as being cloud free. "All-sky" denotes no sub-sampling of specific scenes. The  symbol is used to indicate that a result has been rounded to the displayed digits.
To account for the eccentricity of the Earth's orbit, as well as the uneven sampling of the orbit done by monthly averages,  weights every month by the number of days in that month. We additionally ensure that only full years (total time points multiple of 12) participate in the average. Not doing this can give wrong results. For example,    ( ( )) ( ( )) 0.002 I I     Wm −2 , but using standard averaging of the first 240 months (20 full years currently available in CERES EBAF) instead of weighted gives  0.6 Wm −2 . Such large differences are reported in the literature, e.g. Figure 4 of (Stephens et al, 2016), but are inconsistent with Kepler's second law, which has the direct result that each hemisphere receives the same amount of total insolation over a full orbital period.
By definition,   R I, with  the "albedo". For the energy balance we mostly care about the portion of insolation that is scattered back to space, which we will call effective albedo and define as / R I. It can only be estimated when  0 I . We use the term physical albedo to distinguish the case when  is estimated in some other manner, and thus can be defined also for the case  0 I .

Albedo's Value and Hemispheric Symmetry
The effective albedo is  Clouds thus increase the planetary albedo (on average) by almost a factor of two. The hemispheric difference in R is, on average, ). The same difference for the clear-sky is 6 Wm −2 , 11% of its global average. While the value R = 0.1 Wm −2 appears small, we want to quantify the range of possible values of R supported by the CERES data when one considers statistical variability. This range is useful even in the case that the albedo hemispheric symmetry is by chance, as it can help bound how sensitive clouds are to large scale environmental changes that exist between the two hemispheres. We adopt two approaches to provide distributions of possible values of R, and justification on why these approaches are valid will be made clear in §2.2.
For the first approach, we make use of surrogate time-series (Lancaster et al., 2018;Theiler et al., 1992) to approximate ( ) R  and ( ) R  . To construct the surrogates we adopt the method of Small et al. (2001), which is adapted to periodic data without any correlation structure between cycles. In essence these surrogates sub-sample the months with replacement (i.e., for all Aprils, some are repeated in the surrogates, some are DATSERIS AND STEVENS Abbreviation: TOA, top of the atmosphere.

Table 1
Notation Used in This Paper shuffled around, while some others are skipped entirely). An example is shown in Figure 1a. The hemispheric difference of temporal averages of the surrogates, as estimated from several thousand surrogates of ( ) R  and ( ) R  , gives a distribution of possible R values shown in Figure 1c.
For the second approach, we decompose hemispherically averaged R into a component Y containing the seasonal cycle and a residual for Northern Hemisphere (NH), similarly for Southern Hemisphere (SH), with the seasonal component defined as with  i the chosen frequencies of the decomposition, which are the annual (12 months) and semi-annual (6 months) cycle. The fitting parameters,  Through the decomposition (a),  0, We now estimate the uncertainty range of R that arises simply because CERES provides a finite measurement span (such uncertainty would exist even in the assumption of an active symmetry-establishing mechanism). To quantify this range, we analyze the fluctuations of We model D as an auto-regressive process (Brockwell & Davis, 1996) with  white noise (with same standard deviation as D), M the auto-regressive order and  i the parameter choices. We estimate  i that best fit the measured time-series with the Levinson method of linear predictive code (Levinson, 1946), using  12 M . Figure 1b presents D and a realization of t D . Using t D we simulate a very long autoregressive time-series (that in the limit   t has 0 mean by definition) and sample 242-long subsets of it and calculate their mean. The resulting distribution is shown in Figure 1d.
Both distributions of (Figures 1c and 1d) are very similar and lead us to the same conclusions. First, observing   0 R Wm −2 is a likely scenario, as it is well contained in the 5%-95% quantiles of the distributions. More interestingly though, the distributions provide the range of    0.1 0.28 R Wm −2 as obtained using the 2 of distribution Figure 1. Given that the mean is much smaller than its uncertainty, this establishes the hemispheric albedo symmetry. This property of R was noted in the very first space-based measurements (von der Haar & Suomi, 1971), its systematic and quantitative study only became possible with the advent of the CERES data (Bender et al., 2017;Haywood et al., 2016;Stephens et al., 2015Stephens et al., , 2016Stevens & Schwartz, 2012;Voigt et al., 2013Voigt et al., , 2014.

Dynamics of Residuals
In this and the following subsections we analyze the radiation time-series further in search of an indication of dynamics, or communication, that leads to the hemispheric albedo symmetry, assuming one exists. If there were some process that established the hemispheric albedo symmetry, we might be able to identify it in the dynamics of the residuals  R  and  R  . For instance if one hemisphere responded to internal anomalies that develop in the other hemisphere we would expect to see this in the correlation structure between the two time-series.  1.1 Wm −2 for either hemisphere, very little signal is carried by the hemispheric residuals. Also, apart from a consistent downward trend in both hemispheres, no correlation exists between the time-series of  R in the two hemispheres (checked via a mutual information analysis).
To test for the possibility of dynamics in the time-series of  R we first try to reject the null hypothesis that it can be described by a linear Gaussian (stochastic) process. To do so we first compare the different quantiles of the data with those from a Gaussian distribution fit to the same data. The result falls almost exactly on the diagonal (Figure 3a), which means that the original data can be modeled well from this distribution. Similarly, but not shown, a K-sample Anderson-Darling test (Scholz & Stephens, 1987) also fails to reject the null hypothesis. This establishes that the distribution of the data is Gaussian, but not whether their sequence is correlated, or just noise.
To explore this latter question we create a surrogate time-series for  R  (we find identical results for SH, not shown) following the approach of Nakamura et al. (2006), which is designed to represent fluctuating data with constant trends. The surrogate follows a linear Gaussian process with autocorrelation the same as the input time-series, and same trend, and thus satisfies the null hypothesis by construction ( (Brockwell & Davis, 1996)). We then attempt to distinguish this time-series from the actual time-series of  R  . To do so a 1-step self-mutual-information statistic is adopted as a discriminatory statistic q., as it has been shown to distinguish noise from determinism in small data sets (Lancaster et al., 2018). The distribution of values ( ) p q can be constructed from different realizations of the surrogate. If the real data are significantly outside this distribution, the null hypothesis can be rejected. As Figure 3b shows, the real q is well within the possible q of the null hypothesis, thus failing to reject the null hypothesis.
The apparent lack of dynamics in the residuals of either hemisphere may simply be telling us that we have projected-through averaging-a too high dimensional dynamics onto a too low dimensional space (and thus the result is indistinguishable from noise due to extreme information loss). To test for this possibility we repeated the above analysis on 10°  10° longitude-latitude decompositions and reached similar conclusions. On yet smaller scales, Voigt et al. (2013) showed that albedo features become more strongly correlated with neighboring areas, and thus are no longer independent samples. A principle component analysis also failed to identify dominant patterns of variability.
In summary, we were unable to identify any evidence of dynamics in the residual time-series or any cross-hemisphere correlation, outside the statistically significant trend shared by both hemispheres that we discuss more in §2.4. This result provides insight on a symmetry establishing mechanism, assuming one exists. It seems unlikely that its effect is visible in the yearly timescale (which is what we can meaningfully analyze statistically given the short timespan of CERES). This could be because such a mechanism operates on timescales well beyond the year (see also discussion in §2.4), or, even if it does operate on short timescales, it affects shortwave reflection only indirectly. For example, energy transfer across the equator would have major impact on the fluctuations of air temperature or major circulation patterns (Kang et al., 2008;Voigt et al., 2014), but those do not directly tune the shortwave reflection. Rather, over longer timescales they can tune the overall amount of cloudiness and thus increase or decrease cloud albedo accordingly. And indeed, we did a simple check to see whether the hemispheric difference of the residuals DATSERIS AND STEVENS 10.1029/2021AV000440 5 of 13  of R can be related with those of surface temperatures, but we found no mutual information between the two.

Temporal Variability
The seasonal decomposition as applied in Equation 1 removes seasonal fluctuations of R that are due to seasonal fluctuations of physical albedo (e.g., the melting of ice during summer) as well as those directly due to I. To separate the two effects we calculate the temporal average value of ,   . This is done for each hemisphere, after hemispheric averaging, and the results are shown in Figure 4. In essence, solar R is the R we would observe if we replaced Earth with a completely static, time-averaged version of itself.
Even in this decomposition, which disentangles the internal and solar fluctuations, the insolation accounts for most of the variability of R. Specifically for the NH, 84% of the variance of R is attributed to solar R , only 1% to internal R and 13% to their co-variability. For the SH the numbers are 68%, 3% and 28%. A lack of power in Fourier components other than those corresponding to the seasonal cycle (1 and 0.5 years, Figure 4b) in internal R is consistent with our earlier analysis, which failed to distinguish  R from noise.

Secular Trends and Hemispheric Co-Variances
The lack of signal in  R makes the magnitude of its secular trend surprising. From the best-fit of the residuals we estimate its magnitude to be −0.7 Wm −2 per decade. This can't be explained by a trend in ( ) I  , which is thought to arise from the the 11-year solar cycle and the shortness of the CERES record, which began near a solar maximum and ends near a solar minimum. That trend is only  −0.05 Wm −2 , hence  ( ) I  can only explain about 2% of the observed trend in  ( ) R  . Because of the miss-match of trends in incoming and outgoing solar radiation, we assume that these trends have physical origin and are not due to calibration errors (the possibility of which has been discussed extensively with the CERES team, but for which we have not been able to find any evidence, see also (Loeb et al., 2020)).
To develop a sense of the magnitude of the trend in  ( ) R  , we compare its value to roughly 0.2 K per decade rise in globally averaged surface temperatures since 2000. If the former were attributed to the latter it would imply a positive albedo feedback greater than  0.7 (1-0.02)/0.2 = 3.4 Wm −2 K −1 . This is a factor of four or more larger than assessed values (Sherwood et al., 2020 report a central estimate of 0.75 Wm −2 K −1 ), which is to say it is a large number. The trend, first identified and analyzed by Loeb et al. (2020), was attributed to changes in Northeast Pacific stratocumulus forced by decadal variations in sea-surface temperatures, in a way that models appear to largely capture. What hasn't been previously noted, and what we find difficult to explain given the lack of dynamics in  R , is why a trend in  ( ) R  attributed to changes in stratocumulus in the Northeast Pacific, is so well mirrored by much less spatially coherent (Loeb et al., 2020), but equal, changes in  . Were a substantial (more than a fourth) part of the observed trend attributable to global warming, it would have dramatic consequences. Put more broadly, if usefully quantifying the pace of global warming requires an ability to understand and quantify cloud responses to warming on the order of 0.2 Wm −2 K −1 , then surely more effort, building on prescient work by Loeb et al. (2020), to understand tenfold larger trends, and their symmetry, is warranted.
Despite our inability to distinguish detrended values of  R from noise, the similarity of the trends in  R  and  R  makes a case for hemispheric communication. We further analyze the spatial variability of R on sub-hemispheric scales. We split each hemisphere into two halves (quadrants) separated by the great circle aligned with the longitude . We then calculate the average difference of R between unique (four) combinations of quadrants as a function of . Differences between arbitrary quadrants are, on average, much larger than hemispheric differences as we show in Figure 5. This analysis demonstrates (see also (Voigt et al., 2013)) that large scale differences in R can occur, even after averaging over all latitudes. That the DATSERIS AND STEVENS 10.1029/2021AV000440 6 of 13 differences between semi-spheres within a hemisphere are much larger than those between semi-spheres across hemispheres could arise by chance, but this seems unlikely. We hypothesize that this difference is the signature of large-scale constraints on cloudiness; and while this signature could arise from an intrahemispheric constraint, why it would act independently in each hemisphere to maintain symmetry in R is difficult to understand.

Cloudiness and Albedo
In this section, we quantify the contribution of cloudiness to  to investigate how cloud asymmetries compensate hemispheric asymmetries in  K , the cloud-free albedo.

Defining Cloudiness, C
We begin by creating a "cloudiness" field C that represents physical cloud albedo by combining two independent definitions of cloud albedo. These two definitions qualitatively agree with each other, and give us the confidence that our results do not depend (qualitatively) on the exact definition of cloudiness. Combining them allows us to take advantage of their differing properties: one is energetically consistent with I and R, the other can exist for  0 I .
The first definition for cloudiness is the cloud contribution to atmosphere albedo. Donohoe and Battisti (2011) provide equations that can decompose the TOA albedo into an additive contribution from the atmosphere and the surface, based on the radiative fluxes at the TOA and surface and assuming that the atmosphere can be approximated as a single, uniform layer described by a given albedo and transparency, having the same scattering properties regardless if radiation comes from upwards or downwards.
An advantage of  C , as compared to using the "cloud radiative effect" to define the cloudiness (as  ( ) / R K I), is that it does not conflate cloud with surface variability. Notice that  C is effective and not physical albedo and is valid only for  0 I .
The second definition comes from a cloud albedo parameterization defined as with f the cloud area fraction,  the cloud optical depth, and g the asymmetry factor from the cloud particle phase function. Equation 4 is the same as Equation 19 of Lacis and Hansen (1974), but multiplied with f . CERES provides f , and  , but we have to estimate g (see below). The field  has missing values, at the high-latitudes of the winter hemisphere. In Appendix A we describe a regularization process to fill these values, and also discuss the impact of using different observational data. Using constant  0.9 g for all grid points already gives very good qualitative agreement between   , C C in both temporally averaged maps but also in spatially averaged time-series.
Choosing g so that results in a physical cloud albedo,  C , that is energetically consistent with the time-averaged effective cloud albedo contribution,  C , everywhere, but at the expense of g varying spatially. Because the spatial variations are small and within a physical range (see Appendix A), we adopt this approach, and denote the resultant cloud field by C. The reason for combining both definitions, by creating a spatially varying g, is that it allows us to define an energetically consistent cloud field that is also defined when  0 I . In addition, by testing that our results are qualitatively similar for two independent DATSERIS AND STEVENS 10.1029/2021AV000440 7 of 13  Figure A2.

Mean Cloudiness
Temporally and zonally averaged distributions of C are presented in Figure 6 and compared to the cloud fraction f (normalized so that it has same mean and standard deviation as C). This shows that simply using f to represent cloudiness overestimates the impact of clouds in the deep tropics (equatorward of about 20°). As measured by C, the tropics are substantially less cloudy than the extra-tropics. Despite large latitudinal variations in cloud regimes-with the intra-tropical convergence zone being mostly north of the equator -C varies little (by less than 15%) with latitude within the tropics.
. Hence the SH is cloudier by  0.02 -about 12% of the global average value. This asymmetry in C is (see Figures 6 and 8) largely a result of the SH extra-tropics being much cloudier than the NH extra-tropics. This asymmetry is already evident in early cloud climatologies (Arrhenius, 1896;Brooks, 1927) and was further pointed out by (Bender et al., 2017). Its precise quantification here shows its importance for compensating the hemispheric asymmetry of the surface albedo.
In principle, asymmetries in C need not imply asymmetries in the effective albedo (and hence R). For example, were the asymmetries carried by nocturnal cloudiness they would have no effect on R. Qualitatively however, the asymmetry in cloudiness is sufficient to establish the hemispheric albedo symmetry (i.e., on first order the spatio-temporal characteristics of C don't matter). Because the hemispheric difference of C represents an effective albedo value, and because the clear-sky albedo difference is     ( ( )) / ( ( )) ( ( )) / ( ( )) 0.02 , then the hemispheric difference of C is sufficient to counter that of K. Using the model of Donohoe and Battisti (2011), increasing the atmospheric contribution to TOA albedo by 0.02 units gives a total TOA albedo increase of  0.02 for a wide range of choices of surface albedo, atmosphere albedo and atmosphere transmittance. This does not prove that hemispheric differences in the mean cloudiness compensate asymmetries in K, but suggests that it could.

Correlation With the Annual Cycle
In addition to hemispheric differences in the temporal and spatially averaged cloudiness, another candidate for compensating asymmetries in K is asymmetries in the covariability between I and C. In Figure 7a we show hemispherically averaged time-series of C and in Figure 7b the cross-correlation function of the hemispherically averaged time-series of cloudiness C and insolation I. C is strongly linearly correlated with I (maximum cross-correlation values are  1), which is not surprising since cloudiness has a strong annual cycle. What we did not expect is the different delays d in the maximum of the cross-correlation (which is the phase shift between C and I) between the two hemispheres.
For the SH C has its maximum approximately at the time of maximum insolation as  0.3 d months is small. In the NH,  1.5 d months is larger, and appears to mostly be attributable to clouds in the NH extra tropics (not shown), for reasons that remain unclear. The larger lag of cloudiness in the NH contributes a small but measurable contribution to the compensation of the clear-sky asymmetry. Shifting the cloudiness time-series by 1.2 months (the difference of the delays between NH and SH), increases the cloud reflection by 0.28 Wm −2 -or about 10% of what is required to balance the asymmetry in  K . DATSERIS AND STEVENS 10.1029/2021AV000440 8 of 13 and cloud area fraction f (normalized for same mean and std. as C). We also display C averaged over four equal-area latitude zones (colored areas).

Cloudiness Over Different Surface Types
Having established that extra-tropical asymmetries in C largely compensate for hemispheric asymmetries in  K , here we ask whether these are predominantly associated with clouds over a particular type of surface. For instance, is the NH extra-tropics less cloudy by virtue NH land masses being less cloudy? We investigate this possibility by defining different surface types: O denotes the ice-free ocean area fraction, E the ice & snow coverage and    1 L O E the ice-free land fraction. All three of these are spatiotemporal fields and , O E are accessible from CERES auxiliary datasets. To estimate cloudiness over the three different surface types here we use each one of , , O E L as statistical weights. For each point in time, and for each hemisphere, we perform a weighted spatial average (in addition to the standard weighting by area) of the field C with weight , O E or L. This gives us the average value of C over a specific surface type as a time-series. We then perform a temporal average  and present the results in Figure 8. Differences between the two hemispheres can be directly compared with the average hemispheric difference of C of  0.02. Taking the co-variability of clouds and ice into account mostly impacts the results for C over ice, which makes sense given that E varies much more than , L O.
A trivial explanation for a compensating cloud contribution to the hemispheric clear-sky asymmetry   0.02 K would be that land is less cloudy than ocean, in a way that directly compensates for its increased surface albedo relative to the ocean. To test this hypothesis, we have to calculate the difference of cloudiness under the assumption of same cloudiness for a given land type regardless of hemisphere. Neglecting contributions from ice-covered surfaces, which Figure 8 shows to be almost identical between hemispheres, we need to solve the equation . We also need to explicitly write (approximately true for land, not true at all for ocean). We arrive at the expression , we obtain that, in order for the albedo symmetry to be explained exclusively from the fact that land is less cloudy than ocean, then the difference in cloud albedo between land and ocean should be . This is much larger than the real differences shown in Figure 8.
Instead, the true compensation of surface albedo asymmetries comes from the clouds over SH oceans being more reflective than those of NH oceans (as seen by combining the information in Figures 6 and 8). We found this surprising, and despite it having long been appreciated that the SH extra-tropics are cloudier than the NH, it has generally been assumed that this follows from the storm-tracks in the south occupying a greater area, rather than being brighter per unit area.

Conclusions
We use twenty years of CERES data to study the global properties of Earth's planetary albedo and reflected solar irradiance R. The hemispheric albedo asymmetry, defined as the difference between the temporal and hemispheric averages of R is estimated as 0.10(28) Wm −2 . The uncertainty (two sigma) is estimated using advanced time-series analysis techniques and indicates that the measured asymmetry is indistinguishable from zero and provides a margin of error models should satisfy. Further analysis indicates that the yearly cycle contains most variability of R, while its residuals are indistinguishable from noise, at least in the yearly timescales, and thus not correlated across hemispheres. The residuals of both hemispheres share identical decadal trends of −0.7 Wm −2 per decade, which argues for a communication mechanism at longer timescales.
Our analysis of cloud albedo C (mostly combining Figures 6 and 8) identifies the global tropics (equatorward of 30°) as surprisingly transparent, with a zonally and temporally averaged value of C of 0.12 varying little with latitude. The extra-tropics are nearly twice as cloudy, the SH (0.24) substantially more so than the NH (0.20). Differences between cloudiness in the NH and SH are primarily found over the ocean. Whereas land is less cloudy than the ocean (0.15 vs. 0.17), the differences are insufficient to compensate for DATSERIS AND STEVENS 10.1029/2021AV000440 9 of 13 differences in land versus ocean clear-sky reflectances. Where past work (Voigt et al., 2014) has investigated the role of shifts in tropical clouds as a means of symmetrizing the hemispheric albedo, the CERES data indicates it is the asymmetry between the southern and northern hemispheric extra-tropical storm tracks that is most responsible for compensating hemispheric asymmetries in the cloud-free planetary albedo. The broad dynamical similarity between extra-tropical storms in the two hemispheres (Kodama et al., 2019) makes this difference in cloudiness unexpected. If anything, microphysical arguments would lead one to expect more ready rain initiation and less cloudiness in the SH storms.
In this work, we did not provide a mechanism that establishes the albedo symmetry nor argue that one necessarily exists. Instead, we provide evidence favoring both sides of the argument, list properties that such a mechanism should satisfy, and exclude certain timescales and candidate processes as candidates. Let's assume that a symmetry-enforcing mechanism exists. Then, we would either expect the residuals of R to not be noise and be correlated across hemispheres, or the mechanism operates on longer timescales than the yearly and sub-yearly (only these timescales can be meaningfully analyzed statistically with the methods of §2.2). However, it cannot happen at ultra-long timescales, because then this wouldn't explain the hemispherically symmetric trends, which have a timescale of 10 years. On the other hand, the symmetry could be by chance. But we find this attribution weak, because it is hard to attribute to luck both the mean value of R being hemispherically symmetric, as well as its trends. The cloud analysis makes matters more intriguing. A reasonable first candidate for a mechanism would have it associated with shifts in tropical rain bands (Voigt et al., 2014), which modeling studies show are more closely associated with heat transport between the hemispheres (Kang et al., 2008). That the asymmetry in cloudiness is counter-balancing the clear-sky albedo asymmetry is mostly confined to clouds over the extra-tropical ocean seems to argue against this idea, favoring the hypothesis that a communication mechanism does not exist, or suggesting that if it does exist, energy transfer on longer timescales that connects tropics and extratropics (e.g., by the ocean) might play a role.

Appendix A: Observational Data Processing for Cloudiness
To derive the atmospheric and surface contributions to the albedo  we re-write the model of Donohoe and Battisti (2011)  the surface albedo which in general is different than the surface contribution to the planetary albedo due to further reflections between surface and atmosphere. F simply stands for shortwave radiation, and all necessary instances of F are provided in the same CERES data set used throughout. Note that Stephens et al. (2015) adopted a similar approach, but arrived at slightly different expressions, which appear to be incorrect, or incorrectly typeset. The optical depth field  provided by CERES has missing values in a large portion of its time-series for spatial points near the poles. These missing values make hemispheric averages impossible and thus need to be "fixed". Here, we used a simple sinusoidal continuation, shown in Figure A1. Using the same method as in §2.1, we fit sinusoidals of frequencies 1/year and 2/year to the available (i.e., nonmissing) data. The value of the sinusoidal fit is used to fill in only the entries of  that are missing. Furthermore, Figure A2 shows temporally averaged maps of  , , , C f g.
Different satellites and observational data products can have significant differences in both how they define and measure cloud area fraction f and optical depth  , this will cause differences in the diagnosed values of g chosen to maintain consistency between  C and C. In addition, we found it difficult to assess the degree to which the signals we identify are free of retrieval biases. Specifically, the large differences in C between winter and summer, as well as between tropics and extra-tropics, could be influenced by the influence of sun-sensor geometry on the retrievals. This question merits further investigation with the retrieval experts to ensure that the signals we identify are physical.

Data Availability Statement
Analysis is based on the monthly averaged CERES data. The authors use the Energy Balanced and Filled (EBAF) Ed. 4.1 (Kato et al., 2018;Loeb et al., 2018), data set. For cloud properties we use the SYN1deg data (Doelling et al., 2013;Rutan et al., 2015). We also use the auxiliary data of SYN1deg to get ice-free ocean and ice & snow area fractions. All data have spatial resolution of 1°  1° and span from March 2000 to April 2020, totaling 242 months. All of these datasets are publicly available. This work is also available as a fully reproducible code base online (Datseris, 2021