Empirical Relationships between Remote-Sensing Reﬂectance and Selected Inherent Optical Properties in Nordic Sea Surface Waters for the MODIS and OLCI Ocean Colour Sensors

: The Nordic Seas and the Fram Strait regions are a melting pot of a number of water masses characterized by distinct optical water properties. The warm Atlantic Waters transported from the south and the Arctic Waters from the north, combined with the melt waters contributing to the Polar Waters, mediate the dynamic changes of the year-to-year large-scale circulation patterns in the area, which often form complex frontal zones. In the last decade, moreover, a signiﬁcant shift in phytoplankton phenology in the area has been observed, with a certain northward expansion of temperate phytoplankton communities into the Arctic Ocean which could lead to a deterioration in the performance of remote sensing algorithms. In this research, we exploited the capability of the satellite sensors to monitor those inter-annual changes at basin scales. We propose locally adjusted algorithms for retrieving chlorophyll a concentrations Chla , absorption by particles a p at 443 and 670 nm, and total absorption a tot at 443 and 670 nm developed on the basis of intensive ﬁeld work conducted in 2013–2015. Measured in situ hyper spectral remote sensing reﬂectance has been used to reconstruct the MODIS and OLCI spectral channels for which the proposed algorithms have been adapted. We obtained MNB ≤ 0.5% for a p ( 670 ) and ≤ 3% for a tot ( 670 ) and Chla . RMS was ≤ 30% for most of the retrieved optical water properties except a p ( 443 ) and Chla . The mean monthly mosaics of a p ( 443 ) computed on the basis of the proposed algorithm were used for reconstructing the spatial and temporal changes of the phytoplankton biomass in 2013–2015. The results corresponded very well with in situ measurements.


Introduction
Climate change processes are the most pronounced in the Arctic, where the rate of temperature rise is the strongest [1]. Increasing temperatures lead to the reduction and thinning of the sea ice cover [2][3][4]. During last 10 years, the sea ice extents were the lowest on record [5]. The reduction in sea ice cover area and sea ice thickness is exposing more open water areas to the absorption of light [6] and the dissipation of this absorbed energy into heat. This positive feedback (ice-albedo feedback) may influence upper ocean layer thermodynamics and accelerate warming [7].
The amount of energy absorbed by oceanic water depends solely on its inherent optical properties (IOPs): these are additive functions of the spectral IOPs of clear water alone and the spectral The cold East Greenland Current (EGC) and the Sørkapp Current (SC) marked in blue, and the warm West Spitsbergen Current (WSC) marked in red. Match-ups with chlorophyll a concentration maps from MODIS are depicted by red asterisks (modified [19]).
The ocean optics theory linked IOPs with Apparent Optical Properties (AOPs), and the most common approximate relationship between the remote-sensing reflectance and IOPs derived by Morel and Prieur (1977) [20], was the theoretical foundation of ocean colour remote sensing. Ocean colour remote sensing has been used to monitor IOP dynamics at the ocean surface for more than three decades. From the technical point of view, the quality of ocean colour measurement depends on three main factors: (i) the number of spectral bands, which restricts the level of detail of a distinguishable reflectance spectrum (ii) the positions of the bands, which should focus on characteristic spectral signatures of the reflectance spectra, and (iii) the high-quality vicarious calibrations that ensure instrument stability over time. Numerous sensors dedicated to ocean observation have evolved during the last three decades. For example, the Sea-viewing Wide Field-of-View Sensor (SeaWiFS), placed in orbit in 1997, was followed by the Moderate-Resolution Imaging Spectroradiometer (MODIS) deployed on the AQUA satellite, launched in 2002 and operated by US NASA [21]. In 2002, the European Space Agency (ESA) also launched an Earth observation satellite. The Medium-Resolution Imaging Spectrometer (MERIS), in operation until 2012, was highly successful. The descendant of MERIS, the Ocean and Land Colour Instrument (OLCI), was developed within the framework of Copernicus, the European Union's Earth monitoring programme. The OLCI radiometer is equipped with a wider range of spectral bands, enabling more subtle changes in remote sensing reflectance spectra to be captured. Recent modifications have provided flexibility in defining the spectral bands' placement before launch, and it is now possible to adjust the spectral bands in orbit [22].
The general concepts of the ocean colour algorithm for retrieving the optical properties of surface ocean waters and various biogeochemical variables based on bio-optical models have changed little during the last thirty years [23,24]. They involve three general approaches: (i) empirical band-ratio algorithms based on non-linear relationship between the chlorophyll a concentration and the ratio of the remote sensing reflectances measured in the blue and green spectral bands [25], (ii) semi-analytical algorithms based on the relationship between IOPs and biogeochemical variables estimated from bio-optical models [26,27], and (iii) neural network algorithms [28]. The widely used standard global band ratio OC2-OC3 algorithms [29] have been employed for the accurate retrieval of chlorophyll a levels even in optically complex waters [22]. The technological progress achieved in satellite sensor development, along with the changes to phytoplankton ecology forced by climate change, have stimulated the development of satellite algorithms and led to the revision of existing algorithms. This has been followed by calibration and imagery product validation campaigns [30]. The readjustment of locally valid algorithms is particularly desirable in specific regions, such as high latitude and polar regions [31]. According to Stramska et al. (2003) [32], the NASA global OC and Chlor-a-MODIS algorithms generally overestimated chlorophyll a concentrations in the Nordic Seas by a factor of nearly 2 at low chlorophyll a concentrations (<0.2 mg m −3 ) and underestimated it at higher concentrations (20-50% at 2-3 mg m −3 ). Stramska et al. (2006) [11] also found significant variations in IOPs driven by changes in the photosynthetic pigment composition in different phytoplankton communities. Both studies were based on in situ observations conducted between 1998 and 2003, before the phytoplankton phenology shift took place [11,32].
The principle aims of the current study, based on intensive field campaigns during three consecutive summer seasons from 2013 to 2015 were: (i) to re-adjust the empirical band ratio algorithm for retrieving spectral values of the total absorption coefficient, particulate absorption coefficient and chlorophyll a concentration to matching MODIS and OLCI spectral bands, (ii) to test algorithm uncertainty in the light of recent changes in phytoplankton bloom characteristics and to compare this with routine global OC2-OC3 algorithm performance in the Nordic Seas, and (iii) to determine the major source of uncertainties in the remote sensing optical water constituents.

Study Area
The field studies were conducted in July 2013, 2014 and 2015 in Fram Strait and the Nordic Seas, mainly in the AW-dominated eastern part of Fram Strait and along a transect between northern Norway and southern Spitsbergen. The location of the sampling station is shown in Figure 1, and the field survey dates and sampled parameters are listed in Table 1. This area is one of the most important regions of the global ocean, where the highest proportion of its deep water is formed and a number of different water types mix [33,34]. The eastern part of Fram Strait, a 450 km wide passage between Greenland and the Svalbard Archipelago, is occupied by warm, saline Atlantic Water (AW) that enters the Nordic Seas, transported northwards along the Spitsbergen coast by the West Spitsbergen Current (WSC). The western part of Fram Strait is under the influence of Arctic Water, carried by the East Greenland Current (EGC). Glacial meltwater combined with sea ice melt and water from fjords form a layer of distinct Polar Waters (PW), which overlies the more saline and dense AW [35]. There are strong horizontal and vertical gradients of water properties across the current boundaries that form strong frontal structures visible on satellite sea surface temperature (SST) or chlorophyll a concentration imagery. Heat transported with AW controls water stability in the area. Atlantic waters have accelerated ice melt [36], which influences the freshwater budget, nutrient and carbon biogechemical cycles that govern the primary production rate and composition of phytoplankton communities [35]. The transport of anomalously warm and saline waters into the Nordic Seas with the intensified North Atlantic Current has been recorded ever since 1995 [33,37]. Table 1. Number of surface measurements of R rs (λ), Chla, a p (λ) and a CDOM (λ) made during each of the three cruises that were taken into consideration for the formulation of the algorithm. The numbers of the datasets in Table 4 are indicated by an asterisk *, because the number of match-ups varied for the two atmospheric correction algorithms.

Sample Collection and Processing
Water samples for determining chlorophyll a concentration, CDOM and particulate absorption were collected with a SeaBird SBE32 Carousel Water Sampler equipped with Niskin bottles, a SBE 911 plus CTD (SBE 9plus CTD Unit and the SBE 11plus Deck Unit) and a Wetlabs ECO chlorophyll fluorometer FLRTD-1547. Samples were collected at three depths: near the surface, at ca 2 m depth; at the chlorophyll a maximum, usually between 15 and 25 m; and below the chlorophyll a maximum, between 30 m and 50 m. The exact position of the chlorophyll a maximum depth was estimated from the vertical profile of chlorophyll a fluorescence visualized during a CTD downcast. The water samples were filtered immediately after collection under a low vacuum. All pigment and particulate absorption samples were passed through Whatman (GE Healthcare, Little Chalfont, UK) 25 mm GF/F filters. The filter pads with retained particulate material, were immediately deep frozen in a freezer and thereafter stored at −80 • C prior to analysis.
In 2013 (AREX2013), the samples for CDOM absorption measurements were immediately filtered in two steps: first through acid-washed GF/F filters, then through acid-washed Sartorius 0.2 µm pore size cellulose membrane filters to remove finer particles. From 2014 onwards, CDOM samples were filtered directly from rosette Niskin bottle taps with an Opticap XL4 Durapore flow through filter cartridge of nominal pore size 0.2 µm into acid washed amber glass 200 mL sample bottles. The cartridge filter was kept in 10% HCl solution and was rinsed with ultrapure MilliQ and sample water before CDOM sample collection. In 2013 and 2015, the unpreserved water samples for determination of CDOM absorption were stored in dark, at temperature of 4 • C, and were transferred after the cruise to a land-based laboratory for the spectroscopic measurements. In 2014, all the spectroscopic measurements for determining CDOM absorption were performed in the laboratory onboard r/v Oceania, immediately after sample collection, using a double-beam Perkin Elmer Lambda 35 spectrophotometer in the 240-700 nm spectral range in 10 cm quartz cuvettes. CDOM absorption of samples collected during the AREX2013 and AREX2015 cruises was measured with a double beam Perkin Elmer Lambda 650 spectrophotometer in the same spectral range and the same cuvettes in the laboratory at the Institute of Oceanology, Polish Academy of Sciences, Sopot, Poland (IOPAN). Optical densities (absorbance) OD(λ) were converted into CDOM absorptions by multiplying OD(λ) by 2.303 and then dividing by the optical pathlength (m). Finally, the spectra were corrected for residual scattering as per Kowalczuk et al. (2006) [38]. The CDOM absorption spectrum slope coefficient in the 300-600 nm spectral range was calculated by nonlinear least squares fitting that employed the Trust-Region algorithm implemented in Matlab R2011b [38,39].

Particulate Absorption
After the cruises, particulate absorption samples were transferred in liquid nitrogen to IOPAN. There, spectra of the light absorption coefficient by suspended particles a p (λ) were measured in the 300-750 nm spectral range with a UNICAM UV4-100 double-beam spectrophotometer equipped with an integrating sphere (LABSPHERE RSA-UC-40). The transmission-reflection (T-R) method described by Tassan and Ferrari [40,41] was applied to the a p (λ) measurements with further recommendations by Doerffer et al. (2002) [42] and Mueller et al. (2003) [43]. These techniques allow one to measure the optical density OD(λ) of the suspended material collected on the filter: where A(λ) is the total absorbance obtained from measurements in transmission and reflectance modes. Blank samples were measured for every batch of samples to estimate the correction factors for filter pad light scattering and absorption. The clean filter pads used for the blank reference measurements were rinsed with 30 mL filtered sea water. The path length amplification correction had to be applied to the calculation of particulate absorption coefficients from the OD s (λ) spectra. Light scattering by the suspended matter collected on the filter elongated the optical path length of photons passing through the filters. The correction was applied using the dimensionless path length amplification β-factor, which converts the optical density of particles collected on the filter (OD s (λ)) into the optical density of particles in solution (OD sus (λ)) [44]. The whole data set was corrected with the β-factor proposed by Stramski et al. (2015) [45] for the T-R method: OD sus (λ) = 0.496 · OD s (λ) 2 + 0.388 · OD s (λ) The absorption coefficient of particles was calculated using the formula: where l [m] is the hypothetical optical path in solution, determined as the ratio of the volume of filtered water to the effective diameter of the filter, i.e., the diameter of the coloured area on the filter measured individually for each sample. The coefficient of absorption by non-algal pigments a N AP (λ) was determined from the results of analogous measurements after bleaching off the algal pigments for 20 min with 2% calcium hypochlorite Ca(ClO) 2 [46,47]. The filters were rinsed with artificial sea water (64 g NaCl in 1 L of MiliQ water). The algal pigment absorption coefficient a ph (λ) was determined as the difference between a p (λ) and a N AP (λ) using the value of a ph (λ) at wavelength 750 nm as the null-point correction [11,32,48]. Although the null-point correction at 750 nm of a ph (λ) spectra is recommended in ocean optics data protocols [48], the non-negligible a ph (λ) absorption in the 750-800 nm spectral range, found by Tassan and Ferrari (2003) [41], could be an inestimable source of error in our data set [49]. The total absorption coefficient a tot (λ) is the total non-water absorption coefficient of seawater a nw (λ) (the light absorption coefficient after subtraction of the pure water contribution), which was computed by adding a CDOM (λ) and a p (λ).

Chlorophyll a Concentration
Defrosted samples for spectrophotometric measurements of the chlorophyll a concentration were extracted in 96% ethanol for 24 h at room temperature. The extract was centrifuged for 15 min and decanted. The chlorophyll a concentration Chla was determined spectrophotometrically [50] using Unicam UV4-100 (2013) and Perkin Elmer Lambda 650 (2014 and 2015) instruments. The optical density OD(λ) of the pigment extract in ethanol was measured in a 2 cm cuvette. The raw OD readings at 665 nm were corrected for the background signal in the near-infrared region (750 nm). Chla was calculated from the following formula (Strickland and Parsons (1972)) [51]: where e is the volume of ethanol in cm 3 , OD(665) is the absorbance at 665 nm (after correction for blank ethanol), OD(750) is the analogous absorbance at 750 nm, 83 [L g −1 cm −1 ] is the chlorophyll a-specific absorption coefficient in ethanol, l is the length of the cuvette in cm and V is the volume of filtered seawater in litres.

Radiometric Measurements
The radiometric measurements were synchronized with the water sample collection. All radiometric acquisition and data processing routines were performed in compliance with NASA 2003 protocol recommendations [43]. Two hyperspectral TriOS RAMSES radiometers were used, operating in the 320-950 nm spectral range and collecting information in 194 wavebands, each 3.3 nm wide. A RAMSES MRC-VIS radiometer, mounted on a frame floating at distance of 15-20 m from the ship, was used to measure the upwelling radiance L u (λ) just below the sea surface [52]. A RAMSES ACC-VIS sensor, mounted on the top of the A-frame winch in the ship's stern, well away from the superstructure, simultaneously recorded the incident solar irradiance E s (λ). The remote-sensing reflectance R rs (λ) was computed from the radiometer signals using the formula: where E s (λ) is the downwelling incident solar irradiance measured above the sea surface [W m −2 nm −1 ], which is equivalent to is the upwelling irradiance just below the sea surface, ρ ∼ = 0.025 is the Fresnel reflectance of the air-sea interface, and n ∼ = 1.34 is the refractive index of seawater, so that 1−ρ n 2 is the upward radiance transmittance of the sea surface for normal incidence from below and may be approximated as 0.543 [43,53].
Reflectance spectra were measured in situ at varying wavelengths by both TriOS RAMSES radiometers. The remote sensing reflectance R rs (λ) was unified using a local polynomial regression fitting method, implemented in R as the loess function [54]. MODIS AQUA and OLCI spectral bands were reconstructed on the basis of the Relative Spectral Response function (RSR) (see Figure 2) provided for the MODIS sensor and for each of the satellite sensors, then convoluted based on the following formula: where R k is the reflectance sensed within a satellite channel k, λ min and λ max are the respective lower and upper bounds of the RSF band k filtration function, r i is the RSR of band k at wavelength i, and p i is the interpolated spectral response of RAMSES at wavelength i. The band-weighted reflectances of RAMSES were derived for the MODIS AQUA [55,56] and OLCI sensors [57,58].

Satellite Data Processing
The satellite data subsequently used further for product validation were obtained from the Moderate Resolution Imaging Spectroradiometer mounted on-board the AQUA satellite and acquired from the GSFC NASA service http://oceancolor.gsfc.nasa.gov. All the images were processed from level-0 to level-2 using SeaWiFS Data Analysis System (SeaDAS) software, version 7.3, the default being for the MODIS algorithm OC3M for Chla retrieval (l2gen) [59]. Two types of atmospheric corrections were applied in order to match the MODIS scenes: (i) one based on the dark pixel correction [60] using Near Infrared (NIR) bands as references (748 and 869 nm) [61], (ii) the other discriminating case-1 waters based on the turbid water index T ind (748, 1240) and applying standard NIR correction (748, 869) or Shortwave Infrared (SWIR) bands for turbid waters (1240, 2130) [62].

Statistical Approach
The model development and validation were performed in accordance with well known practices [63,64]. The set of statistical metrics were used for non-parametric quality assessment, based on the recommendation of Greg and Casey (2004) [65]. The mean normalized bias (MNB) (systematic error) as well as the normalized root mean square (RMS) error (random error) were used after Darecki and Stramski (2004) [66]: We also used statistics based on the logarithm of the ratio of the algorithm-derived y i to the measured valueŷ i , which are standard in the ocean colour literature (e.g., [59,66,67]): The fitted models were assessed using the coefficient of determination R 2 and the Pearson correlation coefficient R, computed between predicted and measured in situ values.

Results
The empirical formulas for retrieving the inherent optical properties of sea water from the remote sensing reflectance R rs (λ) were based on the in-water measurements of this parameter. The remote sensing reflectance and the total absorption coefficient spectra measured during the field campaigns in summers 2014 and 2015 over the West Spitsbergen Shelf, with colour scaled according to increasing Chla, are illustrated in Figure 3. The R rs (λ) spectra are U-shaped, characteristic of mesotrophic oceanic waters [20]. In the blue part of the spectrum, there is a superposition of several processes such as the CDOM absorption increasing towards the shortwave, scattering of suspended particles and absorption by phytoplankton pigments with a maximum at 443 nm. When Chla and CDOM input is low, we observe almost flat spectrum in the range 400-500 nm with lower R rs (λ) values only below 400 nm. When Chla rises, a trough appears at 443 nm accompanied by two local maxima around 400-412 nm and 530-550 nm which is clearly visible on the spectra measured in 2015 at Chla >8 mg m −3 . Another distinctive characteristic of R rs (λ) is the third maximum at 670 nm due to chlorophyll a fluorescence emission. The intensity of this local maximum is related to Chla: elevated R rs (670) values were indeed recorded in summers 2014 and 2015, when phytoplankton biomasses were greater than in 2013 [10]. The R rs (670) maximum could also be used as a diagnostic tool to distinguish between remote sensing reflectance spectra that were most likely shaped by factors other than absorption and scattering due to autotrophic protist growth. The blue line indicates the R rs (λ) spectra measured in Chla depleted surface waters, usually occurring very close to the marginal sea ice boundary. The sea ice melt water in this frontal zone usually being deficient in Chla [10,68], the peak at R rs (670) was very small if it was at all visible. In contrast, the suspended particles released from the melting ice contributed to a higher level scattering, as indicated by the increase in R rs (λ) in the 400-500 nm spectral range (Figure 3a-c). The spectra of the total absorption coefficient a tot (λ) were computed as a sum of a CDOM (λ) and a p (λ). We omitted the contribution of pure sea water absorption, a w (λ) to the total absorption coefficient, because a w (λ) is almost constant in the ocean and does not contribute to a tot (λ) variability.
Similarly to the measured R rs (λ) spectra, a significant year-to-year variability in a tot (λ) spectral shape and magnitude were observed. In 2013, when Chla were low, the absorption spectra followed exponential a CDOM (λ) shape (Figure 3d). Elevated Chla in the summers 2014 and 2015 caused higher total absorption coefficient values near 440 nm and 670 nm (Figure 3e,f). The detailed information on absorption properties by most significant optical constituents of oceanic waters and absorption budget in the main water masses in the study area during field campaign in 2013, 2014 and 2015 were extensively described in Kowalczuk et al., (2019) [10].
The close dependence between the spectral shape and magnitude of measured R rs (λ) on Chla motivated us to derive a set of empirical in-water, band ratio algorithms for retrieving chlorophyll a concentrations and spectral values of selected IOPs. To guarantee the statistical significance of these models, prior to regression analysis, we tested the normality of distribution of the optical parameters used to derive these algorithms. The normality of distribution of the in situ measurements was tested using the Shapiro-Wilk test: this was applied to each of the in situ datasets of Chla, a CDOM (443), a p (443), a p (670), a tot (443) and a tot (670). In most cases, the distribution was log-normal. The null hypothesis of the normal distribution of absorption by CDOM at 443 nm and its log-transformed form, had to be rejected, however, as the correlation between log-transformed a CDOM (443) and the remote sensing reflectance band ratio was higher than the linear data set; we therefore applied regression analysis to the log-transformed values of the respective parameters (Figures 4a and 6a). Our dataset contained almost 160 pairs of simultaneous measurements of inherent and apparent optical properties and fell within the recommended threshold number of samples required for algorithm development [69]. We randomly split the datasets into development and validation subsets. Owing to the varying numbers of measurements depending on the water property in question, the sizes of the subsets were also varied in order to maintain the requisite proportion between them. We decided to derive a band ratio algorithm, which yielded satisfactory results at global scale, especially in Case-1 waters [29]. In the first approach, we performed regression analysis to select the OLCI and AQUA spectral bands with the highest coefficient between combinations of band ratios and different water properties. Correlation matrices were computed for the log-transformed data, e.g., log Chla vs. log R rs (λ 1 )/R rs (λ 2 ). The simulation of the satellite radiometer spectral bands from hyperspectral in situ measurements ensured a sufficient match between the identified relationships and the particular satellite sensor.

Empirical Band Ratio Algorithms for MODIS Spectral Bands
The highest correlation coefficients were achieved for the regression between selected IOPs and the same combinations of band ratios for the simulated MODIS bands (B8-B15): log R rs (443)/R rs (531), log R rs (443)/R rs (547), log R rs (488)/R rs (531), log R rs (488)/R rs (547), log R rs (531)/R rs (547) (Figure 4). The highest correlation coefficients were calculated for the relationships between a p (443) and a p (670), and the respective band ratios, which oscillated around 0.9 (Figure 4c,d). The maximum correlation coefficients for the relationship between Chla and a tot (670) and the respective band ratios were close to 0.85 (Figure 4b,f).The relationships between a tot (443) and selected remote sensing reflectance band ratios were weaker than the IOPs mentioned above; the correlation coefficients reached a maximum of 0.8 (Figure 4e). The correlation was the weakest for a CDOM (443) and selected remote sensing reflectance band ratios, and the correlation coefficient did not exceed 0.3 (Figure 4a). The high correlation between a CDOM (443) combinations of band ratios in red and near infra-red (667 nm and 678 nm), as illustrated on Figure 4a, should be regarded as an artifact, beacuse of the negligible CDOM absorption in this spectral region and the cross-correlation effect between those neighbouring bands caused by the broad half width. Thus, no further attempts to fit regression lines between a CDOM (443) and log R rs (λ 1 )/R rs (λ 2 ) were made. As the maximum correlation coefficients were predominantly comparable for each parameter and varied only in the second or third place after the decimal point, we decided to choose the best formula based on two criteria, which was possible on splitting the data into two subsets. One indicator was the R 2 of the regression line fitted to the relationship between the parameter and the band ratio applied to the test subset. The other was the correlation coefficient between the values of a given parameter predicted with the algorithm developed on the basis of test data subset, and the in situ measurements contained within the validation data subset. Such an approach prevented from over fitting of the regression line ( Figure 5). Both statistics were computed for all the best correlating band ratios obtained from the matrices. For almost all the IOPs of water, the best spectral band ratio combination appeared to be R rs (488) vs. R rs (547). Even if the matrices suggested similar results for the two band ratios R rs (443)/R rs (547) and R rs (488)/R rs (547), the regression line fitting and data reconstruction obtained from validation dataset performed better, which was confirmed by repeating the calculations on new, randomly selected subsets. The total absorption at 443 nm a tot (443) was the only one to be retrieved with greater accuracy using the R rs (443) to R rs (547) ratio (Table 2).
Having established the optimal band pair, the order of the polynomial was selected from the first-to fourth-order range. The assessment was performed using six indicators. One of them, the coefficient of determination R 2 , referred to the test subset, the others to the validation subset, i.e., the Pearson correlation coefficient r, the mean normalized bias (MNB) (systematic error) (Equation (7)), the normalized root mean square error (RMS) (Equation (8)), the logarithmic form of the RMS and bias (Equations (9) and (10)), and the average mean squared error LOOCV obtained from the leave-one-out cross-validation method. The second-order polynomial was identified as optimal for all the parameters ( Table 2). The difference between the second-and third-order polynomial was small within the range of the subsets. However, outside the range of measured values, the risk of over-or underestimation increased rapidly owing to the curvature of the fitted line ( Figure 5).  Table 2. Proposed empirical second-order polynomial formulas for retrieving Chla, a p (443), a p (670), a tot (443) and a tot (670) from spectral reflectance ratios adjusted to MODIS AQUA bands.

Parameter Equation
Chla The uncertainty of the derived empirical second-order polynomial formulas for retrieving Chla, a p (443), a p (670), a tot (443) and a tot (670), based on MODIS AQUA spectral band ratio combinations, was calculated by comparing the values estimated using the in situ spectral reflectance and the spectral values of the IOPs contained in the validation data set. In general, we achieved a low level of uncertainty, with MNB ≤ 2.1% and RMS ≤ 27% for all the formulas proposed for retrieving spectral values of particulate and total absorption coefficients (Figure 5a-h). The correlation coefficient between the estimated and measured IOPs was r > 0.9, which also yielded high values of the determination coefficient R 2 > 0.73 ( Figure 5). The empirical formula proposed for Chla retrievals was characterized by a slightly larger uncertainty: although MNB was ≤2.1%, RMS was almost 35%, which is comparable with the accuracy obtained for the open ocean colour remote sensing methods used in oceanographic studies [65].

Empirical Band Ratio Algorithms for OLCI Spectral Bands
We reconstructed sixteen OLCI bands (1-16) from the TriOS radiometric measurements. Correlation matrices depicting correlations between the band-ratio combinations with the parameters revealed that most of the variables, i.e., Chla, a p (443), a p (670), a tot (443) and a tot (670), correlated well with the band ratios containing combinations of the channels R rs (400) and R rs (490) (Figure 6). The correlation coefficient varied between 0.75 for total non-water absorption at 443 nm and 0.89 for absorption by particles at 443 nm. In contrast to MODIS, the OLCI band centred at 443 nm was correlated only with the R rs (412.5) and R rs (490) bands. We also found absorption by particles at 443 nm and Chla concentration to be correlated with R rs (620)/R rs (673.75), with the correlation coefficient close to 0.75. Like MODIS, a CDOM (443) displayed an overall weak correlation with OLCI spectral band ratio combinations. Only the relationship between the ratio R rs (400)/R rs (490) or R rs (708.75)/R rs (767.5) was characterized by a correlation coefficient of >0.3 (Figure 6a). The correlation of a CDOM (443) with the latter band ratio should be regarded as an artifact owing to the null CDOM absorption at wavelengths larger than 700 nm. The results of the regression line fitting were unsatisfactory, so a CDOM (443) was excluded and no empirical formula for estimating this parameter could be established. Identification and verification were conducted in the same manner as before for the MODIS satellite radiometer. On the basis of regression analyses, we found that the most effective band ratio combination to fit all IOPs, except a CDOM (443), was R rs (490)/R rs (560) ( Table 3). The second-order polynomial was identified to be optimal for Chla, a p (443) and a tot (443), as illustrated in Figure 7. Absorption by particles a p (670) and the total non-water absorption at 670 nm a tot (670) were better fitted using linear relationship ( Table 3). The difference was subtle, but higher-order polynomials tended to introduce greater errors at the minimal and maximal data set margins (Figure 7). The thirdand fourth-order polynomials were overfitted. The final versions of proposed algorithm formulas is given in Table 3. Table 3. Proposed empirical second-order polynomial formulas formulations for retrieving Chla, a p (443), a tot (443) and and linear formulas for retrieving a p (670) and a tot (670) from spectral reflectance ratios adjusted to OLCI (Sentinel-3) bands.

Validation of the Proposed Band Ratio Algorithm for Estimating Chlorophyll a Concentration
We conducted a validation study based on Chla values measured in situ during field campaigns in summers of 2013-2015. A total of 24 match-ups were available in 2013 and nearly 20 in summers 2014 and 2015. This number of match-ups produced a dataset that was even bigger than the validation subset. In order to compare the performance our algorithm with standard products, we processed available MODIS ocean colour imagery data using the standard OC3M algorithm [59,70] with two atmospheric correction schemes: the one by Wang and Shi (2007) [62] and the other by Bailey et al. (2010) [61]. We also applied MODIS OC3M and OLCI OC4Me algorithms to in situ remote sensing reflectances to assess their performance in the Nordic Seas without interference from atmospheric disturbances. The chlorophyll a concentration extracted from MODIS AQUA pixels, corresponding to the surface sampling locations that had been registered within a 24-h time lapse, were used for comparison with our algorithm (Figure 8). The difference between the two atmospheric correction methods was negligible, but the atmospheric effects contributed much uncertainty to Chla estimations when compared to in-water algorithms ( Table 4). The Chla estimations using any of the standard MODIS OC2M, or OC3M algorithms applied to in-water measured spectral remote sensing reflectances showed noticeable reduction in uncertainty. The uncertainty of MODIS OC2M was the lowest and OC3M the highest in Chla estimations. There was a small improvement in comparison to the MODIS algorithms when the OLCI OC4Me algorithm was applied to in situ radiometric data. Comparison of Chla concentrations computed from in-water measurements with both of our MODIS and OLCI algorithms, showed them to be clearly superior to standard band ratio algorithms developed for both satellite sensors. In particular, the RMSE and log_RMSE turned out to be lower than for the OC2M and OC3M and OC4Me algorithms. Our algorithms performed better with respect to both bias and RMS (Table 4). Even though they were based on the second-order polynomial, they outperformed the fourth-order polynomial formulations. We did not compare our algorithms for estimating selected IOPs with existing formulas, because we proposed a simple band ratio approach, whereas the eariler algorithms developed earlier were based on semi-analytical models [27,71]. Therefore we relied on the uncertainty budget estimated based on available validation subsets, as presented earlier in the Results section.

Seasonal and Year-To-Year Variability of A P (443) Distribution in the Nordic Seas
We applied our algorithm to the estimation of a p (443), in order to derive the monthly averages in the Nordic Seas, based on MODIS AQUA global coverage remote sensing reflectance data (GSFC NASA Ocean Color Level-3). We used 4 km monthly mosaic data to produce maps of the surface distribution of a p (443) in May, June, July and August in 2013, 2014 and 2015 (see Figure 9). We also calculated averaged values of a p (443) ( Table 5) along three transects (S, Z, and EB2) spanning the Fram Strait between the western Spitsbergen coast and the sea ice edge of Greenland (Figure 1) to quantify seasonal and year-to-year variability and for comparison with in situ data measured in the same area in the same months and years [10,12]. The parameter a p (443) is a very good optical proxy for chlorophyll a concentration and phytoplankton biomass [32,72,73], and its temporal and spatial distribution is largely mediated by phytoplankton growth dynamics. The phytoplankton bloom in the Nordic Seas north of latitude 70 • N start in the eastern part of Fram Strait along the sea ice marginal zone that stretches southwards along the eastern coast of Greenland (Figure 9). In May 2013, this bloom extended northwards along the sea ice marginal zone as far as 78 • N. In May 2014 and 2015 this initial phase of the phytoplankton bloom was contained in an eddy-like structure near the sea ice edge at around latitude 70 • N and longitude 15 • W. As time passed, in June of all the years in question, the phytoplankton bloom was displaced and intensified (increase in a p (443)) towards the northern and central parts of the Greenland Sea and Fram Strait. By July 2013, most of the sea ice tongue drifting out from the central Arctic Ocean with the East Greenland Current (EGC) had melted, and the surface open water contained small amounts of absorbing particulate material. The surface melt water layer was depleted of phytoplankton containing photosynthetic pigments and absorbing detrital particles, and was more transparent than the underlying CDOM-rich EGC waters [7,10,68]. In July 2014 and 2015, there was much more sea ice in the western part of Fram Strait, and the water masses with elevated a p (443)) were located in the central and eastern part of Fram Strait, north of 75 • N, and along the western Spitsbergen Shelf. Elevated a p (443) values were also recorded close to the sea ice marginal zone north of 78 • N near longitude 10 • E. By August of, all three years, the phytoplankton seasonal cycle was in decline, resulting in reduced a p (443) over vast areas of the Nordic Seas. The mean a p (443) values along the three transects in Fram Strait (see Table 5) revealed an interesting pattern: there was a conspicuous south-north gradient in May and June of all three years. In those months, a p (443) was lower in the southern part of the West Spitsbergen Shelf (transect S) but increased steadily towards the north (transects Z and EB2). In the second part of the growing season, mean values of a p (443) were distributed quite evenly along the shelf area, with a slight trend towards higher a p (443) in its southern part (section S). The mean a p (443) values estimated using our algorithm (Table 5) agree very well with ground truth data. We also found the same year-to-year variability pattern as indicated by Kowalczuk et al. (2019) [10] in July of 2013, 2014 and 2015, when most of the in situ data was collected. Table 4. Error comparison of chlorophyll a concentration retrieved from satellite Chla maps corrected using two atmospheric correction methods [61,62]); predictions based on standard OC algorithms applied to in-water measurements and the new algorithm proposed here.

Assessment of Remote Sensing Algorithms for the Retrieval of Inherent Optical Properties, and Chlorophyll a Concentrations in the Arctic Ocean
Ocean colour remote sensing has been used to study the spatial and temporal dynamics of phytoplankton abundance across a range of spatial and temporal scales in the global ocean, regional marine basins and fresh water bodies [23]. However, a number of factors limit the capability of ocean colour remote sensing to provide daily synoptic maps of essential optical and biogeochemical variables at high latitudes. They include low Sun zenith angle, persistent cloudiness and fog, and ice cover that screens the under-ice algal dynamics, as well as further errors in remote sensing reflectance retrievals caused by the so-called ice adjacency effect, the occurrence of a deep chlorophyll maximum and, last but not least, the optical complexity of Arctic Ocean waters, especially over its shelves [31,74]. It has already been found that relationships between remote sensing reflectance band ratios and chlorophyll a concentrations in polar regions differ markedly from those prevailing in lower latitude oceanic waters [75][76][77]. The global remote sensing algorithms were adjusted for Arctic Ocean waters [78] and the new OC4L algorithm, applied to Western Arctic Ocean shelf waters characterized by strong CDOM absorption, retrieved chlorophyll a concentrations with the requisite accuracy for this parameter (RMSE < 35%) [79]. Recently, new regional algorithms have been developed for the Chukchi Sea [80] and Canadian Northeast Pacific and Northwest Atlantic waters for the use of data from the MODIS and VIIRS satellite sensors [81]. Chlorophyll a concentrations retrievals acquired using generic band ratio algorithms were encumbered with a lower uncertainty compared to the standard global ones.
To date there have been few studies for developing a regional algorithm for chlorophyll a concentration retrievals in the polar North Atlantic and European Arctic, except the one by Stramska et al. (2003) [32]. Those authors used a second-order polynomial band ratio algorithm based on the R rs (490)/R rs (555) ratio. Although we followed that approach in our study, we used a significantly larger in situ data set to develop our algorithms. The observed log-normal distribution and formulations of remote sensing algorithms for retrieving IOPs of optical water based on log-transformed data are consistent with previous research [29]. High correlation coefficients with band ratios around 443-490 nm to 550-560 nm were expected. These wavelengths, centred near the chlorophyll a absorption peak in the blue and the minimum in the green, are commonly used for developing remote sensing algorithms in open ocean water [82]. It is interesting that estimations based on the blue channel around 490 nm seem to be more accurate and stable than those based on bands around 443 nm. Our observations are consistent with O'Reilly et al. (2000) [29] who stated that the band ratio based on the bands close to R rs (490)/R rs (555) were the best overall index of chlorophyll a concentration. This is particularly evident in the case of for the OLCI radiometer, which measures radiation in the narrower bands. Intensive absorption leads to a severe decrease in the light level at these wavelengths, which may result in higher signal fluctuations and bigger statistical errors owing to the limited sensitivity of measuring devices. Our samples showed that the a ph (λ) peak was placed around 443 nm and was much higher than at 490 nm. The CDOM contribution to the total non-water absorption a nw (443) in the study area varied between 21% and 51% [10] and decreased rapidly with increasing wavelength, and the R rs (490) waveband was much less affected by CDOM absorption. For this reason, all the algorithms that included the 490 nm waveband performed better than those based on 443 nm. The total non-water absorption studied in this part of the Nordic Seas was dominated by both components of particulate absorption: absorption by phytoplankton pigments, a ph (λ) and absorption by non-pigmented particles a N AP (λ), especially at longer wavelengths. In the near infrared at 670 nm, the absorption budget is controlled almost entirely by phytoplankton pigment absorption alone [10,12]. Because absorption constituents strongly related to chlorophyll a are dominant, we were able to establish reliable bio-optical relationships between the chlorophyll a concentration and absorption by particulates and phytoplankton pigments at 443 nm and 670 nm [10][11][12]32] that can be used in semi-analytical remote sensing reflectance inversion models.
The chlorophyll a concentration estimated from satellite-derived spectral reflectances had a bigger RMSE than the in-water based estimates. This suggests that the atmospheric correction still contributes significantly to the overall error ( Figure 8). The correction using SWIR bands for water masses of high turbidity [62] simultaneously exhibits a smaller bias, but a higher statistical error (Table 4), because MODIS images have a significant noise level in the SWIR part of the spectrum. Water carrying larger amounts of suspended matter scatters in a spectrally-uniform way and raises the radiation levels in the NIR bands, which are commonly used as a zero-level reference to exclude the influence of the atmosphere [61]. In the marginal sea ice zone, the presence of melt water containing large amounts of suspended particles enhances scattering [18,68] and the assumption of negligible values water-leaving radiance in the NIR band may not be fulfilled locally.
The waters of the Nordic Seas are highly productive, especially in the West Spitsbergen Current (WSC), where phytoplankton growth is stimulated by ample supplies of nutrients and the less intensive light attenuation compared to the East Greenland Current EGC [15]. Chlorophyll a concentrations typically peak in the summer from May until July [7,9]. The majority of sampling in the WSC was done in June and July over three years, which could have been responsible for the greater sample homogeneity and the smaller number of extremes. The results obtained with the second-order polynomials were therefore almost as good as with the third-order ones, which may not be true for the whole year. The low level of uncertainty of our algorithms and their consistency in retrieving spectral values of particulate and total non-water absorption coefficients demonstrated that we were right to apply the waveband combination ratio using satellite bands at 490 nm in the empirical formulations. In the cases where there were rapid increases in CDOM absorption, the water-leaving radiance signal recorded by satellite wavebands at shorter wavelengths would be disturbed [22]. New formulations of the empirical relationships defined for the MODIS and OLCI bands displayed certain similarities (Tables 2 and 3). The band ratio indices were based on almost the same wavelengths. Second-order polynomials applied to log-transformed data turned out to be suitable in most cases for both sensors. Even though the satellite sensors have different bandwidths, statistical quality assessment revealed very close levels of algorithms performance. The algorithms for MODIS performed somewhat better, in particular for a p (443) and a tot (443). This is the result of the narrower bandwidths of the OLCI radiometer in relation to MODIS, which seem to be more sensitive to CDOM absorption. With both satellites, moreover, the accuracy of a p (443) and a tot (443) retrieval was less than of a p (670) and a tot (670).
The semi-analytical algorithms for retrieving spectral values of IOPs like GSM [27] or QAA [71] were developed using a dataset in which the majority of data were gathered in temperate and subtropical waters [31]. The accuracy of a semi-analytical algorithm depends closely on parametrizations between phytoplankton pigment absorption and chlorophyll a concentration, and the backscattering spectrum and backscattering phase function. These parameters have changed considerably in the Nordic Seas and European Arctic Ocean. Stramska et al. (2003) [32] reported that phytoplankton communities dominated by diatoms had a higher b b (λ)/b b (555) ratio and that its spectrum was steeper with increasing chlorophyll a concentration compared to phytoplankton communities dominated by dinoflagellates and coccolithophorids. Changes in b b (λ)/b b (555) ratio impacted the blue-green reflectance ratio. Part of the Nordic Seas where phytoplankton communities were dominated by diatoms were characterized by relatively high blue-to-green reflectance ratio, compared to water where phytoplankton communities were dominated by diatoms were characterized by a relatively high blue-to-green reflectance ratio, compared to water in which phytoplankton communities were dominated by dinoflagellates, where the reflectance ratio was much lower. The last two decades have witnessed a significant northward expansion of temperate phytoplankton communities into Arctic Ocean waters [15,16] bringing about significant changes in their bio-optical properties [10,12,83]. The poorer accuracy of semi-analytical algorithms in high latitude oceanic waters, documented by Clay et al. (2019) [81], may be explained by changes in phytoplankton phenology and the associated bio-optical properties in these waters. Figure 3 also clearly shows that reflectance spectra tend to congregate in two clusters with lower and higher values of R rs (λ) in the 400-550 nm spectral range. The cluster of low R rs values is clearly visible in the R rs (λ) spectrum measured in summer 2014 and 2015. Sea surface temperatures in those years were significantly higher than in 2013 owing to the enhanced inflow of warm AW observed in Fram Strait at that time [84]; the intensification of AW inflow contributes to the poleward expansion of temperate phytoplankton in the Arctic Ocean [85]. Stramska et al. (2003) proposed a band ratio algorithm for the retrieval of a tot (442) and a ph (442) [32]. The log-linear parametrization was used to estimate the values of both coefficients, utilizing the R rs (442)/R rs (555) and R rs (490)/R rs (555) band ratios. Our results demonstrated that this band ratio combination was also effective deriving total and particulate absorption at 443 and 670 nm for both the MODIS and OLCI satellite radiometers. The only difference was the second-order polynomial parametrization proposed in this paper for retrieving Chla, a p (443) and a tot (443) from MODIS and OLCI data. The exceptions were the formulas for retrieving a p (670) and a tot (670) from OLCI. For these two cases we suggest linear parametrizations. All the absorption coefficients proposed in this paper are characterized by a low level of uncertainty, and the values of the respective coefficients are consistent with the data reported by Kowalczuk et al. 2017 and 2019 [10,12].
CDOM in AW dominating in the eastern part of Fram Strait has a relatively low level of absorption and is predominantly of marine origin [9,86]. The CDOM absorption coefficient at 350 nm a CDOM (350) is on average up to three times lower than of the water masses influenced by EGC [7]. Statistically significant inter-annual variabilities in the level of CDOM absorption in Atlantic-dominated waters between 2009-2010 were reported by Pavlov et al. (2015) [7] and between 2013-2015 by Makarewicz et al. (2018) [9]. In the Nordic Seas, the IOPs of water on the opposite sides of the Polar front differ significantly as regards the proportions of absorbing constituents [7,13,87]. The CDOM in the western part of Fram Strait and on the East Greenland shelf CDOM originates from the riverine input of Siberian Rivers to the Arctic Ocean [7,13,14]. The CDOM absorption level in PW carried southwards with EGC, though much higher than in AW, is temporally more stable [13,88,89] and could be effectively lowered by fresh, glacier melt waters from the Greenland ice sheet [88,89] or drift sea ice melt [68]. In spite of the broad spatial and temporal variability in CDOM absorption in the Nordic Seas, the accuracy of Chla, a p (443), a p (670), a tot (443) and a tot (670) retrievals in our algorithms was not sensitive to changes in this parameter.

The Role of Sea Ice Cover in the Spatial and Temporal Variability of Inherent Optical Properties in the Nordic Seas
Particulate absorption at 443 and 670 nm is a very reliable proxy of the chlorophyll a concentration in the Nordic Seas [10,12,32] and is a convenient tool for observing phytoplankton dynamics. The monthly averaged a p (443) products from May to August in three consecutive years from 2013 to 2015 ( Figure 9) revealed a distinct temporal and spatial zonation of phytoplankton dynamics in the Nordic Seas. The bloom started in May along the sea ice marginal zone, along the Greenland shelf in the western part of Fram Strait, whereas in the eastern part of Fram Strait and along the West Spitsbergen shelf, phytoplankton abundance was relatively low. As summer progressed, the phytoplankton was displaced northwards along the sea ice marginal zone and eastwards towards the central part of Fram Strait, but the bloom intensity decreased (Figure 9). In late summer (July and August), the phytoplankton concentration declined. It was distributed in the central and eastern parts of Fram Strait, and followed the northern part of the sea ice edge in the Nansen Basin and on the northern Spitsbergen Shelf. This zonation is consistent with the previous observations using the satellite imagery of chlorophyll a concentrations [90] and the modelling studies [91]. Mayot et al. (2020) [92] also found that the intensity of spring sea ice exported from the Arctic Ocean through Fram Strait had a significant influence on Greenland Sea phytoplankton abundance. He reported that in those years when sea ice export was high, stratification was stronger and phytoplankton blooms were more intense in the Greenland Sea. We, too, witnessed such an effect in the monthly averaged a p (443) products when comparing July 2014 and 2015: in July 2014 the phytoplankton bloom in the central Fram Strait was more intense than in July 2015. This is wholly consistent with the in situ observations presented by Kowalczuk et al. (2019) [10]. In 2014, the mostly southerly winds brought about advection of sea ice towards the central Fram Strait, intensifying stratification and producing a stronger phytoplankton bloom. In July 2015 the predominantly northerly winds caused sea ice to converge along the east coast of Greenland: the oceanic waters in the central and eastern parts of Fram Strait were less stratified and the phytoplankton bloom was weaker.
A narrow belt of less absorbing water is usually visible between the ice edge and the areas with elevated a p (443) (Figure 9). This is the effect of multiple processes due to the broad temporal averaging bracket (one month) used in mosaicking. Such a long time leads to strong fluctuations in the ever-moving ice edge position. Another process is the sea adjacency effect, which influences atmospheric correction performance. Last but not least is the vertical distribution of phytoplankton. The presence of melt water in the sea ice marginal zone generates pycnocline, which forces the phytoplankton to aggregate along the density gradient between the melt and oceanic waters, creating a deep chlorophyll a maximum [10,93]: its presence means that part of the phytoplankton biomass is located at depths well below satellite radiometer detection limits [90]. As summer progressed, the sea ice melted and vast areas of the Eastern Greenland shelf waters were exposed, but we did not record any elevated absorption. This must have been due to the optical properties of the melt water, which is low in phytoplankton pigments and containing diluted CDOM compared to the underlying CDOM-rich East Greenland Current waters [13,68]. The mineral particles that melt out from ice and are contained in the melt water made a greater contribution to scattering than to absorption [18].

Conclusions
Our research confirmed that even in the wake of recent changes in Nordic Sea phytoplankton communities, the uncertainty levels of the empirical remote sensing algorithms for retrieving chlorophyll a concentrations and IOPs based on spectral blue-green band-ratio combinations were low. The bulk optical properties in the sampled water masses were governed by phytoplankton dynamics, so the most effective models for retrieving Chla, a p (443), a p (670), a tot (443) and a tot (670) were based on the ratio index using log-transformed bands around the 480-490 nm and 550-560 nm wavelengths. The exponent in the form of a second-order polynomial reconstructed these optical variables in most cases with the very high degree of accuracy: e.g., the chlorophyll a concentration based on OLCI bands was estimated with bias 3% and RMSE 34%. Total non water absorption as well as particulate absorption at 670 nm were best estimated using just the first-order polynomial for the exponent (the bias of a tot (670) was −0.62% and RMS 28.32%, the bias of a p (670) was 0.28% and 28.45% for the OLCI bands). The higher-order polynomials not only failed to improve the results, but tended to cause model over-fitting and a deterioration in accuracy. Our algorithms have been used for the surface distribution of Chla, a p (443), a p (670), a tot (443) and a tot (670) demonstrating their potential operational use, especially for new satellite sensors like the OLCI.

Abbreviations
The following abbreviations are used in this manuscript: