The solar energy flux reaching the ocean surface in the spectral range 400–700 nm, referred to as photosynthetically available (or active) solar radiation (PAR), controls the rate of photosynthesis by phytoplankton and therefore the development of crustaceans, fish, and other consumers (e.g., Ryther, 1956; Platt et al., 1977; Kirk, 1994; Falkowski and Raven, 1997). It ultimately regulates the composition of marine ecosystems. Sunlight absorbed differentially by the upper ocean affects mixed-layer dynamics and oceanic currents (e.g., Nakamoto et al., 2000, 2001; Murtugudde et al., 2002; Sweeney et al., 2005; Ballabrera-Poy et al., 2007), with local and remote consequences on atmospheric temperature and circulation (e.g., Miller et al., 2003; Shell et al., 2003). Absorption by phytoplankton and other water constituents tend to reduce the planetary albedo, i.e., warm the planet (Frouin and Iacobellis, 2002). Knowing the spatiotemporal distribution of PAR over the oceans is critical to understanding biogeochemical cycles of carbon, nutrients, and oxygen and biological-physical interactions (a major uncertainty in coupled climate models) and, therefore, to addressing important global change issues such as the fate of anthropogenic atmospheric carbon dioxide and making accurate projections of future climate (e.g., Frouin et al., 2018a).

Regional and global maps of PAR at the ocean surface can be obtained from a variety of passive Earth-viewing satellite optical sensors. The sensors operating from geostationary altitude provide adequate temporal sampling to deal with cloud diurnal variability but have degraded spatial resolution at high latitudes, and they only cover part of the oceans, i.e., several sensors, optimally positioned are necessary to provide global coverage. Sensors in polar orbits provide the same spatial resolution at all latitudes but pass less frequently over the same target at middle and low latitudes. For ocean primary productivity computations, it is convenient to estimate both PAR and bio-optical variables (phytoplankton chlorophyll abundance, absorption coefficients) from the same sensor. Ocean-color sensors offer this capability, even though they are principally designed to retrieve water reflectance if they do not saturate over clouds. The same data preprocessing is required, i.e., PAR can be produced with little extra effort as part of the same processing line. In this way, the key variables in primary production modeling are provided together at the same resolution, facilitating studies of photosynthesis and ecosystem dynamics.

In this context, a simple yet efficient and fairly accurate algorithm has been developed to estimate the daily mean PAR at the ocean surface from Sea-viewing Wide Field-of-view Sensor (SeaWiFS) data (Frouin et al., 2003) and adapted for application to MODerate resolution Imaging Spectroradiometer (MODIS) data (Frouin et al., 2012), GLobal Imager (GLI) data (Frouin and Murakami, 2007), GOCI data (Frouin and McPherson, 2013; Kim et al., 2016), Medium Resolution Imaging Spectrometer (MERIS) data, Visible Infrared Imaging Radiometer Suite (VIIRS) data, Second-generation Global Imager (SGLI) data, and Advanced Hiwamari Imager (AHI) data with plans for an extension to future ocean color sensors. Daily mean PAR refers to the 24-h averaged planar quantum energy flux from the Sun in the spectral range 400–700 nm. It is expressed in units of Einstein per meter squared per day, i.e., Em⁻²d⁻¹. The global daily mean PAR products from SeaWiFS, MODIS, VIIRS, and MERIS data have been routinely generated by the National Aeronautics and Space Administration (NASA) Ocean Biology Processing Group (OBPG) and made available to the user community from their website (https://oceancolor.gsfc.nasa.gov). Estimated uncertainty, based on comparisons against in situ measurements, expressed in relative root-mean-square (RMS) difference and bias, is typically 10-30% and 4-9%, respectively, depending on satellite sensor and atmospheric conditions (Frouin et al., 2003; Frouin et al., 2012; Laliberté et al., 2016; Ramon et al., 2016; Somayajula et al., 2018). Somayajula et al. (2018) compared satellite-based PAR algorithms used in primary production studies; they concluded that the best overall performance was obtained with the NASA OBPG algorithm. This uncertainty is reasonable for large-scale studies of aquatic photosynthesis (e.g., Frouin et al., 2012; Frouin et al., 2018a), but better accuracy is desirable. Note, in this respect, that primary productivity models depend not only on PAR but also on efficiency factors that are difficult to estimate with uncertainty comparable to (i.e., as low as) that of PAR.

The standard Level-2 and -3 PAR products generated by the NASA OBPG have been used extensively in the science community for a variety of applications. In primary productivity calculations, they have replaced PAR estimates obtained from a clear sky model corrected for cloudiness using fractional cloud coverage or deduced from satellite estimates of total solar irradiance, the treatment applied in Longhurst et al. (1995), Antoine et al. (1996), and Behrenfeld and Falkowski (1997) to obtain the first global maps of seasonal and/or annual oceanic primary productivity from space. Such treatment is limited, because the effect of clouds on PAR does not depend only on fractional coverage, but also on optical thickness, and the relation between total solar irradiance and PAR, rather constant under clear skies (Baker and Frouin, 1987), varies strongly with water vapor and cloud liquid water content (Frouin and Pinker, 1995). Studies using the NASA OBPG PAR products have addressed a variety of topics, including biosphere productivity during an El Niño transition (Behrenfeld et al., 2001), chlorophyll-a and carbon-based ocean productivity modeling (Behrenfeld et al., 2005; Platt et al., 2008), climate-driven trends in productivity (Behrenfeld et al., 2006; Kahru et al., 2009; Henson et al., 2010), phytoplankton class-specific productivity (Uitz et al., 2010), inter-comparison of productivity algorithms (Carr et al., 2006; Lee et al., 2015), and the relation between primary productivity, vertical mixing, and atmospheric input (Tang and Shi, 2012). They have also been used to check the stability of CERES measurements (Loeb et al., 2006).

The parameters governing PAR variability are essentially the Sun zenith angle and the cloud transmittance. Aerosol properties and surface albedo have a smaller impact. Since the Sun zenith angle can be computed precisely, estimating daily PAR from data collected by a single sensor aboard a Sun-synchronous satellite is chiefly limited, in terms of accuracy, by the lack of information about diurnal variability of cloud properties, especially at low and middle latitudes. This variability may be large in some regions, as evidenced by the International Cloud Climatology Project (ISCCP) cloud analyses (Bergman and Salby, 1996; Rossow and Shiffer, 1999) and other studies (e.g., Wang and Zhao, 2017; Zhao et al., 2019; Yang et al., 2020). Consequently, the PAR products from individual polar-orbiting sensors exhibit biases, not only with respect to ground truth but also between themselves, as evidenced in inter-comparison and evaluation studies (Frouin et al., 2003, 2012; Tan and Frouin, 2019). Merging data from several sensors with different overpass times may significantly improve the quality of daily PAR estimates, as demonstrated with MODIS-Terra, SeaWiFS, and MODIS-Aqua, which cross the equator at approximately 10:30, 12:00, and 13:30 local time (Frouin et al., 2012). In generating a long-term PAR time series, however, one must deal with various sensor combinations, and there is a need, for ocean biogeochemistry studies related to climate change to reduce the individual biases against in situ measurements and make the PAR estimates consistent across individual sensors (Frouin et al., 2018a).

The Earth Polychromatic Imaging Camera (EPIC) onboard DSCOVR, operating from the first Sun-Earth Lagrange point (L1) one million miles from Earth (Marshak et al., 2018; https://avdc.gsfc.nasa.gov/pub/DSCOVR/Web_EPIC/), provides a great opportunity to generate accurate PAR products and address issues associated with polar-orbiting sensors. By frequently observing the sunlit part of the Earth, EPIC inherently allows one to account properly for diurnal cloud variability, while maximizing spatial coverage. In other words, EPIC with respect to PAR can do the job of several geostationary sensors with the further advantage that spatial resolution at high latitudes is less of an issue (the L1 orbit is much farther from Earth than the geostationary orbit). The spectral bands centered on 443, 551, and 680 nm, the non-saturation of measured radiance over clouds, and the spatial resolution of 10 km at nadir are adequate for PAR calculations, especially using the NASA OBPG algorithm, which does not require knowing whether the pixel is clear or cloudy, i.e., is applicable to large pixels.

In view of the above, the current NASA OBPG daily mean PAR algorithm has been modified/adapted for application to EPIC data. Algorithm uncertainties have been associated with EPIC PAR estimates on a pixel-by-pixel basis. A full processing line has been created and implemented to generate operationally daily mean EPIC PAR products at the NASA Center for Climate Simulation (NCCS). The data are archived at and distributed by the Langley Atmospheric Science Data Center (ASDC). In Section 2, the methodology to estimate daily mean PAR from EPIC data is presented and the various steps to obtain the surface flux values are detailed. The tasks include integrating atmospheric functions spectrally and temporally during the day (the number of observations in a day varies depending on geographic location), eliminating data contaminated by Sun glint, incorporating ancillary information such as ozone content, sea ice extent (for masking), and aerosol optical properties, and remapping the data to a common grid. In Section 3, a procedure is described to associate algorithm uncertainties (i.e., bias and standard deviation) to each EPIC daily mean PAR estimate as a function of parameters readily available from applying the algorithm, i.e., daily mean clear sky PAR and cloud factor (characterizes the effect of clouds on daily mean PAR). In Section 4, EPIC PAR estimates are compared to in situ measurements routinely collected from long-term fixed buoys. Experimental performance is also compared to that of MODIS PAR estimates. In Section 5, examples of global daily mean PAR products are displayed and examined in view of corresponding MODIS and AHI products, and PAR time series at contrasted locations are presented to illustrate the capability of EPIC to describe PAR seasonal to interannual variability. In Section 6, finally, the EPIC PAR algorithm and its performance against field data and other satellite estimates are summarized, advantages and limitations of using observations from the L1 orbit are pointed out, the significance of the new ocean PAR product in complementing existing PAR time series for a wide range of research applications is emphasized, and a perspective for future work to estimate variables used more directly in primary productivity or water reflectance models, such as scalar PAR, spectral PAR, and average cosine of the light field just below the surface, as well as ultraviolet fluxes, is provided.

The algorithm estimates daily mean PAR reaching the ice-free ocean surface, as defined above. Following Frouin et al. (2003), a budget approach is used, in which the solar flux reaching the surface is obtained by subtracting from the flux arriving at the top of the atmosphere (know) the flux reflected to space (estimated from the EPIC measurements) accounting for atmospheric transmission and surface albedo (modeled). Clear and cloudy regions within a pixel do not need to be distinguished, which is appropriate to the relatively large (i.e., 10 km at nadir) EPIC pixels. This approach was shown to be valid by Dedieu et al. (1987) and Frouin and Chertock (1992).

Based on the previous work, the PAR model assumes that the effects of clouds and other atmospheric constituents are decoupled. The planetary atmosphere is therefore modeled as a clear sky layer that contains molecules and aerosols positioned above a cloud/surface layer, and surface PAR is expressed as the product of a clear-sky component and a transmittance that accounts for cloudiness and surface optical effects. Under solar incidence (zenith angle) θ s , the incoming spectral solar flux at the top of the atmosphere, E T O A c o s ( θ s ) , is reduced by a factor T a ( θ s ) T g ( θ s ) [ 1 − S a A ( θ s ) ] − 1 by the time it enters the cloud/surface system, where T a is the clear-sky transmittance (due to scattering by molecules and aerosols), T g is the gaseous transmittance (essentially due to absorption by ozone), S a is the spherical albedo of the clear atmosphere, and A is the cloud/surface system albedo. As the transmitted flux, E T O A c o s ( θ s ) T a ( θ s ) T g ( θ s ) [ 1 − S a A ( θ s ) ] − 1 , passes through the cloud/surface system, it is further reduced by a factor 1— A . The instantaneous planar flux reaching the ocean surface at any wavelength in the PAR spectral range, E , is then given by: E ( θ s ) = E 0 ( θ s ) [ 1 − A ( θ s ) ] [ 1 − A s ( θ s ) ] − 1 [ 1 − S a A ( θ s ) ] − 1 , (1) where A s   is the albedo of the ocean (surface and water body) and E 0 ( θ s ) = E T O A c o s ( θ s ) T a ( θ s ) T g ( θ s ) is the solar flux that would reach the surface if the cloud layer and the surface were non-reflecting. In clear sky conditions, A reduces to A s and E to E c l e a r = E 0 ( 1 − S a A s ) − 1 . The term [ 1 − S a A ( θ s ) ] − 1 in Eq. 1 represents interactions between the cloud/surface layer and the clear atmosphere.

To compute E , A is expressed as a function of the radiance measured by EPIC in bands centered on 443, 551, and 680 nm, and this for each observation collected during the day (i.e., between sunrise and sunset). These bands do not saturate over clouds, and they sample sufficiently the PAR spectral range. In this range, the scattering properties of molecules and aerosols vary smoothly with wavelength, gaseous absorption is relatively weak, and A is quite constant spectrally. The algorithm works on a pixel-by-pixel basis, i.e., a daily mean PAR estimate is obtained for each pixel of the satellite imagery. The various processing steps and procedures are detailed in the following.

First the bidirectional reflectance of the cloud/surface layer, ρ , is determined from the instantaneous EPIC Level 1b reflectance data at source resolution. This is only accomplished for pixels that are not contaminated by sea ice or Sun glint. For sea ice masking, fractional ice coverage at 25 km resolution from the National Snow and Ice Data Center (NSIDC) is used. If the fractional coverage is greater than 0.1, then the pixel is discarded. For Sun glint masking, the Fresnel signal at the wavy interface surface is calculated from the wind speed according to Cox and Munk (1954), and all pixels with a glint reflectance above 0.05 at 780 nm are eliminated. National Centers for Environmental Prediction (NCEP) Reanalysis-2 wind data at 1^o and 6-h resolution is used after interpolation to the EPIC observation time. Both ancillary data sets are remapped to the source grid prior to masking. Since in the modeling the clear atmosphere is located above the cloud/surface layer, ρ is obtained from the top-of-atmosphere (TOA) reflectance, ρ T O A by inverting the following equation (Tanré et al., 1979): ρ T O A ( θ s , θ v , ϕ ) = T g ( θ s , θ v ) [ ρ a ( θ s , θ v , ϕ ) + ρ T a ( θ s , θ v ) ( 1 − S a ρ ) − 1 ] , (2) where ρ a is atmospheric reflectance and θ v and ϕ are view zenith and relative azimuth angles, respectively. In Eq. 2 T g and T a are transmittances along the Sun-to-surface and surface-to-sensor path. This expression is strictly valid for an isotropic target of reflectance ρ , which is not the case for the cloud/surface layer, but a good approximation for many geometries, except when θ s   is large and clouds are optically thin, i.e., when the albedo of that layer, A , is relatively close to ρ (Zege, et al., 1991). At large θ s , however, E becomes small since modulated directly by c o s ( θ s ) , and when clouds are thin their impact on ( 1 − A ) in Eq. 1 is also small, which tends to minimize effects on daily mean surface fluxes. A better treatment would be to use the following equation (Tanré et al., 1979; Deschamps et al., 1983): ρ T O A ( θ s , θ v , ϕ ) = T g ( θ s , θ v ) { ρ a + [ ρ ( θ s , θ v , ϕ ) − A ( θ s ) ] e − τ a m + A ( θ s ) T a ( θ s , θ v ) [ 1 − S a A ( θ s ) ] − 1 }   , (3) where τ a is the optical thickness of the clear atmosphere and m is airmass, i.e., m = 1 / cos ( θ s ) + 1 / cos ( θ v ) . Eq. 3 reduces to Eq. 2 when ρ is isotropic. The retrieval of ρ would be accomplished through iteration. This would require a first estimate of A , which can be obtained as described below for isotropic ρ .

To compute the atmospheric functions T g , T a , and ρ a , the surface pressure and water vapor amount are extracted from NCEP Reanalysis 2 data, aerosol optical thickness at 550 nm, Angström coefficient, and single scattering albedo at 550 nm from MERRA-2 data at 0.625^o × 0.5^o (Gelaro et al., 2017), and ozone amount from EPIC Level 2 data at source resolution (Herman et al., 2018). Single scattering albedo is assumed constant over the PAR spectral range. For ρ _a , the quasi-single scattering approximation (Q-SSA) used in Frouin et al. (2003), which gives too high values at large zenith angles, is replaced by a parameterization that combines Q-SSA and exact single-scattering approximation (SSA): ρ a = f S S A ( m , τ a e r ) ρ a ( S S A ) + f Q ‐ S S A ( m , τ a e r ) ρ a ( Q − S S A ) , (4) where f S S A and f Q − S S A depend on m and aerosol optical thickness, τ a e r , but are independent of the aerosol model. These functions were obtained from the simulations with the 6S code (Kotchenova et al., 2006, 2007) for atmospheres containing various mixtures of maritime and continental aerosols and angular geometries. Figure 1 displays ρ a at 443, 551, and 680 nm as a function of θ v for θ s   = 70^o. Scattering angle, quasi constant with EPIC, is fixed at 171.5^o. The aerosols are of maritime type and τ a e r at 550 nm is 0.2 and 0.4 (top and bottom panels, respectively). The parameterization works well yielding ρ a values very close to the 6S values except in the blue, where τ a e r   is large (0.4) and the molecular scattering is effective. At 443 nm the coupling between aerosol and molecule scattering is relatively important and the 6S ρ a is substantially smaller than ρ a ( Q – S S A ) , even at θ v < 65^o, which is partly explained by the fact that in Q-SSA the molecule and aerosol contributions to atmospheric reflectance correspond to the atmosphere containing only molecules or aerosols, i.e., the neglected aerosol-coupling term is negative (Deschamps et al., 1983). The 6S and SSA values are quite similar at that wavelength, a consequence of the SSA formulation that also assumes either molecules or aerosols, i.e., the probability of encountering one type of scatterer is enhanced (SSA values are expected to be lower in the presence of both molecules and aerosols). As wavelength increases, this effect is less prominent because the atmosphere is less thick optically contributing to larger differences between 6S and SSA values.

Once the reflectance of the cloud/surface layer, ρ , is determined, it is converted into albedo A . This is accomplished for each observation during daytime and the three spectral bands by applying a cloud bidirectional correction factor F (independent of wavelength) to ρ   −   A s since A   ≈   N A c   +   A s where A c is cloud albedo and N is fractional cloud cover (see Frouin and Chertock, 1992; Frouin et al., 2003), i.e., A   =   F ( ρ   −   A s )   +   A s , (5)

Analytical formulas developed by Zege et al. (1991) for optically thick scattering layers in the non-absorbing medium are used for F . The parameterization depends on the asymmetry coefficient of the cloud indicatrix, g , and the cloud optical thickness, τ c , which are fixed at 0.853 and 15, respectively. Sensitivity of F on g and τ c is an issue, especially when τ c is relatively small (ρ exhibits more directionality), but since E depends on ( 1   –   A ) the relative impact on E computed from Eq. 1 is reduced in such situations ( A is also small). The g and τ c variables can be viewed as tuning parameters, e.g., to reduce biases in estimated fluxes. The surface albedo A s   is parameterized as a function of the Sun zenith angle, aerosol optical thickness at 500 nm, and wind speed following Jin et al. (2004). In the A s calculation, the chlorophyll-a concentration is fixed at 0.1 mgm⁻³, and the wind speed is from NCEP Reanalysis-2 data.

Next, the daily mean PAR, < E P A R > , is computed by integrating over the PAR spectral range, i.e., 400 to 700 nm, and the length of the day: < E P A R >   =   ( 1 24 ) ∫ t [ c o s ( θ s ( t ) ) G ( t ) ] d t w i t h   G   =   ∫ λ   [ E T O A T g   T a   ( 1   –   A ) ( 1   −   A s ) − 1 ( 1   −   S a A ) − 1 ] d λ , (6) where E _TOA is the extraterrestrial spectral solar flux per unit of wavelength corrected for Earth-Sun distance variation during the year and time t is expressed in hour. The integral is calculated for each EPIC observation during the day. In other words, if after masking a surface target (pixel) is observed n times during the day, the algorithm generates n daily mean PAR values. The dependence of A and A s on θ s is accounted for in the time integration, but assuming that the characteristics of the atmosphere and surface are unchanged during the day. This is a crude assumption, especially for cloudiness, but diurnal variability of the atmosphere is considered implicitly in the final step of the algorithm (see below).

In the final step, the individual daily mean estimates obtained on the source grid (number varies from 1 to 13 depending on geographical location, the time during the year, and data availability) are first remapped to an 18.4 km equal-area grid and weight-averaged using the cosine of the Sun zenith angle and then remapped to a Plate Carrée (equal-angle) grid with 18.4 km resolution at the equator. The remapping algorithm is exactly the one used by NASA OBPG to generate a Level 3 binned ocean color products (https://oceancolor.gsfc.bnasa/gov/docs/format/l3bins). Triangular-based linear interpolation is used to fill missing pixels at the edges.

The weighting procedure to obtain the final < E P A R > estimates is accurate, except when only one or two EPIC observations per day is available. Such situations are infrequent; they are due to instrument problems and mission activities (e.g., sensor calibration, spacecraft maneuvers). The accuracy of the weighting procedure was checked from radiative transfer calculations with hourly MERRA-2 data for cloud and aerosol properties and ozone and water vapor contents. In the calculations, E P A R was expressed as E P A R = E c l e a r P A R ( 1 − N A c ) , where A c   is obtained from τ c   according to Fitzpatrick et al. (2004) and E c l e a r P A R   is simulated with the 6S code assuming a chlorophyll-a concentration of 0.5 mgm⁻³ and a wind speed of 5 ms⁻¹. Two dates, i.e., January 1 and July 1, 2018, and ocean locations evenly distributed in the spatial domain within 60°S to 40°N and 180°W to 180°E, i.e., every 5° along latitude, and every 6.25° along longitude, were selected. Table 1 displays the comparison statistics of cos ( θ s ) -weighed < E P A R > obtained from 1, 2, 3, 6, 9, and 12 individual < E P A R > estimates during the day versus the < E P A R > value obtained by trapezoidal integration of all the hourly E P A R values at the time of the MERRA-2 observations (referred to as actual or theoretical < E P A R > ). The individual estimates correspond to hourly observations randomly distributed between sunrise and sunset. In estimating < E P A R > for each observation, the aerosol and cloud properties and gaseous absorber amounts are assumed unchanged during the day, but the dependence of A c and A s on θ s are accounted for in the time integration, as in the EPIC daily mean PAR algorithm. The cos ( θ s ) -weighted < E P A R > values agree well with the theoretical ones when 3 and more observations per day are used, with biases less than 0.10 (0.4%) Em⁻²d⁻¹ in magnitude and root-mean-square difference (RMSD) less than 3.83 (14.2%) Em⁻²d⁻¹. RMSD is reduced to 1.34 (4.5%) and 1.07 (3%) Em⁻²d⁻¹ with 9 and 12 observations per day. The biases are small because of the random sampling of times during the day and the large number of locations considered in the data ensemble. RMSD is noticeably larger when 1 or 2 < E P A R > estimates are used, i.e., 7.00 (25.9%) and 4.86 (18.0%) Em⁻²d⁻¹, respectively. In such situations, the alternative procedure to obtain < E P A R > , i.e., trapezoidal integration of hourly E P A R values, would also give inaccurate results.

A way to reduce the sampling biases in such situations, not yet implemented in the algorithm, is to use MERRA-2 hourly cloud products for the very day of the EPIC observations, as proposed by Tan et al. (2020). If τ c M E R R A and N M E R R A denote the MERRA-2 cloud optical thickness and fractional coverage, A is replaced by   A '   in Eq. 5 as follows:   A ′ = ( A − A s ) [ N M E R R A ( t ) A c ( τ c M E R R A ( t ) ) ] / [ N M E R R A ( t i ) A c ( τ c M E R R A ( t i ) ) ] + A s , (7) where t i is the time of satellite observation and A c is computed from τ c M E R R A at time t and t i . Basically, the quantity A − A s   is adjusted by the ratio of N A c computed from MERRA-2 products at times t and t i In Eq. 6 A , determined at t i , is corrected for variations with θ s , but assuming that cloud properties are unchanged. Figure 2 illustrates the potential improvement expected using Eq. 6. It displays scatter plots of OBPG MODIS-A < E P A R > estimates at source resolution (about 1 km) versus in situ measurements at the Chesapeake Bay CERES Ocean Validation Experiment (COVE) site (36.9^oN, 75.7^oW) in the Western Atlantic during 2003–2014 (see Tan et al., 2020, for details about the in situ data collected at that site). Bias is reduced from 1.72 (5.3%) to −0.05 (-0.1%) Em⁻²d⁻¹, i.e., practically eliminated on average, but RMSD is marginally improved, i.e., 7.02 (21.4%) instead of 7.37 (22.5%) Em⁻²d⁻¹, when using the hourly MERRA-2 cloud products.

Associating uncertainty to each < E P A R > estimate is necessary to use properly the information provided by the estimate. It allows one to verify that variability and trends are actual to merge different datasets optimally, and to constrain adequately model predictions (Frouin et al., 2018a; IOCCG, 2019). Comparisons with contemporaneous and collocated in situ measurements at a few sites, the current approach to quantify experimental performance may not be sufficient to describe uncertainty in the wide range of conditions expected to be encountered. Representing uncertainty on a pixel-by-pixel basis requires modeling the measurement, identifying all possible error sources (e.g., noise in the input variables, imperfect or incomplete mathematical model), and determining the combined uncertainty (JGCM-100, 2008; Povey and Grainger, 2015; IOCCG, 2019). Algorithm uncertainty is due to model approximations and parameter errors (e.g., using plane-parallel atmosphere, decoupling effects of clouds and clear atmosphere, neglecting diurnal variability of clouds, fixing cloud optical thickness to a constant in bidirectional correction) if the input variables (EPIC TOA reflectance at wavelengths in the PAR spectral range) are known perfectly. A complete pixel-by-pixel uncertainty budget should include, not only algorithm uncertainty, but also uncertainty due to measurement noise, radiometric calibration, and preprocessing to Level 1b, a component that is difficult to determine accurately.

The procedure described in Frouin et al. (2018a, b) is used to estimate and provide, for each pixel of the daily mean PAR product, the algorithm uncertainty component of the total uncertainty budget, which is expected to dominate. The bias and standard deviation portions are calculated as a function of clear sky daily mean PAR, < E c l e a r P A R > , and cloud factor, < C F P A R >   =   < E P A R > / < E c l e a r P A R > (characterizes the effect of clouds on PAR and varies from 0 to 1), by simulating for many situations the satellite measurements and corresponding < E P A R > and comparing the latter to the < E P A R > estimated from the TOA reflectance. The simulations were performed with the Atmospheric Radiative Transfer Database for Earth and Climate Observation (ARTDECO) code (Dubuisson et al., 2016; https://www.icare.univ-lille.fr/artdeco/) using as input several years (2003 to 2012) of MERRA-2 hourly data (aerosol and cloud properties). The large number of data points allows one to sample a wide range of geometric configurations for the satellite data and atmospheric variability, in particular many situations of daytime nebulosity for all latitudes. At this stage, the resulting look-up-tables (LUTs) only depend on < E c l e a r P A R >   and < C F P A R > , but other parameters, e.g., angular geometry and latitude, will be considered in the future. Since the final < E P A R > is a cos ( θ s ) -weighted average of individual < E P A R > estimates during the day, the uncertainty on the final < E P A R >   is obtained by weighting the individual uncertainties in the same way.

Figure 3 displays the resulting uncertainty (bias and standard deviation) on individual <E _PAR > estimates as a function of < C F P A R > for several < E c l e a r P A R >   levels, i.e., 12, 35, and 58 Em⁻²d⁻¹. The bias between estimated and simulated < E P A R >   values (Figure 3A) is small (mostly within 1 Em⁻²d⁻¹ in magnitude) and does not vary much with < C F P A R > at low and moderate < E c l e a r P A R >   values (blue and red curves) but reaches 3 Em⁻²d⁻¹ at high < E c l e a r P A R >   values when < C F P A R > is small (black curve). The standard deviation (Figure 3B) is peaked toward intermediate cloud factors (the risk or probability that cloudiness at the time of satellite measurement may not be representative of the conditions at other times during the day is larger), increasing from about 2.2 to 8.7 Em⁻²d⁻¹ when < E c l e a r P A R >   increases from 12 to 58 Em⁻²d⁻¹ and < C F P A R > = 0.5, i.e., about 15% relatively.

As mentioned above, a complete per-pixel uncertainty budget must include errors in the Level 1b data, which may require estimating the sensitivity of < E P A R >   to input reflectance and the covariance in the various spectral bands (JGCM-100, 2008). Noise in the input data is difficult to estimate, especially spectral correlations, but this can be attempted by analyzing the Level 1b imagery using structure functions or variograms (e.g., Curran and Dungan, 1989; Wald, 1989) and considering information gleaned from inflight calibration studies. The approach, therefore, would be to establish the uncertainty LUTs using noisy TOA reflectance simulations. Note, however, that in the case of MERIS comparisons with match-up data at the COVE site have revealed that experimental uncertainties are similar to algorithm uncertainties obtained from modeled data without noise (Frouin et al., 2018b). Similar results are expected for EPIC < E P A R >   uncertainties.

Figures 8A–C displays an example of an EPIC < E P A R > product with associated uncertainties. The date is March 20, 2018 (equinox); the land is in black and the sea ice is in white. The values range from a few Em⁻²d⁻¹ at high latitudes to about 58 E m⁻²d⁻¹ at equatorial and tropical latitudes with atmospheric disturbances modulating the < E P A R > field, especially at middle latitudes (Figure 8A). The bias uncertainty tends to be positive when < E P A R > values are high (overestimation by up to 2 Em⁻²d⁻¹) and slightly negative when < E P A R > values are low (underestimation by up to 1 Em⁻²d⁻¹). The overestimation is relatively larger when < E c l e a r P A R > is high and cloud factor is low, which occurs at low and middle latitudes (Figure 8B). The standard deviation uncertainty is more variable spatially, ranging from 2 to 8 Em⁻²d⁻¹, with highest values obtained at moderate cloud factors and high < E c l e a r P A R > values, as expected from Figure 3, for example off the coast of Chile and Peru at about 25^oS and 80^oW where cloud factor is about 0.5 (Figure 8C).

Compared with the MODIS-Aqua < E P A R > product (Figure 8D), the EPIC product (Figure 8A) is less noisy due to multiple observations during the day and does not have any spatial gaps at low/middle latitudes (gray color). The difference map between the two products (Figure 9A) shows good agreement, with higher differences generally encountered in middle to high latitude regions affected by storm activity. The MODIS-Aqua values cover a slightly larger range, which is expected because more extreme values are likely to be encountered with only one observation per day (Figure 9B). The frequency of values between 15 and 30 Em⁻²d⁻¹ is higher for EPIC reflecting the lower probability of having very low values when several observations during the day are used in estimating < E P A R > in cloudy conditions. The difference histogram indicates that the EPIC < E P A R > values are slightly lower than the MODIS-Aqua values, by 0.7 Em⁻²d⁻¹ on average (Figure 9C). This may result from MODIS-Aqua observing at about 1:30:pm local time, i.e., when cloudiness is usually reduced (e.g., Bergman and Salby, 1996).

The EPIC < E P A R > imagery of March 20, 2018, was also compared with corresponding imagery from AHI onboard Hiwamari-8 (operated by the Japanese Meteorological Agency) over the oceans and seas surrounding East and Southeast Asia and Australia (Figure 10). The AHI < E P A R > product, available at 5 km resolution on the equal latitude-longitude grid from the Japanese Aerospace Exploration Agency (JAXA) was generated from geostationary observations acquired every 10 min using an adapted version of the algorithm described in Frouin and Murakami (2007). It was remapped to the equal-angle grid (18.4 km at the equator) for comparison to the EPIC < E P A R > product. The spatial features are very similar in both products (Figure 10A,B), and differences do not exhibit a distinct pattern, although there is some evidence that larger differences often occur near the edges of low-pressure systems (Figure 10C). This may be due to the different temporal resolution of the two products with AHI capturing more accurately the daily variability of moving disturbances. The range of < E P A R > values is practically the same for EPIC and AHI, but the AHI histogram exhibits a small number of higher values between 5 and 25 Em⁻²d⁻¹ and lower values above 55 Em⁻²d⁻¹ (Figure 10D). Again, this is plausibly attributed to using more observations during the day in the AHI < E P A R > estimation. The < E P A R >   differences are generally small (i.e., mostly between -5 and 5 Em⁻²d⁻¹) but may be as large in magnitude as 20 Em⁻²d⁻¹, resulting in an overall bias (higher EPIC values) of 1.34 Em⁻²d⁻¹ or 3.6% (Figure 10E). Of course, this bias is only valid for the date and area considered in the comparison—it is expected to differ depending on region and period, yet it corresponds to a wide range of atmospheric conditions.

Figure 11 displays the time series of EPIC and MODIS daily and monthly mean E P A R during the entire DSCOVR mission until August 2021 at 50^oN and 30^oW (North Atlantic) and 0^o and 30^oW (Equatorial Atlantic). The MODIS values are averages of MODIS-Aqua and -Terra estimates. The seasonal cycle and day-to-day variability are large at 50^oN, while short-term variability dominates at the Equator. Seasonal E P A R changes are depicted in the same way by EPIC and MODIS, but the MODIS values tend to be slightly higher, as expected. These examples of time series, and the good performance against in situ measurements (Figure 7), illustrate the potential of the EPIC < E P A R > product in ocean biogeochemistry studies.

An algorithm was developed to estimate daily mean PAR at the ice-free ocean surface, E P A R , from EPIC observations in spectral bands centered on 443, 551, and 680 nm. The algorithm, based on Frouin et al. (2003), uses a budget approach that does not require distinguishing whether a pixel is clear or cloudy, which is appropriate for the coarse EPIC pixels. Algorithm uncertainties (bias and standard deviation) were associated to each < E P A R > estimate using LUTs established from RT simulations. A preliminary evaluation showed good agreement with < E P A R > estimates from other satellite sensors (polar-orbiting MODIS and geostationary AHI) and in situ measurements at ocean moorings. Match-up data analysis indicated that the EPIC-derived < E P A R > was less biased than the MODIS < E P A R > , and the EPIC < E P A R > imagery was less noisy, which was explained by using multiple observations during the day with EPIC instead of only one observation with MODIS. The uncertainty of the EPIC < E P A R > product, with biases of −1.5 to 2.2% and RMSDs of 10.0–13.3% depending on the site is lower than the uncertainty of other parameters coming into play in primary production modeling (e.g., phytoplankton absorption and quantum yield). The < E P A R > 2015–2021 time series at selected oceanic locations demonstrated the algorithm ability to capture monthly to interannual variability for investigating the ocean response to temporal changes in available light over a wide range of scales.

The EPIC < E P A R > product is generated routinely by the NASA Center for Climate Simulation (NCCS) and distributed by the Langley ASDC. The data (with documentation) are archived since the beginning of the DSCOVR mission, i.e., June 13, 2015, at https://asdc.larc.nasa.gov/project/DSCOVR/ DSCOVR_EPIC_L3_PAR_01. They are available on equidistant cylindrical Plate Carrée grid with 18.4 km resolution at the equator and in 18.4 km equal-area sinusoidal projection, i.e., on spatial grids that are compatible with MODIS and VIIRS Level 3 NASA OBPG ocean-color products.

The current algorithm can be improved in several ways, i.e., by calculating atmospheric reflectance more accurately at large zenith angles (LUTs may be used instead of approximate analytical representation), by relaxing the Lambertian assumption in the retrieval of the cloud/surface layer reflectance, by improving the parameterization of cloud bidirectional effects, and by including from reanalysis data information about cloud variability, which would provide better accuracy when only a few EPIC observations are available to estimate the daily means. Uncertainty may also be specified as a function of angular geometry and latitude, even region, instead of using an average estimate for all latitudes over several years of MERRA-2 data, and they can be fitted by a generalized additive model with proper auxiliary variables.

Other < E P A R > datasets exist (e.g., OBPG MODIS and VIIRS products), but the EPIC time series is uniquely valuable for several reasons. First, the daily mean PAR estimates, which exploit the multiple EPIC observations from sunrise to sunset (i.e., consider diurnal cloud variability), are more accurate than those from sensors in polar orbit (i.e., typically use one observation per day). Second, coverage is global on a daily time scale at low and middle latitudes, which is currently not achieved with instruments onboard polar orbiters (due to Sun glint, limited swath). Sensors operating from geostationary orbit have limited coverage and reduced spatial resolution at high latitudes, a smaller problem with EPIC. Third, and particularly important, comparisons between daily < E P A R > estimates from EPIC and concurrent polar-orbiting sensors (e.g., MODIS and VIIRS) would determine biases associated with estimates from these sensors, allowing for a consistent < E P A R > time series across sensors, not only during the overlap period, but before and after, i.e., for a long-term (over several decades) science quality < E P A R > record.

The EPIC < E P A R > product is useful to a wide range of research applications, such as primary production and carbon export modeling, ecosystem dynamics and mixed-layer physics, photochemical transformations of dissolved organic matter, and control of stable soluble iron in marine waters. It complements existing < E P A R > datasets and, as mentioned above, may bring about consistency across sensors, allowing a better description of biological phenomena that could lead to new information about temporal variability of biological processes.

The methodology can be easily extended to estimating ultraviolet (UV) surface irradiance using the EPIC spectral bands centered on 317, 325, 340, and 388 nm, especially since ozone content, a key variable governing atmospheric transmittance in the UV, is a standard EPIC product. Furthermore, planar and scalar fluxes below the surface, as well as average cosine for total light (a measure of the angular structure of the light field), variables more directly relevant to addressing science questions pertaining to biogeochemical cycling of carbon, nutrients, and oxygen can also be estimated without major difficulty from the above-surface quantities. Approaches have been identified and procedures devised (Frouin et al., 2018a); they are based on LUTs of clear sky and overcast situations and the derived cloud factor, < C F > ( < C F P A R > for the PAR spectral range), from the EPIC observations (Frouin et al., 2018a). Following Mobley and Boss (2012), the LUTs may only depend on a reduced set of parameters, the most important ones being the location and date, which control the day length and mean Sun zenith angle, then the influence of the clouds which is between null (clear sky) and maximum (100% cloud cover), and finally the wind speed. In other words, the observed < C F > may be used as a proxy for the actual cloud factor in linearly interpolating between clear sky and overcast LUTs as a function of < C F > . The prospects are promising for ocean biogeochemistry applications.