Assessment of the spectral downward irradiance at the surface of the Mediterranean Sea using the radiative Ocean-Atmosphere Spectral Irradiance Model (OASIM)

. A multiplatform assessment of the Ocean– Atmosphere Spectral Irradiance Model (OASIM) radiative model focussed on the Mediterranean Sea for the period 2004–2017 is presented. The BOUée pour l’acquiSition d’une Série Optique à Long termE (BOUSSOLE) mooring and biogeochemical Argo (BGC-Argo) ﬂoat optical sensor observations are combined with model outputs to analyse the spatial and temporal variabilities in the downward planar irradiance


Introduction
The availability of in situ oceanic radiometric data has recently increased owing to the deployment of autonomous robotic profiling platforms called biogeochemical Argo floats (hereafter referred to as BGC-Argo floats; Johnson and Claustre, 2016; please also see Appendix A for a list of abbreviations) equipped with radiometric sensors (Organelli et al., 2016). Such data may be exploited to improve the calibration and tuning of the bio-optical models embedded in threedimensional global and regional physical-biogeochemical coupled models. Furthermore, radiometric data availability will further increase with the development of new autonomous profiling floats dedicated to ocean colour measurements (Leymarie et al., 2018) and with the expanding data streams provided by the Ocean Colour Radiometry Virtual Constellation (OCR-VC) satellite. In turn, these new observations will pave the way towards the next generation of ocean biogeochemical models. State-of-the-art bio-optical algorithms, with the atmospheric component providing multispectral light boundary conditions at the sea-water interface with different levels of complexity, are already integrated into the biogeochemical models developed at the Massachusetts Institute of Technology (MIT, Dutkiewicz et al., 2015), Commonwealth Scientific and Industrial Research Organisation (CSIRO, Baird et al., 2016), and National Aeronautics and Space Administration (NASA, Gregg and Rousseaux, 2017). Additionally, the direct assimilation of optical/radiometric data  yields a higher robustness than traditional assimilation of modelled optical/radiometric data based on the chlorophyll a concentration due to the greater in-depth knowledge of the uncertainties in optical measurements (Dowd et al., 2014). In perspective, the evaluation of the uncertainty of the multispectral light at the ocean-atmosphere interface is important both for the solution of the radiative transfer model within the water column and for the development of assimilation schemes of radiometric parameters.
In this framework, the Mediterranean Sea appears to be a key area for the development of a new multispectral biooptical model to be potentially integrated with the European Copernicus Marine Environment Monitoring Service (CMEMS). A dense network of BGC-Argo floats providing quality-controlled radiometric data has indeed been deployed in the Mediterranean Sea (Organelli et al., 2016(Organelli et al., , 2017, which can be combined with the high-frequency radiometric and bio-optical observations acquired by a dedicated fixed platform (BOUée pour l'acquiSition d'une Série Optique à Long termE -BOUSSOLE; Antoine et al., , 2008. This meets the requirements and high data quality standard expected both for remote system calibration of ocean colour spaceborne sensors  and for the CMEMS biogeochemical operational model system for the Mediterranean Sea (MedBFM; Lazzari et al., 2010Lazzari et al., , 2012Lazzari et al., , 2016Cossarini et al., 2015;Teruzzi et al., 2014Teruzzi et al., , 2018Teruzzi et al., , 2019Salon et al., 2019). This system, recently upgraded to assimilate BGC-Argo float data , is used to produce analysis, forecasts and reanalysis of the biogeochemical state.
Assessment of the quality of the direct and diffuse downward irradiance produced by the multispectral Ocean-Atmosphere Spectral Irradiance Model (OASIM; Gregg and Casey, 2009) at the sea surface is a prerequisite to constrain bio-optical in-water light propagation modelling. This activity has been carried out within the framework of the CMEMS Service Evolution project BIOPTIMOD, which is aimed at the development of a new multispectral bio-optical model that will include the integration of MedBFM with data provided by both BGC-Argo floats and multispectral satellite sensors (e.g. the Ocean and Land Colour Instrument, OLCI, on board Sentinel3-A and Sentinel3-B; Donlon et al., 2012).
In MedBFM, the computation of the downward irradiance at the sea surface will be updated with OASIM, which has been pre-operationally interfaced with European Cen-  (left) and cloudy skies (right): E dir is the direct downwelling irradiance, E dif is the diffuse downwelling irradiance, and ρ dir and ρ dif are the direct and diffuse surface reflectances, respectively (the size of the arrows approximates the relative contributions; the figure is modified from Fig. 1 of Gregg and Rousseaux, 2017).
tre for Medium-Range Weather Forecast (ECMWF) atmospheric products and validated against reference data in the Mediterranean Sea, provided by the BOUSSOLE buoy and BGC-Argo floats.
This work also contributes to the extension of the validation of OASIM in the Mediterranean Sea, an area not covered by the previous skill assessment of Gregg and Casey (2009).
Section 2 presents OASIM, the required input data and the datasets adopted for validation purposes. The results are provided in Sect. 3 and examined in Sect. 4. Conclusions are summarized in Sect. 5. The full name of the abbreviations used in the present paper is provided in Appendix A.
2 Materials and methods

OASIM
OASIM (Gregg and Casey, 2009; hereafter referred to as GC2009) simulates the propagation of the downward spectral radiance in the atmosphere and provides the direct and diffuse irradiance over the ocean surface as the output (Fig. 1) at 33 wavelengths ranging from 200 nm to 4 µm (15 wavelengths at a 25 nm spectral resolution in the near-ultraviolet (UV) and visible regions of the light spectrum, 350-700 nm).
OASIM is based on the RADTRAN spectral model developed by Gregg and Carder (1990) for clear-sky conditions, with an upgraded aerosol parameterization scheme, and on the Slingo (1989) parameterization scheme for the spectral cloud transmittance and considers the spectral absorption and scattering of atmospheric gases (ozone, water vapour, oxygen and carbon dioxide). All model parameters (extraterrestrial solar irradiance, Rayleigh optical thickness and absorption coefficients for the various atmospheric gases) are detailed for each of the 33 bands in Table 2 of GC2009.
In OASIM, gaseous absorption by ozone, oxygen, carbon dioxide and water vapour is resolved before cloud transmittance determination, and aerosol effects are ignored in the presence of clouds. In the clear-sky parameterization scheme, the role of aerosols is described by three parameters: aerosol optical thickness (AOT), single-scattering albedo and asymmetry. These quantities have been applied in GC2009 exploiting Moderate Resolution Imaging Spectroradiometer (MODIS; Remer et al., 2005) data at seven wavelengths ranging from 470 nm to 2.13 µm (for the out-of-range wavelengths, linear extrapolation is performed) from February 2000 to July 2007, extended to 2017 in the present work.
The Slingo model requires four cloud properties as input: cloud cover, cov [%]; cloud liquid water path, LWP [g m −2 ] (ice clouds are not considered in OASIM); spectral cloud optical thickness, τ c (λ) [dimensionless]; and cloud droplet effective radius, r e [µm]. The last three properties are linked by the following expression (Slingo, 1989): where a(λ) and b(λ) are spectral cloud coefficients.
In the original formulation of OASIM presented in GC2009, cov and LWP data are retrieved from the International Satellite Cloud Climatology Project (ISCCP, from June 1986 to July 2002, and the prior climatology), while r e is parameterized using MODIS data normalized with mean values obtained from the literature (please refer to GC2009 for details). Specific modelling of the spectral reflectance of sea foam, affecting the transmittance across the air-sea interface, is also included in OASIM (please refer to the Appendix of GC2009).
In addition to the cloud properties necessary for the Slingo model (cloud cover and cloud liquid water path), OASIM requires the following atmospheric input data: surface pressure, sp [mb]; wind speed, ws [m s −1 ]; relative humidity, rh [%]; precipitable water (absorption by water vapour), wv [cm]; and ozone, oz [DU]. Except for ozone, which was obtained from the multiyear dataset of Total Ozone Mapping Spectrometer (TOMS) sensors (from 1979 to May 1993, hereafter referred to as the TOMS climatology), in GC2009, the other four parameters (sp, ws, rh and wv) were acquired from the National Centers for Environmental Prediction (NCEP 1979-July 2002Kalnay et al., 1996) reanalysis dataset. In the present implementation, most of these input data are extracted from the ECMWF ERA-Interim dataset (please refer to Sect. 2.4).
The spatial and temporal resolutions of OASIM are determined by the forcing datasets (i.e. aerosol optical data, cloud property data and atmospheric surface level data). The standard configuration presented in GC2009 is hereby maintained, with a 1 • horizontal resolution, thereby increasing the temporal frequency to 15 min to resolve the diurnal variability and compare the model output to the temporal resolution of the in situ data considered in the present study. High spatial resolutions and operationally oriented setups are of course possible, provided that forcing data are available: an upgrade in this direction is under investigation based on the ERA5 dataset recently updated and released (C3S, 2017;Hersbach et al., 2020).

The BGC-Argo float network in the Mediterranean Sea
The BGC-Argo float network represents the first ever nearreal-time (NRT) biogeochemical in situ large-scale ocean observing system (Johnson and Claustre, 2016). The technology is relatively recent, and floats equipped with sensors measuring the main biogeochemical and optical parameters have now become operational (Bittig et al., 2019). A BGC-Argo float operates similarly to an Argo float, collecting vertical profiles from 0-1000 m every 1 to 10 d and transmitting NRT data. The profiles are processed with a variable-specific quality control (QC) approach and are available within 24 h after data transmission. The BGC-Argo floats deployed in the Mediterranean Sea since 2012 are equipped with sensors, which, in addition to the temperature (T ), salinity (S) and depth, measure the chlorophyll a (Chl) fluorescence, downward planar irradiance (E d ) at three wavelengths (380, 412 and 490 nm) and downward photosynthetically active radiation (DPAR), which represents the integrated amount of the downward planar irradiance in the visible range, i.e. from 400 to 700 nm. Quality-controlled Chl, T and S NRT data are currently available from the Coriolis data centre (Argo Data Management Team, ADMT).
The photosynthetically active radiation (PAR) must be generally derived by spectrally integrating the scalar irradiance (i.e. the radiance integrated across the complete solid angle) because phytoplankton utilizes light from all directions. In the present case, BGC-Argo sensors measure the downward irradiance integrated in the 400 to 700 nm range: to avoid confusion, the measured quantity is referred to as DPAR.
QC of radiometric data (DPAR and E d ), specifically designed for in situ and remote-sensing ocean colour applications, has been implemented by Organelli et al. (2016). The delayed-mode (DM) QC approach to identify data corruption by biofouling and instrument drift based on tests and proce- dures has been developed by Organelli et al. (2017). In the present work, a BGC-Argo dataset was adopted covering the period between 2012 and 2017 with 3800 profiles (Fig. 2).
Before comparing model values to observations, the irradiance profiles obtained from floats were extrapolated to the surface with an exponential fitting procedure based on the curve_fit tool of the python package scipy. Further, we required profiles to have at least one measurement in the first 1.5 m depth from sea surface and to have at least four measurements in the first 10 m. In addition, any sub-basin (as defined in Fig. 2) and month containing fewer than five profiles was discarded.

The BOUSSOLE mooring buoy
The BOUée pour l'acquiSition d'une Série Optique à Long termE (BOUSSOLE) is a long-term mooring station collecting radiometric and bio-optical properties every 15 min from the topmost 10 m of the ocean layer (plus a reference, namely, the above-water measurement of the spectral downward planar irradiance) since 2003 . It is located in the Ligurian Sea at 43 • 22 N and 7 • 54 E, approximately 32 nautical miles off the French Riviera coast where the water depth is approximately 2440 m (Fig. 2). The measured quantities include the quality-controlled multispectral (nine bands) downward planar irradiance above the sea surface and DPAR, covering the period from 2003-2012.

OASIM input data for Mediterranean Sea applications
The input data required by OASIM in the present Mediterranean Sea application are shown in Fig. 3, together with the datasets used for the validation (BOUSSOLE and BGC-Argo). In the setup specific to the Mediterranean Sea, developed within the framework of the BIOPTIMOD project, the variables of the cloud and atmospheric properties are extracted from the European Centre for Medium-Range Weather Forecast (ECMWF) ERA-Interim reanalysis dataset (Dee et al., 2011); further details on their pre-processing are provided in Appendix B. The investigated period ranges from 2004 to 2017, with the input data provided as daily averages and with a spatial resolution of 1 • . The use of ERA-Interim to force OASIM is motivated by increasing the coherence between the input cloud properties and surface level atmospheric properties, which, in the GC2009 configuration, were provided by two different datasets (ISCCP and NCEP, respectively). Furthermore, considering that the optical module of the MedBFM model will be upgraded as part of the CMEMS operations, ERA-Interim was chosen since ECMWF products are already operational upstream data for the Mediterranean Sea regional production centre of the CMEMS (please refer to Salon et al., 2019). The MODIS data for the period from 2004-2017 are derived in the same way as reported in GC2009 but are extended to 2017.
The notation 0 − indicates quantities just below the sea surface with any reflections at the air-sea interface removed, and 0 + is applied for the radiances just above the air-sea interface (the light measured before air-sea transmission occurs).
The OASIM outputs for the irradiance are expressed in watts per quare metres per waveband, while the in situ data for the same parameter are expressed in watts per square centimetre per nanometre. Therefore, before comparison, the units were standardized to watts per square metre per nanometre. Because OASIM simulates data in the range of 350-700 nm centred in 25 nm bins, a linear interpolation was applied to match with measured wavelengths.
Furthermore, the direct (E dir (0 − )) and diffuse (E dif (0 − )) downward irradiance components simulated by OASIM were summed to compare them to the measurements of the downward planar irradiance sensors installed on the BOUSSOLE and BGC-Argo floats. Figure 3. Availability of the model input (cloud, aerosol and atmospheric data) and radiometric output data (E dir and E dif ) and corresponding observational datasets for the validation in the present work. The acronym definitions are reported in Sect. 2.1, 2.2 and 2.3. DPAR (µmol quanta m −2 s −1 ) was computed from the OASIM output (in standardized units) by integrating the downward planar irradiance as follows: where N A is the Avogadro number, h is the Planck constant and c is the speed of light. The definition of DPAR in GC2009 relied on integration with a lower bound of 350 nm, while the BOUSSOLE and BGC-Argo float sensors integrate from 400 nm. The model outputs were standardized according to the observations. To compute DPAR(0 − ), which is required for a correct comparison to BGC-Argo float sensors, Eq.
Hereafter, E d (0 + ) and E d (0 − ) indicate the sum of the direct and diffuse downward components.

ERA-Interim and MODIS atmospheric data analysis
A comparison of the ERA-Interim monthly averaged surface level atmospheric variables at grid points corresponding to the BOUSSOLE mooring buoy and the available observations (sea level pressure, wind speed and relative humidity) by the nearby Côte d'Azur buoy (Météo-France) for the period from 2004-2012 is shown in Fig. 4. An overall good agreement is observed for the surface pressure between the ECMWF data and in situ observations, while ERA-Interim underestimates the observed relative humidity (D'Ortenzio et al., 2008). The wind magnitude provided by the Côte d'Azur dataset is much larger than that of the ERA-Interim coarse data at 1 • , as previously reported by Stopa and Cheung (2014) using data from the US National Data Buoy Center. Large differences occur more frequently during the cold season (from November to April). Sensitivity tests to meteorological inputs were performed by Gregg and Carder (1990), showing that pressure and mean wind speed produced differences in surface spectral irradiance of less than 1 % in terms of model error over the 350-700 nm range, much less than air-mass type, visibility and total ozone.
In terms of the spatial distribution, the cloud cover follows a clear seasonal cycle (Fig. 5a), with low values during summer in the eastern sub-basins (ION1, ION2, LEV1, LEV2, LEV3 and LEV3) and a high cover during winter in the northern sub-basins (NWM, TYR1, TYR2, ADR1, ADR2 and AEG). The maximum monthly cloud cover during winter reaches approximately 50 %.
However, the aerosol optical thickness exhibits a different spatiotemporal pattern. The highest values are localized in the southwestern Mediterranean (SWM1) between July and August, possibly related to Saharan dust events in the area (Varga et al., 2014). High aerosol thickness values are also found in spring (April and May) in the eastern sub-basins (ION1, ION2, LEV3 and LEV4), which also most likely occurs due to aeolian dust transport (Antoine and Nobileau, 2006).

Validation of OASIM at the BOUSSOLE site
The present OASIM configuration produced data with a 15 min temporal resolution at the global scale on a 1 • mesh. As an example, the simulated daily cycle averaged over March is shown in Fig. 6. As expected, the measured data present a higher variability than that presented by the model, especially when the intra-monthly frequency is considered. The root mean square difference (RMSD) substantially decreases when considering the monthly average.
The year-by-year RMSD of the OASIM vs. BOUS-SOLE relationship remains steady at approximately 0.2 W m −2 nm −1 (Fig. 7, the red lines). The statistics (Fig. 7, the red dots) indicate a seasonal oscillation: in winter, the high RMSDs and low slopes imply a low prediction skill of the model, whereas in summer, low RMSDs and a slope near 1 indicate a high model skill. Furthermore, an averaged day in each month of the time series (average day per month) was computed, both for the data and model, discretized in 15 min temporal intervals. The same statistics were applied to the averaged data, and the results reveal a reduction in RMSD (Fig. 7, the cyan dots). The average-day-per-month filtering result indicates that the day-by-day variability is responsible for a substantial and consistent part of the  uncertainty (please refer to the averaged values shown in Fig. 7).
As an additional step, the averaged diurnal cycle grouping all the data in the same month based on the full time series was considered. For example, all the January data from 2004 to 2012 (9 years) were grouped, and the representative climatological day was derived based on a 15 min temporal resolution. In the case of January, at each wavelength and each 15 min interval, there exists a distribution consisting of 31 × 9 data points (including the data reduction due to sensor failure episodes). In principle, these distributions should be constrained by quite homogeneous conditions, with the same daily zenith angle component, and similar seasonal conditions, i.e. with all the data grouped by the month. Therefore, assuming that a large number of small cumulative perturbations affects the variability, the data should be log-normally distributed.
However, the Kolmogorov-Smirnov test revealed that neither the BOUSSOLE nor the model data distributions are log-normal. In fact, the distances between the accumulated empirical distributions and the reference lognormal distribution were almost always larger than the critical distance (D crit = 1.36/ √ N samples ; Bronshtein and Semendyayev, 2013), corresponding to the 5 % probability threshold to reject the null hypothesis (H 0 = the two distributions are the same). Moreover, analysis of the skewness and excess kurtosis confirmed that the data and model qualitatively exhibit similar distributions: both are negatively skewed, and the tails decay slower than does a Gaussian distribution (images not shown).

Validation of OASIM with the BGC-Argo float network
The quality-controlled BGC-Argo float dataset adopted in this work contains radiometric measurements acquired from 10:00 to 14:00 local time. To compare OASIM with these BGC-Argo measurements, a point-by-point match-up analysis was performed, where the closest model output to the float measurement at the surface was selected. The individual match-ups were then spatiotemporally aggregated based on the climatological months and the defined 16 sub-basins, as shown in Fig. 2, at each of the three wavelengths and for DPAR separately (Figs. 8 to 11).
At 380 nm, the mean values of the model outputs are overall higher than the observations. Apart from the seasonal cycle, a west-east gradient is also observed, with high values in the Eastern Mediterranean, most likely due to the low cloud cover (Fig. 8a, b). Almost all model outputs reveal a positive bias, with the largest differences in the southwestern Mediterranean (SWM2) and southern Adriatic seas (ADR2; up to 0.8 in December), as shown in Fig. 8c. The RMSD exhibits high values in the winter months in the majority of the sub-basins and the lowest values in summer in the Levantine region (LEV1, LEV2 and LEV3; as shown in Fig. 8d).
Similar to the findings at the BOUSSOLE site, the model attains low skills with respect to the BGC-Argo float observations at 412 nm, with a less pronounced west-east gradi- Figure 7. RMSD (left column) and slope (central column) and both indicators (right column) of the relationship between the OASIM and BOUSSOLE downward planar irradiance values (E d (λ0 − )) at the nine wavelengths (λ, in nm) and DPAR. The vertical bars in the left panels of each column are the BOUSSOLE (grey) and model (green) annual averages. The RMSD (W m −2 nm −1 ) and the slope (dimensionless) are marked in red and black, respectively, with the lines indicating the annual averages and the dots indicating the monthly averages. The monthly statistics are further separated by filtering out the day-by-day variability and are shown as the cyan dots for the RMSD and fuchsia dots for the slope. The averages over the time series are annotated in the panels. ent than that indicated by the observations (Fig. 9a, b). Except for the winter months in the southwestern Mediterranean (SWM1), Ionian (ION2) and Levantine (LEV1, LEV2, LEV3) sub-basins, the bias is negative overall (up to −0.4, as shown in Fig. 9c). The RMSD does not follow a clear pattern, but it is generally close to 0.3, with the highest values in December in the Ionian Sea (ION2), as shown in Fig. 9d.
From late spring to autumn, the bias is predominantly negative from the Tyrrhenian Sea eastwards, varying between 0 and −0.2 (Fig. 10c).
At 490 nm, the model values are higher during winter and early spring months, especially in the western sub-basins (Fig. 10a, b). The highest RMSD values are observed in the southwestern Mediterranean (SWM1) in November and December and in the Levantine region in February (Fig. 10d).
The modelled DPAR values are generally higher than the in situ observations (Fig. 11a, b), with the highest bias and RMSD values observed during the winter period (Fig. 11c,  d).

Summary of OASIM skills in the Mediterranean Sea
The absolute bias of the comparison of the OASIM outputs to the BOUSSOLE data (Table 1, upper part) is generally lower than 0.1 W m −2 nm −1 , and the regression slope approaches 1. In particular, the model vs. the observations reveals a total bias below −20 % of the average measured signal, and on average, the RMSD is approximately 30 %-40 %. The regression slopes are lower than 1, indicating a slight underestima- A lower agreement than that between OASIM and BOUS-SOLE is observed between OASIM and the BGC-Argo floats (Table 1, lower part), especially in terms of DPAR. Moreover, in this case, the bias is not always negative, being positive only at the 380 nm wavelength and for DPAR. The correlation (R) for E d λ, 0 − indicates a general good agreement between the model and BGC-Argo data (R ∼ 0.8), but the slope is below 0.8. Moreover, similar to OASIM vs. BOUSSOLE, E d (λ, 0 − ) reveals the highest discrepancy at 412 nm, with an RMSD of 0.3 W m −2 nm −1 , a bias of −0.13 W m −2 nm −1 and a regression slope of 0.51. The comparison to DPAR indicates a positive bias higher than 300 µmol quanta m −2 s −1 , which is high with respect to the BGC-Argo floats than in the BOUSSOLE case. A similar summary analysis performed for E d λ, 0 − was also performed for E d λ, 0 + to estimate the impact of reflection processes at sea atmosphere interface on irradiance (not shown). These processes are regulated by atmospheric pa-rameters shown in Fig. 4: percentual skill metrics indicate that RMSD is only marginally affected, with differences lower than 1 %, while bias for E d λ, 0 + shows differences lower than 5 % with respect to E d λ, 0 − .

Comparison of OASIM, BOUSSOLE and BGC-Argo floats in the northwestern Mediterranean Sea
Intercomparison of the model outputs and data obtained from the radiometric sensors of the BOUSSOLE buoy and BGC-Argo floats was possible only when the latter were located in the vicinity of the BOUSSOLE buoy in the NW Mediterranean Sea. Different spatial aggregations of profiles surrounding the fixed buoy were tested, ranging from 1 • (±0.5 • from the location of the buoy, as shown in Fig. 12) to the whole northwestern Mediterranean sub-basin (NWM), as shown in Fig. 13. In regard to the 1 • aggregation, up to 10 BGC-Argo profiles were available per month (Fig. 12, number not shown), while for the sub-basin analysis, the number of available profiles ranged from 8 in October to up to 100 in March (Fig. 13, number not shown). Monthly climatologies were calculated from 2004-2012 for BOUSSOLE and from 2012-2017 for BGC-Argo, with the model values corresponding to the locations of the instruments limited to the model spatial resolution. The BGC-Argo data were extrapolated to the surface, and the considered profiles followed the required conditions, as described in the previous section. The mean values and standard deviations are therefore shown for four different sets of results (Figs. 12 and 13, respectively): the data from the BOUSSOLE buoy and BGC-Argo floats with the corresponding OASIM outputs at 412 and 490 nm (expressed in W m −2 nm −1 ) and DPAR (expressed in µmol quanta m −2 s −1 ).
A separate comparison of the two different data sources conveyed an overall good agreement, with higher standard deviation values for the floats (Figs. 12 and 13) than those for the models, revealing a high variability range for the BGC-Argo float data.
The float values that matched with their corresponding OASIM outputs indicated high model outputs at 412 nm, with the largest differences during the summer months, consistent with Fig. 9, with a positive bias of up to 0.2 W m −2 nm −1 . The highest bias at the BOUSSOLE site was observed during spring with a similar magnitude (Fig. 13).
The best agreement was generally reached at 490 nm, as is also observed from the NWM column shown in Fig. 10c, especially during the winter months, where the differences between the two assessed platforms decreased to less than 0.1 W m −2 nm −1 (Fig. 13). The largest discrepancy was observed during summer, with the BGC-Argo floats resulting in higher values than those of the BOUSSOLE buoy (up to 0.2 W m −2 nm −1 , as shown in Fig. 13).
Consistent with the results shown in Fig. 11, major discrepancies arose when comparing DPAR, where the model values resulted in values much higher than those obtained from the floats, increasing especially during summer (up to 600 µmol quanta m −2 s −1 in August, as shown in Fig. 13).
Such inter-and intra-comparisons of the atmospheric radiative transfer model and available radiometric measurement platforms could also serve as a useful tool to estimate the range of variability when considering optical data from different sources. The spatial aggregation of measurements to the sub-basin level (i.e. a range of up to 10 • ) reveals the preservation of the irradiance seasonal variations. However, much remains to be explained in terms of the sources of both variabilities, both between the model and float data, espe- cially when considering DPAR, as well as between the two different sources, such as the floats and fixed buoy.

Discussion
An extensive comparison of OASIM's temporal and spatial variabilities was performed using the BOUSSOLE and BGC-Argo optical sensor data as references. The results indicated that, in general, the model reproduced the variability in the spectral downward irradiance in the Mediterranean Sea, which depends on the spatiotemporal scale. In particular, considering the OASIM applications within biogeochemical models equipped with a multispectral in-water optical module, the impact will be different according to the specific scale under investigation. In terms of fine temporal scales, it was observed that a large part of the discrepancy expressed by the RMSD occurred due to the day-to-day variability, similar to the findings reported by Somayajula et al. (2018). When this temporal scale is filtered out, the RMSD decreases by half or more. Consistently, the assessed uncertainty (an RMSD of approximately 10 %, bias lower than 10 % and slope > 85 %) should be considered in multiannual simulations. In fact, a key parameter such as the primary productivity (strongly affected by the irradiance) exhibits a dominant component of the variance at the seasonal scale (Lazzari et al., 2012), while the day-to-day variability and interannual variability seem to be less important (Di Biagio et al. 2019). In contrast, for short-term forecasts (i.e. 10 d lead time; Salon et al., 2019), the day-to-day variability in the downward irradiance could be relevant. Clearly, to fully investigate the high-frequency RMSD variability, simulations should be refined by increasing the spatial resolution of the model in the region surrounding the BOUSSOLE site (e.g. C3S, 2017;Hersbach et al., 2020).
Apart from the diurnal cycle, the cloud cover and aerosols are major drivers in modulating the variability in the downward irradiance, with aerosols being subordinate to the cloud cover. These two input parameters show a high seasonal variability (Figs. 14 and 15). In the Mediterranean area, the cloud cover is high during fall and winter, while the AOT is high during spring and summer (Fig. 5). The winter maximum of the cloud cover corresponds to the maximum RMSD at all wavelengths considered (Fig. 14), while the minimum RMSD occurs in July when the cloud cover is also at its annual minimum. A consistent variability is obtained by computing the regression slope (Figs. 14 and 15). These results are in line with a previous multimodel analysis (Somayajula et al., 2018) indicating that models present a negative bias when the cloud cover is higher than 70 %. Such an underestimation could impact phytoplankton dynamics modelling; thus further analysis should be conducted. The multimodel comparison by Nielsen et al. (2014) revealed that the Slingo liquid-optics model tends to overestimate the cloud optical thickness and therefore underestimates the irradiance.
The presence of a systematic underestimation, following the mechanism described in Nielsen et al. (2014), likely affects all temporal scales. However, in the present study, it was shown that the RMSD and regression slope are greatly improved when filtering out the day-to-day variability. This indicates that high-resolution sampling in OASIM could notably improve the model results. In other words, this implies that reducing the spatiotemporal uncertainty in cloud cover is more important than specific details of the liquid cloud optics parameterizations.
The cloud cover spatial distribution indicates that similar conclusions are obtained based on the model vs. BGC-Argo skills with the results obtained at the BOUSSOLE site; namely, at least for E d (λ = 490, 0 − ) and E d (λ = 380, 0 − ), the increase in the match between the model and data is modulated by the cloud cover (Figs. 5a, 8, 9 and 10). In contrast, E d (λ = 412, 0 − ) seems to be less influenced by clouds with the exception of the minimum RMSD during summer in the eastern area.
To investigate the balance between direct and diffuse irradiance components, we introduced a further diagnostic (IND) based on their fraction, defined as where IND is non-dimensional and varies in the in- indicates perfect balance between the two components: E dir λ, 0 − = E dif λ, 0 − . As show in Fig. 14, IND provides a similar interpretation of cloud cover; in fact the model skill is higher when IND is higher and vice versa. This diagnostic indicator could be used to generalize results in regions outside the Mediterranean Sea characterized by lower sun angles, implying a different balance of the direct versus the diffuse component. In these situations, the effect of clouds in increasing RMSD and bias could be even higher. It is worth mentioning that in the present case, since IND covariates with cloud cover, it is difficult to separate the role of clouds from direct versus diffuse irradiance ratio. Nonetheless, IND at 412 nm is lower than all the other wavelengths, and this could explain, at least in part, the lower skill observed at 412 nm. Table 1. Summary of the model skill compared to the available data from the BOUSSOLE buoy (from 2004 to 2012) and BGC-Argo floats (from 2012 to 2017) for the irradiance (E d ) at the different wavelengths (WLs) and for DPAR. RMSD, bias and y-int are expressed in watts per square metre per nanometre, and in regard to DPAR, the same statistical indicators are expressed in µmol quanta m −2 s −1 , while all the other indicators (regression r and slope) are dimensionless, where N is the number of match-ups between the model and observations. For the BOUSSOLE comparison, the numbers in bold font are derived by filtering out the day-to-day variability (i.e. the intra-monthly variability). Given the large number of samples, all statistics are significant (p value < 0.05). For the RMSD and bias, the percentage values normalized by average data are reported in parentheses. Apart from cloud dynamics, aerosols play a role in landlocked regions such as the Mediterranean Sea (Papadimas et al., 2008;Nabat et al., 2015). The AOT reveals a variability at least an order of magnitude larger than that of the asymmetry parameter or even more than that of the single-scattering albedo. The minimum RMSD observed in July does not correspond to an extreme AOT value, as observed in the case of the cloud cover. The AOT generally decreases with increasing wavelength (Fig. 15). However, a reduced model skill is observed at 412 and 682.25 nm. The MODIS aerosol data are extrapolated from the seven wavelength channels of the satellite sensors (i.e. 470, 550 and 660 nm in the visible region), which could possibly explain the high uncertainty at 412 nm. In terms of the spatial heterogeneity, comparing the model and BGC-Argo floats (Figs. 5b and 8,9,and 10), the role of aerosols in the modulation of the model bias and RMSD does not appear to have a clear interpretation.
The low skill at a specific wavelength could be also explained in terms of the wavelength discretization in OASIM: the current model spectral resolution of 25 nm could be refined near 412 and 682.25 nm to investigate whether this would reduce the bias, especially due to the fact that in the present simulations the 412 nm wavelength is at the interface between the band centred at 400 nm and the one centred at 425 nm.
In addition to the seasonal indicators discussed above, the interannual trends of E d (λ, 0 − ) and DPAR were investigated both for the measured data and model results. Given the strong seasonal cycle present in all the properties considered, low-pass filters (i.e. moving averages) were applied to the data before the regression.
The results for both datasets were biased by the gaps in the acquisitions, which introduced spurious trends. Therefore, the focus was on the model inputs instead (ECMWF products and MODIS aerosol data) and the outputs corresponding to a 14-year gap-free time series both in spatial and temporal terms. The model outputs averaged over the Mediterranean Basin indicated a low interannual variability, or at least the procedures adopted based on the low-pass filtering and sub-

Conclusions
We assessed the performance of the OASIM radiative model on the Mediterranean Sea for the period 2004-2017, comparing the model outputs with multiplatform reference data. The BOUSSOLE mooring station provides a dataset with a high temporal resolution at a fixed point, while the BGC-Argo floats allow the partial resolution of the spatial vari- ability (certain regions contain no floats), although with a lower spectral resolution than that of BOUSSOLE. The atmospheric multispectral input data provided by OASIM are necessary to resolve the multispectral propagation of light within the water column. Evaluating the uncertainties and the quality of the these input data is fundamental for all the future applications involving bio-optical modelling and constitutes an important starting point to develop assimilation schemes based on bio-optical modelling. The results indicate an overall good agreement between the model outputs and in situ references, highlighting a clear seasonal variability in the model capabilities, with a generally high performance during summer.
The analysis conveys that an important source of discrepancy between the model and data is the intra-monthly variability, which can be ascribed to cloud cover dynamics (the RMSD is reduced with a low cloud cover). The coarse resolution of cloud cover likely affects the model skill at both Figure 15. Comparison of the OASIM downward irradiance (E d (λ0 − ) at the nine wavelengths) to the BOUSSOLE data from 2004 to 2012 in terms of the RMSD and regression slope and their relationship with the MODIS aerosol optical thickness. The left section of each panel shows the monthly climatology of the RMSD (normalized by its averaged value, red lines and labels) and regression slope (normalized by its averaged value, black line), superimposed on the monthly climatology of the aerosol optical thickness (green bars and labels). Regression slope thresholds at 1 (dotted black line) and 0.75 (dot-dashed black line) are also shown. The right section of each panel shows the monthly means of the time series of the regression slope (black dots, with the 1 and 0.75 thresholds shown; dotted and dot-dashed black lines, respectively), superimposed on the monthly means of the time series of the aerosol optical thickness (green bars and labels). seasonal and intra-monthly timescales. This is shown by the improvement in model skill at BOUSSOLE when variability between days is filtered out.
Nonetheless, the analysis of the relative contribution of E dir λ, 0 − and E dif λ, 0 − indicates that the model skill is correlated to their ratio, suggesting that improving the physical description of the radiative processes should be considered. To this end, novel atmospheric models with improved mathematical descriptions of cloud dynamics and advanced numerical solvers may better simulate both clear-sky and cloudy-sky components (Hogan and Bozzo, 2018).
In regard to climatological simulations, the high skill in terms of the monthly averaged irradiance is probably sufficient to properly constrain biogeochemical dynamics, whilst attention should be paid in the case of short-term simulations, when biogeochemical processes such as chlorophyll acclimation exhibit the same timescale as the relatively large fluctuations observed in the RMSD.
In this case, in addition to improved cloud parameterizations, the use of daily resolved aerosol data could possibly reduce the model uncertainties, for example, the EU-METSAT polar multisensor aerosol optical property product (PMAP), available since 2014, or other products provided by the Copernicus Atmosphere Monitoring Services (CAMS). For recent applications see also Bozzo et al. (2020) and Rontu et al. (2020).
Among the wavelengths considered for the downward planar irradiance, 412 and 681 nm appear to result in large uncertainties. However, at this stage, it is not possible to ascertain the reasons for such different skills, and further investigations are required. In the present application in the Mediterranean Sea, certain input variables are not directly available from the ERA-Interim dataset in the form required by OASIM, so further processing is needed, as detailed in Table B1. Table B1. Correspondence between the input variables required in OASIM (left column) and the implementation in the Mediterranean Sea (right column), with specific pre-processing steps.

OASIM variable Present implementation
Cloud  Gregg and Carder (1990), the precipitable water [cm] is computed as a function of the saturated water vapour pressure, surface pressure and sea level atmospheric pressure (Garrison and Adler, 1990) Code availability. The OASIM Fortran code is publicly downloadable, together with the reference input data, from http://gmao.gsfc. nasa.gov/research/oceanbiology (Gregg, 2020).
Author contributions. PL and SS designed the study and the numerical experiments, and PL carried them out supported by WWG. PL and ET performed the model data comparison. VV was in charge of the BOUSSOLE mooring buoy maintenance and performed the quality control of the optical data. DA provided data and funding through the BOUSSOLE project. EO compiled the BGC-Argo database and performed the quality control of the optical data. FDO provided funding and contributed with float deployments. All authors contributed to the interpretation of the results. PL prepared the paper with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. The ERA-Interim dataset are published under a Creative Commons Attribution 4.0 International (CC BY 4.0), https: //creativecommons.org/licenses/by/4.0/ (last access: 12 November 2020). ECMWF does not accept any liability whatsoever for any error or omission in the data, their availability, or any loss or damage arising from their use.
Special issue statement. This article is part of the special issue "Biogeochemistry in the BGC-Argo era: from process studies to ecosystem forecasts (BG/OS inter-journal SI)". It is a result of the Ocean Sciences Meeting, San Diego, United States, 16-21 February 2020.