Satellite Radiation Products for Ocean Biology and Biogeochemistry: Needs, State-of-the-Art, Gaps, Development Priorities, and Opportunities

lidars with ocean proﬁling capability offer the best hope for improving our knowledge of sub-surface ﬁelds. To maximize temporal coverage, space agencies should consider placing ocean-color instruments in L1 orbit, where the sunlit part of the Earth can be frequently observed.

Knowing the spatial and temporal distribution of the underwater light field, i.e., the spectral and angular structure of the radiant intensity at any point in the water column, is essential to understanding the biogeochemical processes that control the composition and evolution of aquatic ecosystems and their impact on climate and reaction to climate change. At present, only a few properties are reliably retrieved from space, either directly or via water-leaving radiance. Existing satellite products are limited to planar photosynthetically available radiation (PAR) and ultraviolet (UV) irradiance above the surface and diffuse attenuation coefficient. Examples of operational products are provided, and their advantages and drawbacks are examined. The usefulness and convenience of these products notwithstanding, there is a need, as expressed by the user community, for other products, i.e., sub-surface planar and scalar fluxes, average cosine, spectral fluxes (UV to visible), diurnal fluxes, absorbed fraction of PAR by live algae (APAR), surface albedo, vertical attenuation, and heating rate, and for associating uncertainties to any product on a pixel-by-pixel basis. Methodologies to obtain the new products are qualitatively discussed in view of most recent scientific knowledge and current and future satellite missions, and specific algorithms are presented for some new products, namely sub-surface fluxes and average cosine. A strategy and roadmap (short, medium, and long term) for usage and development priorities is provided, taking into account needs and readiness level. Combining observations from satellites overpassing at different times and geostationary satellites should be pursued to improve the quality of daily-integrated radiation fields, and products should be generated without gaps to provide boundary conditions for general circulation and biogeochemical models. Examples of new products, i.e., daily scalar PAR below the surface, daily average cosine for PAR, and sub-surface spectral scalar fluxes are presented. A procedure to estimate algorithm uncertainties in the total uncertainty budget for above-surface daily PAR, based on radiative simulations for expected situations, is described. In the future, space-borne lidars with ocean profiling capability offer the best hope for improving our knowledge of sub-surface fields. To maximize temporal coverage, space agencies should consider placing ocean-color instruments in L1 orbit, where the sunlit part of the Earth can be frequently observed.

INTRODUCTION
From the point of view of biology and chemistry, solar radiation in the photosynthetically active range (roughly 400-700 nm), referred to as PAR, controls the growth of aquatic plants (e.g., Ryther, 1956;Platt et al., 1977;Kirk, 1994;Falkowski and Raven, 1997). It ultimately regulates the composition and dynamics of marine ecosystems. Solar radiation in the ultraviolet (UV) (280-400 nm), by damaging cellular constituents, may stress phytoplankton and inhibit their growth (e.g., Cullen and Neale, 1994;Häder et al., 2011). UV light, via photooxidation of colored dissolved organic matter (CDOM), may increase the bioavailability of nutrients (Sulzberger and Durisch-Kaiser, 2009). In the process, absorption by CDOM is reduced, increasing light penetration. Knowing the distribution (spectral, spatial, and temporal) of UV and visible solar radiation in the upper ocean is critical to understanding biogeochemical cycles of carbon, nutrients, and oxygen, and to addressing climate and global change issues, such as the fate of anthropogenic atmospheric carbon dioxide (CO 2 ).
From the point of view of physics, sunlight absorbed by phytoplankton and other water constituents (CDOM, mineral particles, etc.) heats the upper ocean and distributes heat horizontally and vertically, affecting mixed-layer dynamics and oceanic circulation (e.g., Nakamoto et al., 2000Nakamoto et al., , 2001Ballabrera-Poy et al., 2007). These changes in turn influence atmospheric temperature and circulation, with remote effects (Miller et al., 2003;Shell et al., 2003). Solar radiation diffusely reflected by the ocean also affects the outgoing radiative flux from the planet (planetary albedo), with climate consequences (Frouin and Iacobellis, 2002). In order to make predictions for future conditions, we need to get some idea of how the phytoplankton concentrations and optical properties will evolve with changing conditions. Many processes and feedbacks in which solar radiation absorption plays a role are involved and difficult to untangle, and a large fraction of the uncertainties in projections of future climate is associated with physicalbiological interactions (Friedlingstein et al., 2006). This article reviews operational satellite radiation products, the user needs and gaps, and it provides a scientific roadmap for the use of, and priorities for improving existing products and developing new products, i.e., for closing the gaps in ocean biology and biogeochemistry, including studies of biologicalphysical interactions and feedbacks. This roadmap emerged from the presentations (oral and poster) and discussions during the Color and Light in the Ocean from Earth Observation (CLEO) workshop at ESRIN, Frascati, Italy on 6-8 September 2016. The following questions are addressed: (1) Do existing shortwave downward flux products meet the requirements of the dynamics and bio-geochemical communities? What can be done to serve better the needs of the user community in general and the modeling community in particular? (2) What additional products should be added to the processing streams to increase their usefulness? What should be the characteristics of these products in terms of temporal, spatial, and spectral resolution, spectral range, and accuracy? (3) What are the needs in terms of harmonization between sensors, methodologies, ancillary data and radiative transfer tools?

Current Products
The underwater light field is defined at any point in space by the spectral radiance (W/m 2 /sr/nm) from all directions. Useful properties can be derived from radiance, i.e., planar and scalar irradiance, average cosine, reflectance, and vertical attenuation coefficient (see, e.g., Mobley, 1994or Kirk, 1994. Only a few of these properties are presently inferred reliably and operationally from space, namely daily above-surface downward planar irradiance integrated from 400 to 700 nm (known as "daily PAR product"), above-surface planar spectral UV irradiance at noon, spectral reflectance of the water body, and diffuse attenuation coefficient at 490 nm (derived from water reflectance). These products are generated for each ocean color mission individually (note that UV products are not available in standard ocean-color missions).
Examples of level-3 daily, weekly, and monthly above-surface PAR products at 9 km resolution from MODIS-Aqua data are displayed in Figure 1. Dates are March 22, March 22-29, and March 1-31, 2010, respectively. The NASA Ocean Color Biology Group (OBPG) in Greenbelt, Maryland generates and archives these products operationally. Typical uncertainty (RMS) is ±6.5 (19%), ±4.2 (12%), and ±2.6 Em −2 d −1 (7%) for daily, weekly, and monthly estimates . The daily maps exhibit missing values, especially at low latitudes, due to the limited spatial coverage of the instruments, but the weekly and a fortiori monthly maps are completely filled, except at latitudes where the Sun zenith angle at the time of satellite overpass is above 75 degrees, since the data is discarded in the OBPG oceancolor processing line. The weekly and monthly PAR fields exhibit similar patterns, but the monthly product is smoother (lower variability at small scales), which is expected when the averaging period is longer.
The corresponding maps of level-3 daily, weekly, and monthly diffuse attenuation coefficient (K d ) at 490 nm, also produced by the OBPG, are displayed in Figure 2. RMS uncertainty on logtransformed instantaneous estimates is about ±0.1 (e.g., Morel et al., 2007). Since K d is only retrieved in clear sky conditions, the daily map has many gaps (only about 10-15 % of the observed pixels typically pass through the strict glint and cloud filters), and the weekly product show many areas with no information. On a monthly time scale, information is still missing at low latitudes. The spatial gaps limit considerably the utility of the K d retrievals for propagating light below the surface. In the open ocean, global coverage every 3-5 days is necessary to resolve variability associated with seasonal biological phenomena such as phytoplankton blooms. In coastal waters, wind forcing create "events" (e.g., upwelling) that occur every 2-10 days, and 1day coverage is the requirement for resolving the event time scale. Figure 3 provides examples of level 3 daily, weekly, and monthly maps of noon surface UV irradiance at 1 • resolution and 324 nm from OMI-Aura. Irradiance at 305, 310, and 380 nm as well as erythemally weighted daily dose and erythemal dose rate are also available. Dates are the same as those of Figure 1. The data are routinely processed at the OMI Science Investigatorled Processing System (SIPS) Facility in Greenbelt, Maryland, and are archived at the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC). Overall uncertainty of UV irradiance estimates ranges from ±5 to over ±30%, depending on atmospheric conditions and geolocation (e.g., Arola et al., 2009). The UV irradiance and PAR fields have similar spatial coverage (Figures 1, 3), with missing values due to instrument swath on a daily time scale, except that UV irradiance estimates are obtained at high latitudes. Variability patterns are also similar, since UV irradiance variability is also governed by Sun zenith angle and cloudiness, although ozone absorption plays a bigger role in modulating the surface values (but gradients of total ozone content remain mostly latitudinal). Spatial resolution is coarser than for the MODIS products (OMI sensor footprint is 13 × 24 km 2 ), , and no daily-averaged values (only daily values at noon or overpass time) are generated.
The situation regarding spatial coverage is summarized in Figure 4, which displays the percentage of the ocean surface covered by PAR, K d , and UV irradiance products on daily, weekly, and monthly time scales (imagery of Figures 1-3). In the equatorial region, percent coverage is 65, 80, and 5% for daily PAR, UV irradiance at 324 nm, and K d , respectively. It is increased to 100% for weekly and monthly PAR and UV irradiance, and to 30 and 70% for weekly and monthly K d . In middle latitude regions, daily coverage is almost 100% for PAR, about 80% for UV irradiance, and 15-20% for K d . Weekly and monthly coverage is 100% for PAR and UV irradiance, and reaches 70-75 and 100% in sub-tropical regions for K d . At high latitudes (>70 • ), monthly coverage is <40% for PAR and K d . This lack of coverage is limiting in view of the large productivity of high latitude marine ecosystems (Southern ocean, Arctic ocean). Furthermore, monthly K d products may only contain estimates during a few days, and therefore may not represent accurately actual monthly values in dynamic regions (IOCCG, 2015).
Multiple satellites can improve the daily ocean coverage, especially for K d (can only be retrieved in clear sky conditions). For example, three instruments of MODIS type, flying in a constellation on satellite orbits differing by the mean anomaly (angular distance from pericenter), would increase the daily spatial coverage of water reflectance, therefore K d , from 15 to 25% over 1 day and from 40 to 60% over 4 days (Gregg et al., 1998). Figure 5 shows, for March 22, 2010, the increase in daily ocean coverage obtained by combining estimates from MODIS-Aqua (overpass at 13:30 local time) and -Terra (overpass at 10:30 local time) instead of using MODIS-Aqua only. For PAR, the increase is from 60 to 80% in equatorial regions, and complete coverage is reached at sub-tropical latitudes. For K d , the ocean coverage is more than doubled at most latitudes (e.g., 35-50%, instead of 15-20% in the sub-tropics). Combining PAR estimates from instruments orbiting at different times not only increases spatial coverage, but perhaps more importantly, also takes into account cloud diurnal variability, yielding more accurate estimates.
Satellite instruments in geostationary orbit, by observing the same target multiple times during the day (e.g., every 30 min for GOCI onboard COMS and 10 min for AHI onboard Hiwamari-8) offer an efficient way to account for diurnal changes in cloudiness in daily PAR products. Figure 6 displays daily PAR imagery obtained with GOCI data acquired on April 5, 2011 at 3:16 GMT and at 00:16, 01:16, 02:16, 0.3:16, 04:16, 05:16, 06:16, and 07:16 GMT. In clear-sky regions (north of Japan), the values are close using one or 8 observations, which is expected since the governing parameter is the Sun zenith angle. In cloudy regions (South of Japan), the PAR spatial field is smoother and the range of values smaller since cloudiness changes are accounted for in the daily average. The lowest value is about 12 Em −2 d −1 with 8 observations instead of 5 Em −2 d −1 with one observation. The two types of estimates compare well, with a bias of 0.11 Em −2 d −1 (0.2%), i.e., slightly higher values using one observation, and a root-mean-squared (RMS) difference of 5.92 Em −2 d −1 (13.5%), largely due to differences in cloudy situations. Figure 7 displays an example of AHI daily PAR product (July 20, 2011) and the corresponding MODIS-Aqua product. Observations every 10 min were used to estimate daily PAR from AHI data. The patterns of variability are similar in both products, but as for GOCI, the AHI PAR imagery is smoother and contains less extreme values. The AHI product, unlike the MODIS product, does not exhibit spatial gaps. In terms of comparison statistics, the AHI values are lower by 1.38 Em −2 d −1 (4.8%) on average, and the RMS difference is 6.46 Em −2 d −1 (22.7%).
In summary, existing satellite products generally do not cover the global open oceans (e.g., retrievals limited to Sun zenith angles < 75 • ), except for UV irradiance, and they do not provide information below sea ice, where significant blooms may develop (e.g., Arrigo et al., 2012). In addition, cloud diurnal variability FIGURE 5 | Percentage of ocean coverage by daily PAR and K d products when using MODIS-Aqua data only (red curves) and combining MODIS-Aqua and -Terra data (black curves). Date is March 22, 2010. Coverage is almost total (>80%) in low to middle latitude regions for PAR (left) and at least doubled at most latitudes for K d (right).
FIGURE 6 | Daily above-surface PAR distribution for April 5, 2011 obtained from 8 hourly GOCI observations during the day (right) and a single GOCI observation at 03:16 GMT (left). In the stormy region South of Japan, the range of PAR values is smaller when using 8 observations due to diurnal changes in cloudiness. Reproduced with permission from Springer Nature (After . Frontiers in Marine Science | www.frontiersin.org is not accounted for in daily PAR calculations when using polar orbiting satellites. Geostationary satellites observing frequently during the day account properly for cloudiness changes, but the drawback is a decreased spatial resolution at high latitudes and potentially large uncertainties for slanted viewing geometries. Propagation of surface radiation to depth currently assumes that the ocean is homogenous, neglecting potentially important effects of stratification on the absorption of solar radiation. In sum, our view of the underwater light field from space is limited. Nevertheless, the operational radiation products have been used to address a variety of topics related to aquatic photosynthesis, for example biosphere productivity during an El Niño transition (Behrenfeld et al., 2001), phytoplankton class-specific productivity (Uitz et al., 2010), chlorophyll 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), and inter-comparison of productivity algorithms (Carr et al., 2006;Lee et al., 2015). They have also been used to check the stability of CERES measurements (Loeb et al., 2006).

Users Needs
User needs vary widely in terms of products, spectral, spatial, and temporal resolution, and acceptable uncertainties, depending on the scientific or societal subject of interest. Applications requiring knowledge of radiomeric quantities and apparent properties (radiance, irradiance, average cosine, attenuation coefficients) are multiple and diverse, including phytoplankton phenology, carbon inventory, heat budget and ocean dynamics, fisheries and ecosystem management, toxic algal blooms, and eutrophication (see National Research Council, 2011 for a comprehensive list). Observational and uncertainty requirements (satellite products) generally range from 1 h to 1 day, 0.1 to 50 km, and ±5 to ±20% for PAR and 0.1 to 10 km, 1 to 7 days, and ±10 to ±25% for spectral K d (Malenovsky and Schaepman, 2011), but higher resolution may be needed in some cases, for example rapidly changing phenomena occurring in small water bodies. For applications that need analyzing long-term records (e.g., associated with climate), the products need to be sensor independent, consistent, and continuous across satellite missions.
The satellite products should be defined unambiguously and completely, they should be easily accessible, and they should have associated ATBDs with detailed protocols including description of all ancillary data used and their sources. For example, defining a PAR product merely as downward quantum flux at the surface in the 400-700 nm spectral range is insufficient. One needs to precise whether the flux is just above or just below the surface, whether it is instantaneous or time-averaged (e.g., over 24 h for "daily PAR"), and to indicate spatial resolution. Advantages and limitations should be specified (e.g., product not valid at Sun zenith angles above 75 • , or over sea ice, or in the presence of Sun glint), and uncertainties assessed, preferentially provided on a pixel-by-pixel basis, to make sure that observed changes are interpreted correctly. Computer codes used to derive the products should be available to users with proper documentation, as well as standardized processing tools (e.g., open-source toolboxes).
A survey about satellite PAR observations was conducted in 2015 by Plymouth Marine Laboratory (PML), requesting feedback from the user community on adequacy of available products, including importance, usage, and accuracy, and additional features one would like to see. Figure 8 summarizes the results about desired PAR attributes and new products and acceptable PAR uncertainties. Most respondents answered it was very or extremely important to have uncertainties associated with PAR products. About 50% of the respondents indicated that ±10-25% uncertainty was acceptable, and about 40% wanted uncertainty better than ±10%. A substantial majority of FIGURE 7 | Daily above-surface PAR for July 20, 2011 obtained from AHI (observations every 10 min) and MODIS-A data at 9 km resolution (left and right, respectively). Patterns of spatial variability are similar for the two sensors, but the AHI PAR field is smoother, contains less extreme values, and does not exhibit spatial gaps (Courtesy of H. Murakami, Japanese Aerospace Exploration Agency). respondents ranked PAR below the surface and the fraction of PAR absorbed by phytoplankton about equally as the top two new products to generate. Based on this survey and in view of the extensive list of research and societal applications, current and potential, a list of required (new) products was compiled, including products that may be challenging to generate from space: • Sub-surface planar and scalar irradiance (as opposed to above the surface). • Fraction of PAR absorbed by phytoplankton (APAR).
• Diffuse fraction of total irradiance (average cosine of light field just below the surface). • Spectral planar and scalar irradiance.
• Spectral diffuse attenuation coefficient for downward irradiance. • Surface albedo (ratio of planar upward irradiance to downward planar irradiance just above the ocean surface). • UV-A, UV-B scalar irradiance (with photon and energy units).
• Products without gaps (in space/time) to provide boundary conditions to models. • Upper-ocean heating profile.
• Diurnal distribution of PAR and its attenuation.

Gap Analysis
Some new products (such as below surface planar and scalar PAR) can be easily implemented while others require development (e.g., APAR, vertical attenuation of PAR). For some products (e.g., under-ice light fields), the readiness level, in terms of methodology, is low. The state-of-the-art, however, is such that the strategy to obtain the new products described above is known and summarized below.
Sub-surface planar PAR and scalar PAR depend essentially on the sunlight transmission across the air-water interface, therefore the angular distribution of radiance incident at the surface and surface roughness (Mobley and Boss, 2012). The 24 h-averaged quantities, as well as the average cosine for total light, can be parameterized as a function of latitude and daily cloud factor (i.e., the ratio of actual PAR and clear-sky PAR) and wind speed. This may require look-up tables for clear sky and overcast quantities. Details about procedures are provided in section Examples of New Products, where examples of subsurface products are presented. Note that in the OBPG approach to estimating above-surface daily planar PAR, it is actually easier (more direct) to compute the downward flux below the surface. This "penetrative" flux is obtained by subtracting from the incident extraterrestrial solar irradiance the reflected flux and the flux absorbed by the surface/atmosphere system. The "penetrative" flux is then corrected by 1/(1−A s ), where A s is the surface albedo, to yield the incident flux onto the surface (e.g., . This second step introduces uncertainty, but allows for extensive evaluation at PAR measuring sites, an activity that cannot be performed easily for sub-surface fluxes (lack of data and difficulty to measure sub-surface fluxes accurately).
The above variables, integrated over the PAR spectral range, can be calculated without difficulty for the spectral bands of the ocean color sensors (i.e., visible to near infrared), in which cloud absorption is negligible. Providing the information at regular spectral intervals (e.g., every 5 or 10 nm), a requirement of some primary production models (e.g., Sathyendranath et al., 1989;Antoine et al., 1996), is straightforward since cloud optical properties (extinction coefficient, asymmetry factor) are similar, i.e., cloud albedo can be assumed constant in the entire spectral range, and the coupling between molecules and cloud droplets/crystals is relatively small (i.e., fairly unique relation between cloud transmittance at different wavelengths). In spectral regions of strong gaseous absorption, however, uncertainties may be introduced due to the coupling between absorption and scattering processes (depends on the unknown vertical distribution of the absorbers and scatterers). Extending the calculations to the ultraviolet (e.g., UV-A, UV-B) from measurements in the visible is more complicated, but definitely feasible. The complication is not due to cloud optical properties (they remain similar to those in the visible), but to the coupling between molecules and cloud droplets/crystals, which is effective. In other words, the relation between cloud transmittance (in the presence of molecules and aerosols) in the ultraviolet and visible depends on the type of clouds and their location in the vertical. This indicates that de-coupling the clear atmosphere from clouds, as it is done in the OBPG PAR algorithm, may introduce significant errors in the estimation of ultraviolet irradiance. Suitable modeling of the relation between atmospheric transmittance at visible and ultraviolet wavelengths is therefore required, which can be accomplished via radiation transfer calculations to various levels of accuracy depending on user needs (may require additional information on clouds).
From the spectral downward irradiance at the surface, minimally affected by photons reflected by the surface and backscattered by the water body (spherical albedo of the atmosphere is small, i.e., about 0.15 in the visible) and the spectral upward irradiance at the surface, one may compute the spectral surface albedo. The computation, however, can only be done in clear-sky conditions since water optical properties are necessary but generally not retrieved in cloudy conditions from ocean-color sensors. Estimating the spectral upward irradiance at the surface requires retrieving spectral water reflectance (i.e., the signal backscattered by the water body), which is routinely achieved by satellite project offices, and modeling the bidirectional reflectance function of the surface (e.g., using Cox and Munk, 1954) and the water body (e.g., using Gentili, 1996 or Park andRuddick, 2005).
Although daily averaged quantities are often used (therefore required) in applications, instantaneous quantities (i.e., determined at time of satellite overpass) can be easily provided. In fact instantaneous products, unlike 24 h-averaged products, do not require assumptions about diurnal changes in atmospheric and oceanic properties. Diurnal variability can be described well using sensors onboard geostationary satellites, such as GOCI and AHI, see section Current Products, all the more as these sensors have ocean-color capabilities. One limitation, however, is the reduced spatial resolution at high latitudes, and managing data from different instruments, which may be operated by different space agencies, to achieve global coverage. Observations from the Earth-Sun the Lagrangian-1 (L1) point, 1.5 million kilometers from Earth, such as those made by the EPIC camera onboard DISCOVR (1-2 h temporal resolution, 21 km spatial resolution), provide a suitable alternative to several sensors on geostationary orbit. Because of the high orbit, spatial resolution at high latitudes is much less an issue. Diurnal information may also be obtained from a constellation of Sun-synchronous instruments with local overpass times spread during the day, for example MODIS-Terra at 10:30 am, SeaWiFS-SeaStar observing at noon, and MODIS-Aqua observing at 1:30 pm. From a unique instrument in Sun-synchronous orbit, ancillary data about variability of clouds and aerosols is necessary, for example Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) products, available at a 1/2 × 2/3 degree every hour for the day of the satellite observation.
Propagating fluxes vertically below the surface requires knowledge of the vertical profile of diffuse attenuation coefficient, K d . From space, this can only be achieved in clear sky conditions with passive optical sensors, since K d is deduced from water reflectance (empirical algorithms) or inherent optical properties (e.g., IOCCG, 2006). The K d estimates are actually weighted averages over one optical depth (from which most of the photons from the water body reaching the satellite sensor originate); no vertical information is obtained. In first approximation, one may propagate light below the first optical depth (shallower than the depth of the euphotic zone) by assuming no vertical variation in spectral diffuse attenuation. However, this is generally inaccurate, as many oceanic regions (e.g., oligotrophic provinces) exhibit maximum chlorophyll concentration well below the first optical depth. One has therefore to rely on statistical relations between concentrations of oceanic constituents at the surface and below (e.g., Morel and Berthon, 1989) to estimate the K d depth profile, or to use outputs of predictive coupled physical-biogeochemical numerical models. In the future, with the advent of space-borne polarization lidars such as CALIOP onboard CALIPSO or the Aerosol/Cloud/Ecosystems (ACE) lidar (being designed), one will be able to profile K d in both clear and cloudy conditions, day and night, up to 3 optical depths in the green (532 nm) at a vertical resolution of 3 to 30 m (e.g., Lu et al., 2014;Behrenfeld et al., 2016). The sub-surface fluxes and the vertical profile of diffuse attenuation coefficient give access to fluxes at the bottom (important for studies of shallow coastal ecosystems), average fluxes in the mixed layer (requires knowledge of mixed-layer depth, e.g., from ocean circulation models), and the upper-ocean heating rate profile. Knowing mixed-layer fluxes and vertical heat distribution is especially useful to characterize the role of solar penetration and biological-physical interactions on ocean circulation and climate (Olhmann et al., 1996;Shell et al., 2003; see section Introduction).
In a homogeneous ocean, APAR depends on the ratio of the spectral absorption coefficient by live phytoplankton, a ph , and total absorption (water, yellow substances, non-algal particles), a tot , in the PAR spectral range and the spectral planar irradiance just below the surface. In a vertically heterogeneous ocean, the vertical distribution of those quantities plays a role, as well as the vertical distribution of the diffuse attenuation coefficient for downward irradiance, K d . Computing APAR from space, therefore, requires estimates of spectral planar irradiance below the surface, and vertical profiles of a ph , a tot , and K d . Subsurface spectral irradiance and its vertical attenuation can be obtained as discussed above, and absorption coefficients using various techniques (e.g., IOCCG, 2006). Some of these variables are difficult to retrieve with good accuracy, in particular a ph (requires partitioning a tot into its components), and vertical information in the euphotic zone is generally not directly available (except from future space-borne lidars, see above). Consequently, uncertainties on APAR computations based on satellite estimates of individual variables may be large. Since APAR strongly depends on the spectral ratio of sub-surface reflectance, R(0 − ) and pure seawater reflectance, R w (0 − ), one may envision approximating APAR by a linear combination of R(0 − )/R w (0 − ) in the PAR spectral range (Frouin et al., 2014). Since APAR is expressed linearly in water reflectance, the method is applicable to average values of water reflectance (e.g., spatially averaged), which may reduce the impact of water reflectance noise on the APAR estimate.
Primary production under sea ice is considerable, as evidenced from in situ measurements of phytoplankton concentrations and suggested by numerical model simulations. In the Arctic Ocean, it may constitute more than 30% of the total production (Popova et al., 2010). As a consequence, the seasonal cycle is shifted, with maximum primary production occurring in July, not in August-September (case of open waters). There is a need, for primary production modeling and studies of underice phytoplankton blooms, to determine light fields under the ice (Laliberté et al., 2016). This is quite difficult to realize from space, since knowledge of sunlight transmission though sea ice is required, and this parameter is quite variable depending on ice type and thickness, and the presence of snow and melt ponds. Transmission models are now becoming available (Arndt and Nicolaus, 2014), and they can make use of sea-ice thickness and age, snow depth, and melt pond fraction observations from microwave and optical sensors (SIRAL on Cryosat-2, AMSR-2 on CGOM-W1, ATLAS on IceSta-2, MODIS on Terra and Aqua). One issue is separating the contribution of clouds to the planetary albedo, in order to access the surface albedo. Algorithms are not mature for generating operationally shortwave fluxes under sea ice.
The possible methodologies, difficulties, opportunities, and readiness level for developing and creating the new products were discussed above. The following recommendations regarding these products are made: 1) Spectral fields should be provided at the sensor resolution with protocols (and codes) describing how to interpolate and extrapolate for obtaining other spectral distributions (e.g., 1 nm irradiance field from a multi-spectral sensor). 2) Vertical propagation of products requires an appropriate attenuation coefficient from which other products can be derived (e.g., euphotic depth, mixed-layer depth, isolume depth). Clear guidelines on how to produce the derived products using the attenuation should be provided. 3) Horizontal/temporal gap filling is necessary for certain applications (e.g., ecological forecasting). This can be done using merged products across sensors and/or interpolation schemes (using known de-correlation scales or models). Many techniques are available (e.g., Pottier et al., 2006;Alvera-Azcarate et al., 2007;Krasnopolsk et al., 2016). 4) Products should have associated uncertainties that are consistent with those obtained when validating estimates against in situ measurements, taking into account uncertainties in the in situ data. This requires a calibration/validation program. The product protocol should provide a description of how the uncertainty was derived. It is desirable to provide a per-pixel uncertainty. It is recognized that the level of effort to obtain a very accurate uncertainty estimate can be very large and therefore some trade-offs may need to be done, in consultation with user requirements. 5) Data access should be tailored to users need. For example, modelers will use multithreading with simultaneous access to associated error fields. In contrast, the Earth observation community will, typically, want to access data using FTP. Most users do not care about the satellite mission from which a product was derived, but rather care about the products being continuous in time and consistent across missions. 6) Cross-agency efforts should be made to homogenize their respective products so it is easy for users to use these products (e.g., the definition of a PAR product should be the same). For climate relevant products, it is critical to merge them (and debias) across missions so that models to not experience secular jumps as they assimilate such data.

Radiative Transfer Code
To develop new shortwave radiation products from satellite data (such as those discussed above) and assess accuracies, one needs to simulate the TOA radiance measured by a given sensor and the variable to retrieve (e.g., planar or scalar spectral irradiance below the surface, average cosine). For this, we use the Speed up Monte-Carlo Advanced Radiative Transfer using GPU (SMART-G) radiative transfer code (Ramon et al., 2017). This code, based on the Monte-Carlo method, is fast and massively parallel. It computes the complete light field (i.e., radiance, including polarization) in the ocean and atmosphere. The code simulates the transfer/propagation of solar radiation in a 1-dimensional coupled ocean-atmosphere system with a wavy interface. It accounts for absorption and scattering by molecules, aerosols, and hydrosols, and Fresnel reflection/refraction at the interface. Polarization properties of the various atmospheric and oceanic constituents and the surface are explicitly considered. Inelastic processes (i.e., Raman scattering, fluorescence) are omitted in the current version. The ocean can be infinitely deep or bounded by a reflective bottom at finite depth. The computations are made in either plane-parallel or spherical geometry. The four components of the Stokes vector can be obtained at any wavelength of the solar spectrum and any level of the coupled ocean-atmosphere system.
Gaseous absorption is treated either by correlated kdistribution (Kato et al., 1999) or an equivalent like REPTRAN (Gasteiger et al., 2014;Emde et al., 2016). In the multispectral mode, each photon is assigned a wavelength. All optical properties of the medium are pre-calculated for these wavelengths. The spectrum is computed in one pass but is under the influence of Monte-Carlo noise. For higher spectral resolution computations without spectral noise, for example in order to handle line-by-line (LBL) gas absorption, the ALIS method (Emde et al., 2011) is also implemented in SMART-G. This allows for the calculation of spectra by tracing the photon paths once for all wavelengths (by absorption/scattering decoupling).

Simulations
In the Monte-Carlo code, photon packets are carrying planar irradiance (Wm −2 ) perpendicular to their direction of propagation. Therefore downward or upward spectral planar irradiances are obtained by performing the weighted sum of irradiances crossing a unit area of horizontal surface located above or below the air-sea interface. For spherical irradiances, each photon is assigned an additional weight of 1/cos(θ i ) where θ i is the photon's zenith angle arriving on the detector. The runs were optimized for the calculation of spectral daily fluxes above and below the ocean surface. For 1 day, one run was executed by injecting a large number of photons whose wavelength was chosen randomly as well as the injection angles at TOA and in the atmosphere depending on the hour of the day. Typical runtime is 1 s to reach ±1% uncertainty for clear sky conditions and for 1 day. It becomes 30 s for a totally overcast situation with a cloud optical thickness of 50.
The aerosol and cloud optical properties are taken from the OPAC database (Hess et al., 1998) and distributed within the libradtran software package (www.libradtran.org). Aerosols are supposed to be spherical. The clouds are supposed to be composed of liquid water droplets with a varying effective radius. Rayleigh optical depth is computed according to Bodhaine et al. (1999). The gaseous absorption is parameterized as a correlated k-distribution with spectral intervals of 10 nm from line-by-line calculations using the Py4Cats code (Schreier and Gimeno Garcia, 2013) and HITRAN 2012(Rothman et al., 2013 absorption parameters for H 2 O and O 2 . Ozone and NO 2 smooth absorption coefficient are taken from Bogumil et al. (2003). The wind-roughened sea surface is modeled as an ensemble of uncorrelated facets with a slope distribution (Cox and Munk, 1954) and a uniform azimuth distribution. The reflection and transmission coefficients are then computed using Fresnel's laws. Ocean bulk optical properties correspond to Case-I waters with a chlorophyll-a concentration of 0.5 g.m −3 . The ocean's bottom is black and is located at a depth of 50 m. The ocean phase function is represented by a Fournier-Forand function with a varying truncation angle.
An example of outputs of the SMART-G code in the context of PAR simulations (both spectral irradiance and reflectance) is given in Figure 9 for various levels in the atmosphere-ocean system (top, just above the surface, just below the surface, and at the black bottom of the water column) and several atmospheric conditions (clear and cloudy situations with different Sun zenith angles). The TOA reflectance in the MERIS bands is also indicated. The various graphs show the importance of Sun zenith angle and cloud optical thickness in controlling the downward PAR above and below the surface. They also illustrate the usefulness of the SMART-G tool for algorithm development.
For the year 2011, and for 14 latitudes ranging from −64.5 to 64.5, the various spectral irradiances listed in the Annex were computed at 11 times regularly distributed throughout the day. The typical spectral TOA reflectance for a MERIS-like instrument measuring at 10:30 local time was also computed. The atmospheric content was changing according to the MERRA-2, hourly, gridded datasets of water vapor and ozone contents, aerosol optical depths of black carbon, dust, organic carbon, sea salt, and sulfates aerosols at 550 nm, and cloud optical thicknesses. Individual clear sky and overcast plane-parallel radiative transfer calculations were then mixed using the Independent Pixel Approximation (IPA) using MERRA-2 cloud cover variable as the mixing value.
The following quantities were then computed: PAR 0− o (t) and PAR 0+ d (t) (to describe diurnal variation of sub-surface fluxes), PAR 0+ d 24h (the main product), PAR 0− o 24h (the key product for primary production and photo-chemical processes), PAR 0− o 24h (λ i ) (spectral scalar flux, as requested by modelers, at a resolution of 10 nm), and µ 0− 0 24h (to characterize the angular structure of the light field). Definitions are provided in the Annex.

Rationale
We define a Cloud Factor (CF) as the deviation from a pure clear sky daily averaged PAR above the surface. It is a "measure" of the influence of clouds on the daily PAR: One typical day of simulations is displayed in Figure 10 for a latitude of 55.5 • N. The Sun zenith angle and the length of the day mainly drive the diurnal cycle of both clear sky PAR and partly cloudy PAR. For that particular day, the influence of cloudiness is important (and somehow stable) reducing the daily PAR by a factor CF 24h = 0.34. The average cosine µ 0− 0 is very stable for cloudy conditions, with a value around 0.82 throughout the day. A clear sky has a variable µ 0− 0 culminating at noon. In both cases the spectral variation of µ 0− 0 , not shown here, is weak. The 24 haveraged value of µ 0− 0 should be representative of the direction of propagation of the largest fraction of the daily PAR. That is why we define the 24 h-averaged µ 0− 0 24h as the ratio of the 24 h-averaged net and scalar PAR (and not the 24 h average of their instantaneous ratios): The normalized spectral PAR, obtained by dividing the spectral PAR by its value at 675 nm, is defined in the same manner as: It gives the spectral shape of the PAR and it is very stable and close to the TOA solar irradiance spectrum in most cases. The spectral shape of the clear sky PAR is slightly influenced by the mean Sun zenith angle because the absorption of ozone in Chappuis bands becomes more and more effective. The spectral dependence of the cloudy PAR is smaller. It increases in the blue part of the spectrum as cloud influence increases.
Following the approach of Mobley and Boss (2012), 24 h-averaged secondary radiative quantities may be obtained from 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. The chlorophyll content FIGURE 9 | Examples of outputs of the SMART-G RTC code for 4 situations. On the left axis is reported the "spectral PAR" in units of light quanta per unit time, unit area and unit wavelength for planar and spherical upwelling ("u") or downwelling irradiances ("d") at several levels in the ocean-atmosphere coupled system: at TOA, just above (0 + ) or below (0 − ) the surface, or at the black bottom (superscript B) of the ocean located here at a depth of 50 m. On the right axis is plotted the TOA spectral reflectance at the same resolution, and at the center of MERIS wavebands. The top two graphs are for a clear-sky situation for two Sun zenith angles (SZA) and the bottom two graphs a liquid water cloud, located between 2 and 4 km with droplets of effective radius = 11 mm whose of optical thickness 10 is added. The clear atmosphere model is the US62 standard atmosphere, with maritime polluted aerosols as described in the OPAC database with an AOT at 550 nm equal to 0.1. The ozone column is 300 DU and the precipitable water quantity is 2 g/cm 2 .
of the water is of minor importance for calculating the scalar and net PAR just below the surface. Figure 11 displays the coefficients to be applied to PAR 0+ d 24 h in order to obtain PAR 0− o 24h , or PAR 0− d 24h , and also shows how µ 0− 0 24h and P AR 0− o 24h (λ i ) vary for various latitudes and wind speeds and for the two extreme cases, i.e., clear and totally overcast. Figure 12 displays the result of the application of the coefficients described above for estimating the scalar PAR below the surface. We processed one full year of global simulated data (2011, see section PAR Simulator) for 3 wind speeds (0, 7, and 15 m.s −1 ) and checked the quality of the regression vs. the "actual" (or prescribed) values, i.e., the values obtained by running the Monte Carlo code with the various input variables (MERRA-2 hourly data, see above). These values provide the reference in calculating algorithm performance statistics. The regression is excellent with a residual bias and a R.M.S. difference of about 0.5 mol.ph m −2 .day −1 . Wind speed or cloudiness do not apparently impact the quality of the regression. Figures 13, 14 depict the normalized spectral PAR and average cosine for the year 2011, for a latitude of 55.5 • N and a wind speed of 7 m.s −1 . For both parameters, the clear sky and totally overcast situations are also plotted. They constitute an envelope wihin which the actual values are included, and the deviation from the clear sky value seems proportional to the actual cloud factor CF 24h . This suggests a method to derive P AR 0− o 24h (λ i )

Equations
We propose to use the observed cloud factor CF obs as a proxy of CF 24h and then linearly interpolate between clear sky and overcast look-up tables as a function of CF obs . We have: The error will be large when cloudiness changes a lot during the day and thus we may suspect that CF obs deviates substantially from CF 24h and when the values of the parameters in clear or totally cloudy conditions are substantially different.

Look-Up Tables for Clear Sky and Overcast Situations Models
The clear sky model is based on the AFGL US 62 standard atmosphere with a surface pressure of 1012.15 hPa, an O 3 vertically integrated content of 300 DU, and a H 2 O vertically integrated content of 2 g.cm −2 . The aerosol model is the maritime clean model from the OPAC database with an AOT of 0.1 at 550 nm. The air-sea interface is a wind-roughened surface. The ocean bulk optical properties correspond to Case-I waters with a chlorophyll a concentration of 0.5 g.m −3 . The ocean bottom is black and is located at a depth of 50 m. For the totally overcast model we added a permanent cloud layer between 5 and 10 km consisting of water droplets with r eff = 11 µm and a cloud optical thickness of 50. The aerosols, hydrosols and cloud phase matrices are computed at 550 nm and are assumed spectrally invariant between 400 and 700 nm.

Computations
For latitudes between −90 • and 90 • by step of 10 • , for every 30 days along the year, and for 3 wind speeds: 0, 7 and 15 m.s −1 , we computed the following quantities:

First Results
For the MERIS sensor, first examples of the new parameters (daily and monthly global products for May 15 and May 1-31, 2011) are displayed in Figures 15-17. Data at high latitudes were masked using ESA CCI Sea Ice Concentration products v2.0 (cci.esa.int). The scalar PAR below the surface (Figure 15) follows the planar PAR above the surface (not shown here), but the values are somewhat higher, as expected. The spatial coverage of the daily MERIS product is less than for the MODIS products (Figure 1), because of the narrower swath of MERIS and the glitter mask (more glint in the MERIS imagery). The average cosine product (Figure 16 smooth the product, which becomes a simple function of latitude. The normalized spectral PAR at 405 nm (Figure 17), like the mean cosine, exhibits a contrast between clear and cloudy sky situations that is increasing toward the high latitudes. However on a monthly time scale the latitudinal gradient is very weak and the dispersion of the product is also very weak with a mean of 0.63 and a standard deviation of 0.024 over the global ocean.
Algorithm Uncertainties for PAR 0+ d 24h Associating uncertainties to the satellite radiation products, preferentially on a pixel-by-pixel basis, is obligatory to quantify their quality. This is important to ensure that variability or trends detected in scientific analyses are of geophysical nature, i.e., that the data are interpreted properly in view of their strengths and limitations, and to merge different data sets. This is also essential for data assimilation, a primary application of the products (large uncertainties will have little impact on model runs, small uncertainties will constrain the model to behave like the data). Expressing uncertainties requires modeling the measurement, identifying all possible error sources (e.g., noise in the input variables, imperfect/incomplete mathematical model), and determining the combined uncertainty, as described in JGCM-100 (2008) and subsequent publications. In the following, algorithm uncertainties associated with PAR 0+ d 24h are considered, i.e., those due to model approximations and parameter errors (e.g., decoupling effects of clouds and clear atmosphere, neglecting diurnal variability of clouds, using aerosol climatology) assuming that the input variables (TOA reflectance at wavelengths in the PAR spectral range) are known perfectly. A procedure is provided to estimate and provide, for each pixel of a product, this uncertainty component of the total uncertainty budget, which is expected to dominate. The uncertainty characterization has been done using an extended simulation dataset covering the 2003-2012 time period still using 1 hourly MERRA-2 input data. The large number of data points allows one to sample well the atmospheric variability and in particular many variations of daytime nebulosity, for all latitudes.
FIGURE 12 | Estimated scalar PAR below the surface from the planar PAR above surface and the clear sky "transmission factors" shown in Figure 11 for all the days of 2011 and 3 wind speeds. The color of the symbols is a function of the cloud factor. Figure 18 displays the result of the uncertainty analysis for the whole dataset. The bias and the standard deviation of the daily PAR estimates are plotted as a function of the clearsky PAR (which depends itself mainly on latitude and date) and the actual cloud factor for that day. The Monte Carlo calculations, assumed accurate, provide the reference. This is justified in view of the bias and standard deviation values. The bias exhibits a slight dependence on the clear sky PAR, suggesting an overestimation reaching 2.5 E.m −2 .day −1 (∼4%) for the maximum clear sky PAR values, and a slight underestimation 1 E.m −2 .day −1 (∼6.5%) when the clear sky PAR reaches a low value of 15 E.m −2 .day −1 . When looking at the dependence upon CF 24h , a slight overestimation exists between 1 and 1.5 E.m −2 .day −1 . This overestimation is independent of the cloudiness of the day, with the exception of totally clear days ( CF 24h = 1), for which the bias drops to 0.5 E.m −2 .day −1 . The bias is also quite small for totally overcast situations. A bias correction can be considered, based on the clear-sky PAR value.
The standard deviation (SD) is peaked toward intermediate cloud factors, where the risk of deviation of cloudiness at the time of the satellite measurement and the mean cloudiness of the day is maximum. When the clear sky PAR is large and CF 24h is about 0.5, SD reaches 8 E.m −2 .day −1 (∼11%). But SD drops to 1 E.m −2 .day −1 (1.5%) for clear sky situations ( CF 24h = 1). SD is drastically reduced when dealing with monthly PAR estimates. Whatever the cloudiness, SD is lower than 2 E.m −2 .day −1 . The main feature is that SD seems to be proportional to the clear-sky PAR value, and thus we can consider associating an uncertainty to each pixel of the product from an estimate of the clear-sky value, in addition to the cloud factor. This model uncertainty component could be extended to situations with multiple satellite measurements per day, as it is often the case for high latitudes with polar orbiting sensors like MERIS or VIIRS. In that case we anticipate a reduction of the standard deviation of the daily PAR product. For a complete per-pixel uncertainty budget, the uncertainty associated with TOA reflectance noise (radiometric and due to vicarious calibration) should be included, which may require evaluating the sensitivity of the daily PAR to the TOA reflectance and the covariance of the input reflectance in the various spectral bands (since the measurements are correlated).

SUMMARY AND RECOMMENDATIONS
Studying and understanding the chemical, physical, geological, and biological processes that govern the composition of the marine environment requires knowledge of the underwater light field. Ideally, one wants to figure out and monitor the spectral and angular structure of the radiant intensity at any point in the water body. From space, by means of remote sensing, only limited information about the radiative properties of a water body can be obtained, but the advantage is global and repetitive coverage. Existing satellite products are restricted to planar PAR and UV irradiance above the surface, and diffuse attenuation coefficient (average in the upper layer). These products, despite their drawbacks (e.g., no information at high Sun zenith angles and diurnal variability poorly described for PAR, no retrieval in cloudy conditions for the diffuse attenuation coefficient), have been useful to many studies of aquatic photosynthesis, heat budget, and chemical effects of light. There is a need, however, for other products, i.e., sub-surface planar and scalar fluxes, average cosine, spectral fluxes (UV to visible), diurnal fluxes, APAR, surface albedo, vertical attenuation, and heating rate, and for associating uncertainties to any product on a pixel-bypixel basis. Possible approaches and methodologies to generate these new products, in view of state-of-the-art knowledge and current and planned satellite missions, were discussed, including difficulties, assumptions, and readiness level. A strategy and FIGURE 14 | Mean cosine PAR for all the days of 2011 (dots), for a particular latitude and wind speed, as well as predicted values for clear sky (solid line) and totally overcast situations (dashed lines). The color of the dots is a function of the cloud factor. a roadmap, with development priorities and opportunities to obtain the new products, could be established. Examples of new products, i.e., daily scalar PAR below the surface, daily average cosine for PAR, and spectral scalar fluxes below the surface, and their algorithms, were presented. A statistical way to estimate uncertainties for each pixel, based on radiative transfer simulations for expected clear and cloudy situations, was proposed.
In the short term, the focus should be on improving/completing existing products from satellite ocean-color sensors, i.e., extending calculations of abovesurface fluxes to high Sun zenith angles and accounting, at least statistically, for diurnal variability of clouds. One should also work on how to compute similar products from different missions and their likely uncertainties, to obtain platformindependent global products. On the other hand, with existing scientific knowledge, the following new products (see sections Current Products and Users Needs) should be derived from past and current missions: planar and scalar PAR below the surface, average cosine for PAR below the surface, diffuse attenuation coefficient for PAR, and respective spectral quantities (visible to near infrared), spectral albedo, and APAR. These products should be in the form of daily averages, except for diffuse attenuation coefficient. New processing lines will require links with other ocean and atmosphere products and/or ancillary data (e.g., reanalysis products from observations and models, such as MERRA-2). A calibration/validation program should be planned to evaluate the new products and their uncertainty over a representative set of conditions.
In the medium term (longer-term effort, not lower priority), with specific efforts and new avenues in algorithm development, the aim should be to generate the following products, not only from past and current missions, but also from future missions: spectral fluxes below the surface and diffuse attenuation in the UV (or integrated over UV-A and UV-B ranges), in both photon and energy units, especially using TROPOMI on Sentinel 5P/5, average mixed-layer PAR (will require mixed-layer depth fields from Argo-assimilated circulation models), upperocean heating profile, and under-ice light fields (in conjunction with cryosphere missions and modeling). Diurnal-cycle resolving measurements, combining different satellites over-passing at different times with geostationary satellites, should be pursued  FIGURE 17 | Same as Figure 15 but for the ratio of spectral scalar PAR below the surface at 405 and 675 nm.
FIGURE 18 | Error budget of the daily "MERIS" PAR product above the surface (estimated PAR-actual PAR) from simulations for the period 2003-2012 using 1-hourly resolved MERRA-2 input data. For each day, a set of MERIS spectral reflectance data is simulated for a typical observation at 10:30 UT local time and several viewing geometries (nadir and 20 • view zenith angle (VZA) with relative azimuth of 0, 90, and 180 • , with sun glint avoidance). (Left) 2D plot of the error bias as a function of the clear sky daily PAR above surface (x axis) and cloud factor (y axis), along with 1D (marginal) distribution. (Right) Same as left but for the error standard deviation. For the 1D marginal distribution as a function of the clear-sky PAR is also reported the monthly PAR error bias and standard deviation in magenta.
to describe hourly changes in radiation fields and improve dailyintegrated values (e.g., due to clouds). The products should also be generated without gaps (applying gap-filling techniques) to provide boundary conditions for general circulation models. At this stage, evaluation of the products and their uncertainty should be ongoing on a continuous basis.
In the long term (future vision), for significant improvement of sub-surface light fields, space lidars could be used (e.g., the CNES MESCAL), as they can resolve the vertical distribution of material in the ocean (while all the products described above, derived from passive optical sensors, assume a homogeneous upper ocean). This is particularly critical in high-latitude regions (near the ice) and near land. Approaches using hyper-spectrally resolved sensors such as SCHIMACHY to retrieve the availability of light in the ocean (depth-integrated scalar irradiance) from the vibrational Raman scattering effect of water molecules (Dinter et al., 2015) should be explored. Finally, satellites missions with instruments in L1 orbit (such as the NASA EPIC onboard DISCOVR) would offer the opportunity to continuously observe the sun-lit part of the ocean, maximizing the temporal coverage. Space agencies should consider exploiting using such orbit for an ocean-color satellite mission.

AUTHOR CONTRIBUTIONS
All the authors contributed to the conception and organization of the study, to the analysis of the various products and results, and to the writing of the findings. RF, DR, and EB arranged the contents of the various sections, and wrote the first version of the manuscript. HB provided insight into light and biogeochemical processes. JT and BF carried out the work related to current MODIS and OMI products. DR, DJ, MC, and RF developed algorithms for new products and devised ways to associate uncertainties. DR, DJ, and MC generated examples of new products. TJ, TP, and SS organized the survey about satellite PAR observations and compiled and synthesized user responses. All the authors critically revised the manuscript and helped to bring about its final version.

FUNDING
The National Aeronautics and Space Administration (NASA) and the European Space Agency (ESA) provided funding for the study under Grant NNX14AL91G and Contract 40001 1 3690/1 4/I-LG, respectively.
-Spectral downward planar irradiance cos (θ ) L 0± (λ, t, θ , φ) sin (θ ) dθ dφ -Spectral upward planar irradiance where λ i is the central wavelength of any spectral interval λ in the PAR spectral range, h is Plank's constant, c is the speed of light, and N A is the Avogadro number.

Average Cosines (Dimensionless)
-Spectral average cosine 24h-Averaged Quantities where N sec = 86400 s.day −1 and t is expressed in day.