Evaluation of the NASA OBPG MERIS ocean surface PAR product in clear sky conditions

: The operational MEdium Resolution Imaging Spectrometer (MERIS) daily mean photosynthetically available radiation (PAR) product generated by the NASA Ocean Biology Processing Group (OBPG) was evaluated in clear sky conditions against in-situ measurements at various sites in the northwestern Mediterranean Sea (BOUSSOLE buoy), the northwestern Paciﬁc (CCE-1 and -2 moorings), and the northeastern Atlantic (COVE platform). The measurements were ﬁrst checked and corrected for calibration errors and uncertainties in data processing by comparing daily means for clear days (i.e., no clouds from sunrise to sunset and low aerosol abundance) with theoretical values from an accurate Monte Carlo radiative transfer code. The OBPG algorithm performed well when sky was completely cloudless during daytime, with a bias of 0.26 E/m 2 /d (0.6%) and a RMS diﬀerence of 1.7 E/m 2 /d (4.0%). Using satellite-derived aerosol optical thickness (AOT) and Angström coeﬃcient instead of climatology slightly degraded the results, which was likely due to uncertainties in the aerosol retrievals. A sensitivity study to aerosol properties indicated that climatology may not work in some situations (e.g., episodic dust, pollution, or biomass burning events), suggesting that it is best to use actual aerosol estimates in clear sky conditions. The analysis also revealed that specifying aerosol properties, therefore atmospheric transmittance, from AOT and Angström coeﬃcient, even retrieved from the satellite imagery, may not be suﬃcient in the


Introduction
Photosynthetically available radiation (PAR) at the ocean surface, i.e., the solar energy flux reaching the surface in the spectral range 400-700 nm, controls the growth of phytoplankton The paper is organized as follows. Section 2 gives a brief overview of the OBPG PAR algorithm. Section 3 describes the in-situ and satellite evaluation datasets used in the study. Section 4 details the procedures to select clear (i.e., non-cloudy) days based on in-situ measurements and presents the comparison against in-situ measurements of OBPG MERIS-derived clear-sky daily PAR estimates obtained using either aerosol climatology or satellite observations of aerosol properties. Section 5 discusses the impact of aerosol properties and satellite overpass time on the daily PAR estimates and proposes a methodology to minimize errors associated with clouds occurring before or after the satellite observation. Finally, Section 6 summarizes the collection of evaluation studies and the methodology for improving the current OBPG PAR product with respect to diurnal variability of clouds, not only for clear sky situations at time of satellite overpass, but for all cloudiness conditions.

Overview of OBPG PAR algorithm
The OBPG PAR algorithm [10] estimates daily mean (i.e., 24-hour averaged) PAR reaching the ocean surface in Einstein per meter squared per day, i.e., E/m 2 /d. It is based on plane-parallel theory and assumes that the effects of clouds and other atmospheric constituents can be decoupled. Therefore, the planetary atmosphere is modeled as a clear-sky atmosphere positioned above a cloud/surface layer and the surface PAR is expressed as the product of a clear-sky component and cloud transmittance. A budget approach is used, in which cloud transmittance is essentially obtained by subtracting the flux reflected to space (estimated from the satellite measurements) from the incoming extraterrestrial flux. One advantage of this approach, compared with approaches that require a knowledge of cloud properties (i.e., fraction, optical thickness), is that clear and cloudy regions within a pixel do not need to be distinguished. The scheme proceeds as follows.
Under solar incidence θ s , the top-of-atmosphere (TOA) solar flux, E 0 cos(θ s ), is reduced by a factor T g T a (1 -S a A) −1 by the time it reaches the cloud/surface system. In this expression, E 0 is the solar constant corrected for Earth-Sun distance, T g is the transmittance due to absorbing gases (essentially ozone), T a is the clear sky total (direct + diffuse) transmittance due to aerosols and molecules, S a is the spherical albedo of molecules and aerosols, and A is the cloud/surface system albedo. As the flux passes through the cloud/surface system, it is further reduced by a factor A. Therefore, the instantaneous solar flux reaching the ocean surface at any wavelength is given by where A s is the albedo of the ocean. In this expression, all variables depend on time t and wavelength λ, which is not indicated for clarity. The term (1 -A s ) in Eq. (1) is introduced to correct the flux penetrating the surface, the quantity naturally obtained, for reflection by the surface and water body. Equation (1) would be exact under the de-coupling assumption if the cloud/surface layer were reflecting isotropically, but this is generally not the case. In order to compute E d , A is expressed as a function of the radiance measured by the satellite sensor in the PAR spectral range. i.e., bands centered on 412.5, 442.5, 490, 510, 560, 620, and 665 nm for MERIS. These spectral bands do not saturate over clouds. First a land mask and a sea ice mask (generated from daily NSIDC data) is applied. Then for each pixel not contaminated by Sun glint, the measured radiance, after transformation into reflectance and correction for gaseous absorption, is corrected for molecular and aerosol effects to yield the reflectance of the cloud/surface layer. For this, quasi single scattering approximation is used, and aerosol optical thickness and Angström coefficient are specified from a monthly satellite climatology based on three years (1998)(1999)(2000) of Level 3 SeaWiFS data. Aerosol single scattering albedo and phase function are determined from the aerosol model selected in a base of standard aerosol models [24][25] via the Angström coefficient. This procedure is justified because aerosol properties cannot be determined when the pixel is cloudy and, in general, aerosol effects on E d are secondary compared with cloud and θ s effects.
Once the reflectance of the cloud/surface layer, ρ, is determined, it is converted into albedo A. This is accomplished for all the spectral bands by applying a cloud bidirectional correction factor F (independent of wavelength) to ρ -A s (since A ≈ NA c + A s , where N is fractional cloud coverage and A c cloud albedo; see [26] and Frouin et al. [10]), i.e., A = F(ρ -A s ) + A s . Analytical formulas developed by Zege et al. [27] for non-absorbing and optically thick scattering layers are used for F. The surface albedo, A s , is parameterized as a function of solar zenith angle and aerosol optical thickness at 500 nm, following Jin et al. [28]. In the look-up tables, however, the dependence of A s on wind speed and chlorophyll concentration is averaged. If ρ > A s , a formula for completely diffuse sky is used.
The daily mean PAR is finally obtained by integrating over the length of the day, i.e.
where t is expressed in hour, and E PAR denotes the instantaneous flux. The dependence of A, A s , T g , and T a on solar zenith angle is taken into account in the integration (S a does not vary with angular geometry), but the cloud properties are assumed to be constant during the day and to correspond to those of the satellite observation. This assumption is crude, and consequently the accuracy on daily-averaged quantities is degraded in regions where clouds exhibit marked diurnal variability.

In-situ measurements
Long-term measurements acquired at four mid-latitude oceanic sites were used in this study to evaluate the performance of the OBPG algorithm in estimating daily mean PAR. Geographic location of the sites, referred to as BOUSSOLE, CCE-1 and -2, and COVE, is displayed in Fig. 1, and main characteristics of the data sets are given in Table 1. Since these sites are located between 33 and 44°N, situations for which solar zenith angles stay large during the day, such as in polar regions, were not sampled, but atmospheric conditions exhibited large variability.  . The solar flux between 400 and 700 nm was directly measured every minute using two Li-Cor LI-190 sensors. When PAR data from both sensors were simultaneously available, they were averaged to yield the final measurement.

Satellite observations
Level-2 daily mean PAR imagery from MERIS data at about 1 km spatial resolution were acquired from the NASA OBPG website. For each in-situ site, a 3 × 3 box centered on the pixel containing the in-situ site was first extracted from the imagery, then the daily mean PAR values at the exact site location was obtained via nearest-neighbor interpolation. Since consecutive pixels in 3 × 3 boxes had approximately the same size the center pixel was always selected. Finally, all the pixels flagged cloudy for atmospheric correction were eliminated. (Pixels contaminated by Sun glint were already screened in the NASA OBPG PAR processing.) This match-up procedure is different from the one used by Harmel and Chami [15], who averaged all cloud-and glint-free pixels in a 5 × 5 box to obtain the satellite PAR estimate. It ensured that all pixels containing the in-situ site were clear.
In addition, information about chlorophyll-a concentration, aerosol optical thickness at 865 nm, AOT(865 nm), and Angström coefficient was obtained in a similar way using Level-2 products from Sea-viewing Wide Field-of-view Sensor (SeaWiFS), MODerate resolution Imaging Spectrometer (MODIS) onboard both Aqua and Terra, MERIS, and Visible Infrared Imaging Radiometer Suite (VIIRS) onboard S-NPP. At each site, the nearest-neighbor pixels from the various sensors (acquired between about 10:00 and 13:30 local time depending on the sensor) were selected and arithmetically averaged to represent average conditions during the day. This treatment was used in several studies: 1) to compare radiation transfer codes (Section 4.2), 2) to adjust biases (Section 4.3), and 3) to evaluate the impact of using actual aerosol information during the day of the satellite observation instead of climatology on the daily mean PAR estimates (Section 4.3).

Selection of clear days
To evaluate the performance of OBPG PAR model in clear sky conditions, i.e., without clouds during the entire day, the first task was to select such clear days. This was accomplished by fitting the in-situ data with the theoretical E PAR curve, expressed with a good approximation, except at very slanted solar angles, as E PAR = a cos(θ s ) exp[-b/cos(θ s )] [29]. If the sky is clear during the entire day, and aerosol diurnal variability is not too large, the in-situ measured PAR should fit this equation closely. The following selection criteria, based on the parametric regression above and other considerations, were used. Using these criteria, a total of 229 clear days were finally selected for the four sites during the time periods indicated in Table 1, which include many days after the end of the MERIS mission in May 2012. Figure 2 shows examples of clear-day selections at individual sites. At the CCE-1 and -2 sites, in particular, observations were generally not available during the early and late day time hours, and it is possible that clouds might have been present during those hours, but the impact of such observations (characterized by large solar zenith angles) on daily means is expected to be secondary in most cases. Satellite-derived aerosol optical thickness, determined from multiple satellite observations (see Section 3.2), was generally low during those days, i.e., AOT(865 nm) < 0.2 in most cases, but Angström coefficient was in the range 0-2, with a substantial amount of values above 1 (Fig. 3). This suggested some influence from land (pollution-type, biomass burning, or continental aerosols), in agreement with Knobelspiesse et al. [30], who found similar statistics within a 300km distance from land, but from a large data set of sunphotometer measurements made onboard ships during research and opportunity cruises.

Bias adjustment
Instrument calibration errors (and other errors, e.g., due to the data processing) might exist, affecting the quality of the in-situ measurements. It is important, therefore, to check the calibrated  data sets and remove eventual biases before evaluating the MERIS-derived daily mean PAR estimates. This was accomplished by comparing measured values with simulations from a Monte Carlo radiation transfer (RT) code, i.e., SMART-G (Speed-Up Monte Carlo Advanced Radiative Transfer using GPU; Ramon et al. [31]). The comparisons were made for clear sky situations, selected as described above. Simulations are expected to be accurate in such conditions, especially when aerosol loading is small, since the flux arriving at the surface is then essentially controlled by the solar zenith angle and the molecular optical thickness, quantities that can be determined with little uncertainty.
SMART-G has been evaluated for various benchmarks [31], yet in the present study it was compared with two other codes, namely ARTDECO (Atmospheric Radiative Transfer Database of Earth and Climate Observations, Dubuisson et al. [32]), and 6S (Second Simulation of the Satellite Signal in the Solar Spectrum; Vermote et al. [33]). The comparison was made for the clear sky days at the COVE site (128 days). In the calculations, the instantaneous spectral surface irradiance in the PAR range (code output) was spectrally and temporally integrated from sunrise to sunset to yield daily mean values. Surface albedo was specified according to Jin et al. [28] as a function of solar zenith angle, satellite-derived chlorophyll concentration and aerosol optical thickness (see Section 3.2), and wind speed (ancillary data associated with the satellite ocean color products). Ancillary ozone and water vapor content was also used in the calculations. Aerosol optical thickness was assumed to decrease exponentially with altitude (scale height of 2 km), and the aerosol model was specified as a fraction of maritime and continental models based on Angström coefficient. As shown in Fig. 4, the SMART-G and ARTDECO values are in very good agreement, with a small bias of -0.45 E/m 2 /d (about -1%) and a root mean squared (RMS) difference (of only 0.81 E/m 2 /d (∼2%). The relative bias and RMS difference between the SMART-G and 6S simulated data are larger, i.e., about -2.5% and 5%, respectively. This may be due to approximations in the 6S calculation of the atmospheric transmission functions and to the decoupling of gaseous absorption and molecular scattering. In SMART-G and ARTDECO, on the other hand, absorption and scattering processes are fully coupled, and the scattering treatment is accurate at small to large zenith angles (assuming the atmosphere is plane-parallel, the case of the simulations). Using SMART-G and the above procedures, the daily mean PAR was calculated for the clear days at the other sites as well, i.e., 229 days total, and compared with the corresponding in situ data. Days with AOT(865 nm) < 0.08 were selected, since for those situations the aerosol influence is minimized and, therefore, the calculations are most accurate. Figure 5, left panel, displays theoretical versus measured values, showing differences increasing with daily mean PAR values, which suggests that changing the calibration coefficient would rectify most of the discrepancies. Table 2 gives the comparison statistics for individual sites, indicating relative bias and RMS difference of approximately 5 to 8% and 5 to 9%, respectively, depending on the site. The measured daily clear sky PAR values were then regressed against the calculated values from the Monte Carlo code, and corrections were applied to the measurements based on the best fits. An order 2 polynomial (y = ax 2 + bx) was used. In addition, any noticed temporal dependence in the corrected data was removed. After correction, the relative bias between measured and simulated values dropped to negligible, and the relative RMS difference also decreased, to a maximum of 3.6% and a minimum of 1.4% depending on the site (Table 2). Figure 5, right panel, shows the improvement after bias adjustment, made for each site separately. The overall bias decreased from 2.64 E/m 2 /d (6.1%) to 0.01 E/m 2 /d (0.01%) and the RMS difference decreased from 3.74 E/m 2 /d (8.6%) to 1.42 E/m 2 /d (3.1%). We performed the same analysis using all the clear sky cases, i.e., without selecting smallest AOT(865 nm) values and obtained similar results (Fig. 6). The bias and RMS difference after correction were slightly larger, i.e., -0.14 E/m 2 /d (-0.3%), and the RMS difference was 1.61 E/m 2 /d (3.3%). This indicates that relatively high AOT(865) values are acceptable for adjusting the data, but since the number of observations with AOT(865 nm) < 0.08 was sufficiently large, we used the results obtained with the reduced data set. The correction of the in-situ data was then also applied to all the cases considered in the study, including when clouds were present before and after overpass time.   6. Same as Fig. 5, but all AOT situations. The regression obtained for the reduced dataset, i.e., AOT(865 nm) < 0.08, see Table 2, are applied in the adjustment of the measured values at the various evaluation sites.

MERIS performance
For the 4 sites, match-ups of level-2 MERIS-derived daily mean PAR data and corrected in-situ measurements were determined as described in Section 3.2. A total of 43 match-ups were obtained in clear sky conditions (i.e., no clouds during the entire day). The standard OBPG MERIS PAR product was used for the satellite estimates. Note that this product uses aerosol climatology. Figure 7, Fig. 5, center panel). These authors, however, did not check the in-situ measurements for calibration and data processing errors and, perhaps more importantly, did not exclude situations with clouds occurring before and/or after the MERIS observation time. In such situations, one expects the OBPG PAR algorithm to yield higher values since the atmosphere, clear at the time of the MERIS measurements, is assumed cloudless during the entire day. The impact of cloud presence during part of the day, but not at the time of satellite overpass, will be addressed in Section 5.2. Harmel and Chami [15] further indicated that taking into account actual aerosol properties, i.e., retrieved from the MERIS data instead of using aerosol climatology, should improve the accuracy of the daily mean PAR estimates, and they provided evidence of such improvement (their Fig. 5, right panel). The NASA OBPG algorithm was therefore run with the AOT and Angström exponent determined from multiple satellites for each clear sky match-up (see processing details in Section 3.2). The results (Fig. 7, right panel) show, however, a slightly degraded performance, with bias of 0.56 E/m 2 /d (1.3%) and RMS difference of 1.91 E/m 2 /d (4.4%). This may be due to uncertainties in the aerosol retrievals, which are relatively large when aerosol optical thickness is small. The results again point out that for such aerosol situations (i.e., AOT(865) < 0.1), using aerosol climatology is sufficient to obtain accurate daily mean PAR estimates. However, long-term climatology cannot capture transient events such as dust storms, pollution outbreaks, biomass burning events, or volcanic emissions, and using estimates of actual aerosol properties at the time of satellite overpass (or during the day of the observations) is definitely the way to proceed, as noted by Harmel and Chami [15].

Sensitivity to aerosols
The influence of the type and loading of aerosols on the daily PAR was further examined theoretically using multiple 6S simulations. The simulations were performed for the BOUSSOLE site on two dates, i.e., January 1 and July 1 of 2000. Four aerosol models (maritime, continental, desert, and biomass burning) were used and aerosol optical thickness at 550 nm, AOT(550 nm), was varied from 0 to 1. Aerosol scale height was 2 km. Wind speed was set to 5m/s, chlorophyll concentration to 0.1 mg/m 3 , and ozone content to 0.34 atm.cm. Figure 8 shows how daily mean PAR and relative daily mean PAR difference change with AOT(550 nm). For situations generally encountered over the oceans (AOT(550 nm) < 0.2; see Knobelspiesse et al. [30]), the relative daily mean PAR difference varied within 10% on January 1 and within 5% on July 1 depending on the model, which indicates that the effect of aerosol type, although noticeable, is relatively small. The larger influence in winter is due to larger solar zenith angles during the day, therefore airmass, amplifying the sensitivity of atmospheric transmittance to aerosol type. As AOT(550 nm) increases, the impact of aerosol composition on daily mean PAR becomes increasingly large. When AOT(550 nm) is 1, the relative daily PAR difference is ∼20% on January 1 and ∼15% on July 1. A parameterization of atmospheric transmittance only based on aerosol optical thickness and Angström coefficient, as in the OBPG PAR algorithm, may not distinguish properly the effect of maritime and dust aerosols, for example, the latter absorbing yet exhibiting similar AOT spectral dependence. The good performance of the OBPG PAR algorithm displayed in Fig. 7 is likely due to the lack of (or small) occurrence of absorbing aerosols at the evaluation sites.
With increasing dust-like component and soot in the aerosols (e.g., continental aerosol model), the relative reduction in daily mean PAR reached ∼50% on January 1 and ∼30% on July 1 when AOT(550 nm) increased from 0 to 1. Therefore, while aerosols are a minor issue handled adequately by the OBPG PAR algorithm in most cases, substantial errors can be made on daily mean PAR estimates using climatology when aerosol abundance is large, all the more as the aerosols are absorbing (e.g., dust events). This indicates that it is best to use satellite-derived information about aerosols in clear sky situations, as done with the OLCIPAR algorithm proposed by Harmel and Chami [15], while keeping climatology for cloudy regions where the impact of AOT is relatively small compared with the effect of clouds. This procedure, however, might not be effective when clouds occur before or after satellite overpass, and such cases require special consideration and treatment (see below). Furthermore, specifying aerosol properties at the satellite observation time neglects diurnal variability (time scale of dust or pollution events may be as short as a few hours).

Impact of overpass time
As indicated before, the computation of clear-sky daily mean PAR using the OBPG PAR algorithm assumes that the sky, clear at the time of satellite overpass, is clear during the entire day. Clouds, however, may exist during part of the day, and their properties exhibit strong diurnal variability. Figure 9 shows examples of in-situ measured PAR at the COVE site when the MERIS observation occurs in clear sky conditions, but clouds exist at other times. The clear sky conditions are identified by the MERIS cloud masking used to generate chlorophyll-a concentration. Since the daily mean PAR is estimated from one MERIS observation per day, the cloud diurnal variability is not taken into account, resulting in estimates much higher than the measured values, by 27.7, 17.   The satellite estimates are biased high on average, by 2.8 E/m 2 /d (7.3%) due to the inability of the PAR algorithm to detect cloudiness at times other than the time of overpass, and the RMS difference is relatively large, i.e., 6.0 E/m 2 /d (15.8%). These results are similar to those reported by Harmel and Chami [15] for comparisons at the COVE site (349 match-ups), except that the bias in Fig. 10, left panel, is much smaller (2.8 instead of 4.4 E/m 2 /d). These authors suggested that the single-scattering approximation used to correct the effect of the clear sky atmosphere in determining A (see details in Section 2) might explain the bias, but the likely explanation, as mentioned above, is that clouds occurring before or after the time of satellite overpass reduced the daily mean PAR. Note that in situ measurements, biased low when compared with Monte Carlo calculations (Fig. 5 and Table 2), were not corrected in the Harmel and Chami [15] study, suggesting a much better performance by the OBPG algorithm, with an overestimation less than 1 E/m 2 /d.
To reduce biases due to diurnal variability of clouds, one needs to account for changes in the relevant cloud properties, i.e., fractional coverage and optical thickness. As suggested by Frouin et al. [8], this can be accomplished using reanalysis data such as the Modern Era Retrospective Reanalysis for Research and Applications, Version 2 (MERRA-2) data [34]. Cloud information is provided hourly at 0.5°x0.625°spatial resolution. The grid size of the MERRA-2 products is indeed coarse compared with the size of a MERIS pixel, therefore the information may not represent well actual conditions in some cases (sky may be clear locally, but cloudy away within a grid cell), but statistically one may expect improvements by taking this information into account. If τ c−MERRA and N MERRA denote the MERRA-2 cloud optical thickness and fractional coverage, A will be replaced by A' in Eq. (3) as follows: where t sat is the time of satellite observation and t the time of the MERRA-2 data. Basically, the quantity 1 -A + A s is adjusted by the ratio of 1 -NA c computed from MERRA-2 products at times t and t sat . The parameterization of Fitzpatrick et al. [35] is used to obtain A c from τ c−MERRA and solar zenith angle. In Eq. (4) A, determined at t sat , is corrected for variations with solar zenith angle, but assuming unchanged cloud properties. This assumes (as done to transform reflectance of the cloud/surface system to albedo, see Section 2) that for a pixel characterized by N, A s , and A c , A is approximated by NA c + A s (second and higher orders terms in A c and A s are neglected). The right-hand side of Eq. (4) reduces to A(t sat ) when t becomes t sat , i.e., the treatment keeps into account the estimated A at the time of satellite observation. Figure 10, right panel displays the estimated versus measured daily mean PAR at the four sites after the treatment for cloud variability, i.e., after Eq. (4) is applied. In the calculations, N was estimated from the MERRA-2 fractional coverage of low, middle, and high clouds, N i, using the maximum overlap assumption (e.g., [36]) as N = i (1 -N i ). The cloud variables of the grid cells that included the geographic location of the sites were used without interpolation. The bias and RMS difference are reduced to 1.6 E/m 2 /d (4.2%) and 4.8 E/m 2 /d (12.4%), respectively. Very little change was obtained for situations of clear sky during the entire day when using MERRA-2 data (points depicted in red in Fig. 10), indicating that the relatively coarse grid size was not introducing significant degradation in performance for those situations. Despite the general improvement, some estimates were still too high after treatment, which is plausibly explained by the MERRA-2 data not representing properly the cloud properties at the sites. Note that the Fig. 11. Cloud properties, i.e., N (red lines), A c (blue lines), and 1 -NA c (i.e., basically cloud transmittance, green lines), as a function of time during the day (sunrise to sunset at the COVE site for the selected days of Fig. 9, as determined from MERRA-2 data. Time of MERIS overpass is depicted by a vertical line. The MERIS-derived daily mean PAR obtained by taking into account diurnal cloudiness is indicated with the corresponding in-situ value. Agreement with in-situ measurements is closer using the MERRA-2 data. results were practically the same (i.e., bias and RMS difference smaller by only 0.1%) when using, for the cloud variables, a linear interpolation of data in a 3 × 3 MERRA-2 cell box. Figure 11 displays, for the Fig. 9 examples, the time evolution of N, A c , and 1 -NA c during the day, from sunrise to sunset. The variable 1 -NA c represents the transmittance (due to clouds) of a partly cloudy sky. Except for March 26, 2005 (Fig. 11(b)), the MERIS-derived daily mean PAR values are in close agreement with the in-situ measurements, i.e., within a few E/m 2 /d. For October 3, 2004, in particular, the dramatic overestimation of 27.7 E/m 2 /d (198%) is reduced to 1.9 E/m 2 /d (9%). The MERRA-2 data captured well the increase in N and A c and, therefore, the reduction of 1 -NA c , after the MERIS overpass. A similar analysis explains the improvement obtained for December 4, 2009 andSeptember 7, 2011 (Figs. 11(c)and (d)). For March 26, 2005, on the other hand, the improvement is relatively small, i.e., only 2 E/m 2 /d. During that day, the MERRA-2 data indicated reduced cloudiness after the MERIS overpass until mid-afternoon, which is unrealistic since sky was clear at the time of overpass, therefore compensating a decrease in PAR in the afternoon due to increased N and especially A c . This illustrates the limitation of using reanalysis data at coarse resolution, which may not capture properly the cloud properties at the scale of a MERIS pixel. Nevertheless, the significant improvement (in the statistical sense) revealed in Fig. 10, right panel, justifies the approach. It should be further noted that the treatment using MERRA-2 data for diurnal cloud variability should work not only for clear satellite pixels, but also cloudy pixels, i.e., is applicable generally, and should be considered for operational processing. The drawback, however, is that MERRA-2 products are generated with some delay, even though relatively short (i.e., one month), but it is conceivable to refresh regularly the PAR products generated by the standard NASA OBPG algorithm as the MERRA-2 data become available.

Summary and conclusions
The operational MERIS PAR product generated by the NASA OBPG was evaluated in clear sky conditions, i.e., for situations without clouds during the entire day or at the time of satellite overpass. Long-term in-situ measurements obtained at four sites (BOUSSOLE, CCE-1, CCE-2, and COVE) were used in the comparisons with MERIS estimates. Selection of completely clear sky situations from sunrise to sunset was accomplished by fitting the measurements by a theoretical curve expressing PAR as a function of solar zenith angle. Possible errors in the field measurements, due to radiometric calibration and data processing, were checked by comparing the daily mean PAR values of the completely clear days with those computed by a Monte Carlo radiative transfer code. The impact of aerosols was minimized by selecting cases with low AOT(865 nm). For clear days with AOT(865) < 0.08, the bias between Monte Carlo modeled and in-situ measured daily mean PAR was 2.6 E/m 2 /d (6.0%), i.e., the measured values were significantly lower, and the RMS difference was 3.7 E/m 2 /d (8.6%). After correction, the bias decreased to 0.01% and the RMS difference to 3.3%.
Match-ups of MERIS and in-situ observations for completely clear sky situations during daytime showed that the NASA OBPG algorithm performed very well, with a bias of 0.26 E/m 2 /d (0.6%) and RMS difference of 1.7 E/m 2 /d (4.0%). This performance is much better than the performance reported by Harmel and Chami [15], but these authors did not exclude situations of cloudy sky before and/or after satellite overpass and did not account for potential calibration errors and data processing uncertainties in the in-situ measurements. Since the default NASA OBPG algorithm is based on aerosol climatology, it is expected that using actual aerosol properties should yield more accurate daily mean PAR estimates, as suggested by Harmel and Chami [15]. However, when using satellite-derived AOT and Angström coefficients as input to the algorithm, the PAR estimation was slightly degraded, with the bias of 0.6 E/m 2 /d (1.3%) and RMS difference of 1.9 E/m 2 /d (4.4%). This is most likely due to uncertainties in aerosol retrievals and the generally low AOT values during the selected clear days at the evaluation sites. Nevertheless, the use of satellite aerosol properties is recommended, especially to account for episodic aerosol events (e.g., dust plumes, biomass burning emissions, pollution outbreaks) not captured in aerosol climatology, and could be easily implemented in the NASA OBPG PAR algorithm since ocean-color and aerosol variables are retrieved concomitantly.
The study also examined the impact of aerosol abundance and composition on the estimation of daily mean PAR. Theoretical simulations showed that for most cases encountered over the oceans (i.e., AOT(550 nm) < 0.2) the relative daily mean PAR difference varied by less than 10% for different types of aerosols at the BOUSSOLE site. However, as AOT(550 nm) increased, aerosol type impacted significantly the daily mean PAR, with relative difference reaching about 20% on January 1 and 15% on July 1. It was also noticed that the daily mean PAR was reduced by 30 to 50% when AOT(550 nm) increased from 0 to 1, the larger reduction being obtained for absorbing aerosols. This indicates that specifying aerosol properties, therefore atmospheric transmittance, from AOT and Angström coefficient, as done in the NASA OBPG PAR algorithm, may not be sufficient in the presence of such aerosols, especially when loadings are important, even if these variables are retrieved from the satellite imagery. Since aerosol properties cannot be retrieved in cloudy conditions, using aerosol climatology is appropriate in those conditions, all the more as cloud effects generally dominate.
In addition to aerosol properties, sensor observation time has an impact on the daily mean PAR estimation when clouds are present before and/or after satellite overpass, a consequence of using only one MERIS observation per day and assuming that atmospheric properties do not change during the day. For all cases that exhibited clear sky at the time of satellite overpass (436 cases at the four sites), the statistical performance against in-situ measurements was degraded, with a bias (overestimation) of 7.3%, in better agreement with the results obtained by Harmel and Chami [15] at the BOUSSOLE site, although their reported bias was more than twice as large. Adjusting the in-situ measurements, biased low against "exact" Monte Carlo calculations, would have reduced the differences substantially. Our results suggest that the key to reducing biases when sky is clear at the time of satellite overpass but contaminated by clouds at other times during the day resides not so much in better specifying aerosol properties or treating atmospheric radiation transfer more accurately, but rather in taking into account properly the diurnal variability of clouds. A methodology was proposed to account for such variability, that utilizes hourly MERRA-2 data (cloud fractional coverage and optical thickness). Despite the difficulty of accurately depicting cloud properties at the size of a MERIS pixel from the coarse resolution reanalysis data, the bias was reduced to 4.2%, with negligible change for situations of clear sky during the entire day. The treatment is applicable operationally, not only to clear pixels but also cloudy pixels, and should be considered in a future version of the NASA OBPG PAR algorithm, although delay in the availability of the MERRA-2 products, even relatively short, would require refreshing the standard PAR products accordingly.