The EPIC images constitute a time series of the sunlit part of the Earth, including the atmosphere, starting from June 2015 and continuing at the present time. There are usually about 22 images per day during Northern Hemisphere summer, and 13 during winter. Each multispectral image has 10 wavelength channels between 317 and 780 nm. The channels have full-width-at-half-maximum (FWHM) values between 1 and 3 nm (Ohtake et al., 2010; Ohtake et al., 2013; Geogdzhayev and Marshak, 2018; Herman et al., 2018).

Each image pixel represents the radiance reflected by the corresponding area of Earth. This radiance is a function of the reflective properties of the area, and the solar and satellite angles. The DSCOVR satellite is located in the first Lagrange point between the Earth and the Sun, therefore the radiance is observed close to the backscattering geometry with satellite and solar angles being almost equal. To convert measured narrowband radiances at backscattering into integrated broadband albedo values at top-of-atmosphere (TOA; see Narrowband-to-Broadband and Radiance-to-Radiosity Transformations for more details), we apply the ADMs provided by the CERES project (Su et al., 2015a; Su et al., 2015b). The ADMs are provided for several surface types and cloud fractions, and we combine these into three types for our analysis: clear land, clear ocean, and cloud-covered areas. The reason for this is that the angular resolution of the tabulated ADMs available at the CERES project web pages ² is only 10°, and we want to interpolate with a finer resolution in the backscattering direction and combining similar surface types gives us more robust estimation (see Scattering Geometry and Interpolation of Angular Distribution Models). Additionally, estimating temporary cloud fractions from EPIC images is challenging, a more robust method can be developed for simple clear versus cloud-covered area estimation. Land or ocean surface classification is available from the International Geosphere–Biosphere Programme, and temporal cloud coverage we estimate from the EPIC images using a logistic regression model with input from EPIC channels at 325, 551, and 780 nm (see Cloud Classification in EPIC Images).

The TOA albedo values for each pixel over the Earth’s sunlit disk are averaged for albedo of the Earth at each narrowband wavelength channel. Finally, the narrowband albedos are summed with weights from incident solar spectra at each channel, obtained from NOAA’s climate data record of solar spectral irradiance (Coddington et al., 2015; see Narrowband-to-Broadband and Radiance-to-Radiosity Transformations). The solar spectrum is taken as constant in this study with only the distance between the Earth and the Sun influencing the total flux level.

The daily averaged time series has been filtered to remove outliers. For some EPIC images, not all the wavelength channels are present. If there is not a proper coverage of usable images over the day, some parts of Earth are not present, and this would introduce a bias to the mean albedo of that day. Therefore, these days are left out from our time series data. There is a period between mid-2019 and early 2020 when the EPIC camera was not operational.

We derive the mean yearly shortwave spherical albedo of the Earth by grouping the daily time series values by the day in the year, averaging per day, and finally averaging over days in a year.

The EPIC images contain 10 narrowband channels, namely at 317, 325, 340, 388, 443, 551, 680, 688, 764, and 780 nm. The FWHM values for the channels are between 1 and 3 nm. The channels are measuring the directional radiance. On the other hand, we are interested in the broadband value of the TOA integrated radiosity, i.e., reflected flux integrated over all directions, and the broadband albedo, i.e., the ratio of the reflected flux to the incident flux.

The ADMs provide the conversion from directional broadband radiance to TOA integrated radiosity. They are applied for example in Su et al. (2018), Su et al. (2020) for broadband radiance that is first derived with a narrowband-to-broadband conversion model from the EPIC images. In this work we first apply the ADMs and second the narrowband-to-broadband conversion, but we would like to point out that both of these conversations can be described together with a single linear transform and that the order of narrowband-to-broadband and radiance-to-radiosity is interchangeable.

After the radiance-to-radiosity conversion with the ADMs we combine the narrowband channels into a broadband estimate. We use 8 of the 10 EPIC channels for this conversion. Channels at 688 and 764 nm are measuring the O₂ absorption and are left out of the conversion, since they do not describe the reflectance outside these narrow channels. The solar spectral irradiance (Coddington et al., 2015) is convoluted with the response curves of the eight EPIC channels to give weights for the narrowband contributions that are summed together for the broadband estimate.

As the EPIC channels range from 317 to 780 nm, strictly speaking this is the coverage of our broadband estimate. The portion of the solar input in this range is 52.7%. In ultraviolet (UV) below 317 nm there is only 0.5% of the total energy, so we can safely extend the validity of our broadband range into UV. The CERES ADMs, on the other hand, are derived for shortwave broadband range of 300–5000 nm. The use of ADMs to the broadband with wavelength range of 317–780 nm only can be a source of bias in our estimate. However, Su et al. (2018) reports successful usage of ADMs with only three EPIC channels (443, 551, and 680 nm) for narrowband-to-broadband transformation.

The algorithm for cloud classification in the EPIC images was trained with a manually prepared data set of spectral radiances together with subjective cloudy/clear label decision. For image data, we collected images from the first week of every month in 2018. We converted the original HDF5 library data for one multi-channel image into PNG image stack of all the wavelength channels. We opened the image stacks in ImageJ ³ , and manually selected areas of either clear or cloud-covered land and ocean surfaces. The spectral radiance values of pixels within these areas were recorded together with the label (1) clear land/(2) clear ocean/(3) cloud-covered area. Altogether, we produced 338 areas of clear land, 331 areas of clear ocean, and 481 areas of cloud-covered areas. Each area contains several tens of pixels, the mean values over the areas were employed in the analysis. The ice sheet covered areas, such as Antarctica, are indistinguishable from cloud-covered areas.

A logistic regression model was separately fitted into data with clear land and cloud-covered areas, and data with clear ocean and cloud-covered areas. A logistic regression model has a linear function modeling the log-odds of the probability of an event. In our case, the event is that the surface is covered with clouds. The model can be written as p = 1 1 + e − ( β 0 + β 1 x 1 + ⋯ + β k x k ) where p is the probability of cloud-cover, and β 0 + β 1 x 1 + ⋯ + β k x k is the linear function with unknown coefficients β i and known radiances x i from channels 1 , … , k .

We executed forward-selection stepwise regression to optimize the model with the best value of the Bayesian Information Criterion (BIC) statistics of the model. We ended up with models with a constant coefficient β 0 and three coefficients for the EPIC wavelength channels at 325, 551, and 780 nm. For land vs. clouds model, all the four coefficients were tested significant with p-values less than 0.1%. For ocean vs. clouds model, all other but the coefficient for the channel 780 nm were significant with the same p-value limit. The values for the model coefficients are given in Table 1. The final classification is done with the probability limit of 1 / 2 : if p ≤ 1 / 2 the pixel is classified as cloud-covered. The channel values in this stage are in the counts/second units and are not yet converted into directional radiance.

We found it realistic to successfully classify pixels of clear or cloudy land or ocean areas from the EPIC images, but not different cloud types (water or ice) or levels of cloud cover that are present in CERES ADMs. ADMs are also divided into different wind speeds for ocean and into different vegetation types for land. To exploit these, wind speeds on ocean pixels would need to be connected from weather models into our analysis together with up-to-date land cover information. We feel that this would complicate the albedo estimation too much at this stage and could be a topic for future development. We decided to emphasize more the interpolation of the scattering geometry inside the ADM models over having many ocean/land/cloud subtype ADMs.

The scattering geometry in the EPIC observations is such that the observation is always done close to backscattering. The phase angle between the Sun and the DSCOVR spacecraft, as seen from Earth, is below 12° except for few rare cases in years 2020–21. On average, phase angle has been 8.2°. After the break in the DSCOVR operations in 2019, the minimum phase angle has gradually started to decrease (see Figure 7). It was never below 4° before 2020, but now the minimum value is 1.8°.

While in terms of phase angle EPIC is observing close to backscattering, the local solar zenith angle on each Earth pixel varies between 0° and 90°. The ADM tables, on the other hand, have 10° bins for all the three angles (solar zenith, viewing zenith, azimuth). So, our backscattering geometry implies that we need ADM table values inside cells where the solar and viewing zenith angles are in the same bin and azimuth angle is in the bin 0–10°. Thus, the 3-dimensional (in angle) ADM tables will reduce for our purposes to 1-dimensional tables that can be parametrized with the solar zenith angle only.

We find the reduced 1-dimensional table binning of 10° too coarse as it would introduce inaccuracies when applied to pixels having the solar zenith angle anywhere between 0° and 90°. Therefore, we interpolate the reduced ADM tables over the 0° and 90° range using a cubic spline interpolation. We use these interpolated coefficients when converting the TOA radiances measured by EPIC into integrated TOA albedos. We find that the reduced 1-dimensional tables show relatively similar behavior for all clear ocean surfaces regardless of the wind speed, for all clear land surfaces, and for all cloud-covered areas regardless of the cloud type or coverage. Thus, our final set of ADM coefficients come from interpolated data on these three categories, averaged over their subtypes.

The annual variation of the Earth’s spherical (Bond) albedo is described for the first time in detail. Previous investigations mainly describe the variations on the incident solar flux due to the Earth’s elliptic orbit (see Figure 5 in Kandel and Viollier, 2010), as for us, this effect is considered and properly accounted for Su et al. (2018), Su et al. (2020) are using also EPIC images, but they do not reduce the albedo but produce outgoing fluxes. The accurate annual albedo behavior enables the detailed analysis of local geographic and atmospheric effects on the albedo during a year.

The Earth contains regions of low albedo (e.g., cloud-free ocean areas, vegetated land areas) and high albedo (e.g., clouds, ice and snow surfaces). Diurnal and annual variations in the albedo time series (Figure 1) result from a modulation of these two aspects, depending on the apparent longitude (diurnal variations) and latitude (annual variations). Especially, the cloud-covered ocean areas increase the albedo. The daily cloud-covered ocean fraction, estimated from the EPIC images, has the correlation coefficient value of 86% with the daily albedo value (see Figure 7).

In annual albedo variations, the main albedo maximum occurs in December, around Southern solstice, when the Antarctic ice sheet, sea ice, and snow cover are visible entirely. At this time, cloud formations of the mid-latitude cyclones over the Southern Ocean are well pronounced. Albedo is further enhanced by the shallow convective cloud cover over the subtropical oceans and the relatively small areas of cloud-free ocean areas visible, on average, at this time of the year.

The secondary albedo maximum occurs in June, around Northern solstice, when the Greenland ice sheet and sea ice in the Arctic Ocean are well exposed. Cloud formations of the mid-latitude cyclones of the North Atlantic and Pacific storm tracks are active, although reduced from their winter maxima. Deep convective clouds of the inter-tropical convergence zone are on the Northern Hemisphere at this time of the year. Again, relatively little cloud-free ocean areas are visible.

The main albedo minimum precedes the main maximum in August–September when neither of the Polar regions is visible. At this time of year, there are relatively few shallow convective clouds over the subtropical Northern Atlantic and Pacific, as well as over the Mediterranean. North Atlantic and Pacific storm track activities are near to their minima. Also, vegetation is at its maximum extent in the Sahel region. Monsoon is in active phase over the Indian peninsula and South-East Asia, thus enhancing the albedo. The secondary minimum precedes the secondary maximum in March–April.

The solar energy input at the Earth’s average distance from the Sun is 1361  W / m 2 for the sunlit disk of the Earth over all wavelengths (Coddington et al., 2015). The EPIC filters range from 317 to 780 nm in wavelength and the portion of the solar input between these wavelengths is 52.7%. We can make a rough estimate that the albedo derived here is valid for somewhat longer wavelengths, say to 60% of the reflection or absorption of the total solar input. With this estimate, the input energy that is absorbed into the area of sunlit Earth’s disk and atmosphere is about 4.9  W / m 2 more during the albedo minimum in September than on average. During the albedo maximum in December, about 8.7  W / m 2 less energy is absorbed.

An exceptional albedo occurred on 13 December 2020 and persisted for about a week before and after the peak value. This anomaly coincides with the main annual maximum in December. Our main explanation is a short-lived and exceptional maximum in the amount of shallow convective clouds over the subtropical Indian Ocean and the Pacific east of Australia, see Figure 5. The development of exceptionally wide shallow cloud cover coincides and is favoured by a constellation of several mature mid-latitude cyclones, feeding moist air from the subtropics towards the Antarctic. At the time of the maximum albedo, all relevant large-scale atmosphere-ocean indicators were in close-to-neutral stages ⁴ : the Madden-Julian oscillation (MJO) phase centre was in the Oceanic continent, but the amplitude was weak; the Indian Ocean Dipole (IOD) was neutral; El Niño–Southern Oscillation (ENSO) was in a weak La Niña phase (El Niño 3.4 index) and returning to neutral.

In December 2020, the EPIC-derived daily albedo was elevated to 0.334 from its typical December value of 0.312 in 2015–2019. Assuming these to represent well the planetary reflectivity, the emission temperature of the Earth would decrease by 2.1 K through a 7.8  W / m 2 reduction in the absorbed solar radiation. For comparison, the global mean sensible heat flux from the Earth surface to the atmosphere amounts to about 20  W / m 2 and the absorbed solar radiation in the atmosphere to about 80  W / m 2 (Hartmann, 2015). It is therefore fair to treat the short-term albedo peak in December 2020 as quite remarkable.

The phase angle of EPIC observations varies between 2° and 12°. After the DSCOVR operations break in 2019, the phase angle variation has increased, and the minimum phase angle has decreased from previous 4°–1.8° (Marshak et al., 2021). These decreasing phase angles introduce a possible source of bias in our albedo estimation. The resolution in the binned data of the ADM models that we use is 10° in all three angles, indicating that we do not have exact information on how possible backscattering effects (self-shadowing, coherent backscattering) are behaving on small phase angles inside this bin.

Currently, we can only estimate the possible bias resulting from small phase angles. We can do this by computing the deviations of daily albedo values from the average-year albedo values, and by modeling how they are changing as a function of the observation phase angle. From Figure 6 we can see that there is an obvious correlation between these, the linear correlation coefficient being −0.59. Decreasing phase angles are introducing a positive bias, on average, into the albedo values. A second-degree linear fit has the p-values of all the model coefficients significant below 0.1% level, and a slightly better value (0.37) of adjusted R 2 measure of the model than with the first-degree model (0.36). According to this model, there can be up to 0.01 positive average bias coming from the phase angle effect if the phase angle approaches 2°.

This same bias can also be roughly estimated from the small upturns in daily albedo values in November 2020 and March 2021, see Figures 1, 7. These upturns are simultaneous with local phase angle minima of about 2°. Graphical estimation of possible excess in albedo values in these upturns also suggests a value of 0.01 at most.

The exceptional albedo values in December 2020 were at most 0.023 larger than the average December value. This exceptional value was received at the time of local phase angle minima of the DSCOVR spacecraft, 2.1°. We conclude that the effect of the phase angle can explain, at maximum, about half of the difference, and that December 2020 is exceptional in our albedo time series even when taking the possible phase angle effect into account.