An Empirical Ocean Colour Algorithm for Estimating the Contribution of Coloured Dissolved Organic Matter in North-Central Western Adriatic Sea

The performance of empirical band ratio models were evaluated for the estimation of Coloured Dissolved Organic Matter (CDOM) using MODIS ocean colour sensor images and data collected on the North-Central Western Adriatic Sea (Mediterranean Sea). Relationships between in situ measurements (2013–2016) of CDOM absorption coefficients at 355 nm (aCDOM355) with several MODIS satellite band ratios were evaluated on a test data set. The prediction capability of the different linear models was assessed on a validation data set. Based on some statistical diagnostic parameters (R2, APD and RMSE), the best MODIS band ratio performance in retrieving CDOM was obtained by a simple linear model of the transformed dependent variable using the remote sensing reflectance band ratio Rrs(667)/Rrs(488) as the only independent variable. The best-retrieved CDOM algorithm provides very good results for the complex coastal area along the North-Central Western Adriatic Sea where the Po River outflow is the main driving force in CDOM and nutrient circulation, which in winter mostly remains confined to a coastal boundary layer, whereas in summer it spreads to the open sea as well.


Introduction
Ocean colour imagery collected over the past decades from many satellites has the potential for providing synoptic retrievals of bio-optical properties in the world's oceans [1,2].Unfortunately, the optical complexity of coastal waters affects remote sensing products and their accuracy.Thus, the pressing need of validating large-scale remote sensing coastal products makes attractive the availability of cross-site consistent in situ observations performed applying comprehensively assessed and standardized measurement protocols.Improving retrievals of water borne constituents in coastal systems requires the optimization of algorithms using regional in situ measurements of optical and biogeochemical water properties.
In this context, the Coloured Dissolved Organic Matter (CDOM) represents the optically active fraction of DOM in natural waters and plays various roles in physical and biogeochemical processes.CDOM absorption is characterized by an exponential absorption decrease from ultraviolet (UV) to visible wavelengths.It can dominate the inherent light absorption at the blue wavelengths in surface waters of the coastal seas (20%-70% at 440 nm; [3,4]), pelagic seas and oceans (>50% at 440 nm; [5]), hence competing with phytoplankton for photosynthetically active radiation [6,7].In many coastal areas, CDOM absorption is several times that of chlorophyll and confounds the retrieval of the latter from ocean colour satellite observations due to overlapping absorbance spectra at the blue wavelengths [8][9][10].Furthermore, CDOM has proven to be a useful tracer not only for carbon but also as a proxy for mixing in a wide variety of environments [11][12][13].
Primary sources of CDOM are rivers and groundwater near coastlines, which carry CDOM primarily from soils, but coastal waters, can contain plankton-derived CDOM produced in rivers and estuaries, as well as anthropogenic compounds from runoff, sewage discharge and other effluents [14,15].However, biological processes such as phytoplankton growth, zooplankton grazing and microbial activity can also contribute to marine-derived CDOM in continental margins and pelagic ocean [16][17][18][19].On the other hand, photobleaching is the dominant process for CDOM removal from natural waters [20], while microbial decomposition is of less importance [21,22].The balance between sources and sinks controls the patterns of CDOM distribution.
Historically, many ocean colour algorithms have been developed to retrieve biogeochemical properties (e.g., chlorophyll-a, total suspended matter and CDOM) based on direct empirical relationship with the remote sensing reflectance, R rs (λ), or ratios of R rs at various wavelengths.Many of these models estimate the absorption coefficient of CDOM (a CDOM ) and detrital (non-pigmented) particles as a single parameter (a CDM ), because CDOM and detritus have similar spectral responses in the visible spectrum [23][24][25][26][27][28][29][30][31][32][33].The retrieval of a CDM by these algorithms yields reasonable results in open ocean regions.However, such algorithms typically do not work well in coastal waters due to the optical complexity of turbid coastal waters (high levels of CDOM, coloured detrital particles, and phytoplankton; [34]).
During the last years, many authors [32,35,36] have attempted to separate a CDOM from a CDM and to derive the CDOM spectral slope in various spectral regions [37][38][39].Other authors [34,40] described how retrievals of bio-optical properties with a spectral references band near 640 nm showed significant improvement over the standard 555 nm band because total light absorption near 640 nm is primarily greatest for seawater itself rather than for other constituents in the water column.Anyway, the primary limitation to rigorous validation is the lack of sufficient data of coincident field measurements and satellite observations that are independent from the data used to develop the algorithms.
Satellite ocean colour optical sensors currently operating are: (a) SEVIRI (Spinning Enhanced Visible and Infrared Imager) onboard the Meteosat Second Generation (MSG-3; 2005-present) geostationary platforms with a Ground Sampling Distance (GSD) of 1-4 km; (b) Moderate Resolution Imaging Spectroradiometer on board the polar-orbiting Terra/Aqua satellites (MODIS-A; 2002-present) and the Visible Infrared Imager/Radiometer Suite (VIIRS; 2011-present) on the NPOESS Preparatory Project (NPP) satellite mission, with a GSD of 250-1000 m; (c) the Operational Land Imager (OLI; 2013-present) on the polar-orbiting Landsat-8 satellite with a GSD of 30 m; and (d) Sentinel-2A (launched on 23 June 2015) and Sentinel 3A (launched on 16 February 2016 and available from 17 November 2016) polar-orbiting multispectral high-resolution imaging ESA missions (Copernicus Services) with a GSD of 10-300 m.The data recorded by these satellite optical sensors cover the range of temporal (15 min to 15 days) and spatial (10 m to 4000 m) resolutions that are currently free available for the monitoring of sea surface dynamics.
In this study, we evaluate and analyse the performance of the MODIS-A-derived CDOM estimations using the Copernicus Marine Environment Monitoring Service CMEMS Remote sensing reflectance (R rs at wavelengths of 412, 443, 488, 531, 547 and 667 nm) for the North-Central Western Adriatic Sea (NCWAS) in Italy, which is a complex system influenced by the Po River, one of the largest Mediterranean rivers [41].Ocean colour satellite retrievals of a CDOM could have many potential applications and in this work, some regional CDOM empirical algorithms for the NCWAS were tested and validated on MODIS-Aqua-measured R rs (λ) and a CDOM in situ data were collected from 2013 to 2016.After validation, the optimal model was used to produce a CDOM maps of the NCWAS to describe his distribution throughout the basin.The results of the work are drawn from field measurements supported by multiple research projects in order to collect sufficient coincident satellite data with field observations.

Study Site
The North Adriatic Sea (Figure 1) is the northernmost part of the Mediterranean Sea, extending as far North as 45 • 47 N.The area is strongly influenced by river floods [42][43][44][45], which affect both the circulation through buoyancy input and the ecosystem by introducing a large amount of matters [46][47][48].
The Po River provides the major buoyancy flux with an annual mean freshwater discharge rate of about 1500 m 3 •s −1 [41,49].On the annual basis, the Po River alone carries 28% of the total Adriatic Sea runoff (5500−5700 m 3 •s −1 ), 45% comes from the eastern coast (about half of which from Albania), 19% from the northern coast and 8% from the entire western coast [41,50].The southward coastal flow, i.e., the Western Adriatic Current (WAC [50][51][52][53]), is driven by the Po River buoyancy flux (low-salinity waters) and northeastern Bora winds that characterize this region during the winter months.Bora winds cause elevated sea surface height along the western coasts, producing the downwelling and transport of coastal dense waters toward the open sea [54].

Study Site
The North Adriatic Sea (Figure 1) is the northernmost part of the Mediterranean Sea, extending as far North as 45°47′N.The area is strongly influenced by river floods [42][43][44][45], which affect both the circulation through buoyancy input and the ecosystem by introducing a large amount of matters [46][47][48].The Po River provides the major buoyancy flux with an annual mean freshwater discharge rate of about 1500 m 3 •s −1 [41,49].On the annual basis, the Po River alone carries 28% of the total Adriatic Sea runoff (5500−5700 m 3 •s −1 ), 45% comes from the eastern coast (about half of which from Albania), 19% from the northern coast and 8% from the entire western coast [41,50].The southward coastal flow, i.e., the Western Adriatic Current (WAC [50][51][52][53]), is driven by the Po River buoyancy flux (lowsalinity waters) and northeastern Bora winds that characterize this region during the winter months.Bora winds cause elevated sea surface height along the western coasts, producing the downwelling and transport of coastal dense waters toward the open sea [54].The Po River freshwater may flow southward along the shelf strengthening the Western Adriatic Current (WAC) or cross the shelf extending toward Istria in the central basin [55].In particular, during periods of weak stratification, in absence of wind forcing, vortices constrain fresh Po River waters to a southward flow along the Italian shelf [56].In contrast, in periods of stratification, particularly in spring and summer, the Po River fresh water plume spreads across the basin to the Istrian coast to form a front that divides the northern basin [57,58].The extent of Po riverine flow and wind forcing, therefore, modulate fresh water penetration in the North Adriatic [59].The Po River freshwater may flow southward along the shelf strengthening the Western Adriatic Current (WAC) or cross the shelf extending toward Istria in the central basin [55].In particular, during periods of weak stratification, in absence of wind forcing, vortices constrain fresh Po River waters to a southward flow along the Italian shelf [56].In contrast, in periods of stratification, particularly in spring and summer, the Po River fresh water plume spreads across the basin to the Istrian coast to form a front that divides the northern basin [57,58].The extent of Po riverine flow and wind forcing, therefore, modulate fresh water penetration in the North Adriatic [59].
Eutrophication, which periodically occurs along the Italian coast of the Adriatic Sea [60][61][62][63], is mainly caused by the discharge of organic and inorganic matter from the Po River and detrital coming from other rivers flowing into the NCWAS and contributes to make this area a complex optical system due to simultaneous presence of CDOM, coloured detrital particles, and phytoplankton [34].Therefore, in this work, only CDOM-based algorithms (a CDM -based algorithms works better for other study sites where detrital particles concentrations are not so high) were considered, as CDOM is the variable of main interest for the study area; moreover, detrital particles and phytoplankton were not measured as extensively as CDOM.

Field Measurements of Seawater
As part of numerous national and international projects, several oceanographic cruises in NCWAS were conducted between November 2013 and July 2016.One hundred seventy-two water samples, corresponding to clear-sky days during which MODIS-A satellite has acquired images over the selected study area, were collected (Figure 1; Table 1).For each of these water samples, CDOM and salinity were measured at the sea surface level (in situ CDOM).Salinity (S) of the surface waters was measured using a SeaBird Electronics 911-plus CTD (Conductivity Temperature Depth).The CTD data were processed according to UNESCO standards [64].
Seawater samples for analysis of CDOM absorbance spectra were filtered through sterile Whatman GD/X 0.2 µm filters, under low pressure.The filters were washed with 100 mL of sample before its collection in order to avoid any CDOM contamination.Samples were: (a) filtered immediately after the sampling in order to avoid alteration due to microbial activity; (b) collected in amber glass bottles, previously acid-soaked (10% HCl); (c) rinsed with Milli-Q water; and (d) rinsed three times with the sample before its collection.The filtered water was kept at 4 • C until the optical measurements were performed, at the latest within four weeks.
CDOM absorbance was measured throughout a dual beam UV-VIS spectrophotometer (SHIMADZU 2600 Series), using a 10 cm path length cells with ultraviolet oxidized Milli-Q water as the blank and reference [65].Instrument scan settings were as follows: 250-750 nm wavelength scan range, fast scan speed; 1 nm sampling interval; 0.5 nm slit width.The quartz cuvettes for blanks and samples were acid-soaked for one hour and then rinsed with Milli-Q water and sample aliquots.The absorbance (A) was converted into absorption coefficients (a CDOM λ) by Equation (1): where A λ is the absorbance at wavelength λ and L is the path length expressed in meters.Before determining a CDOM λ, absorbance data were corrected by subtracting the mean absorbance from 680 to 690 nm from each wavelength.The absorbance measurements ranged from 250 to 450 nm representing wavelengths where the signal was sufficiently high for a reliable estimation of a CDOM λ [48,66].Specifically, the absorption coefficients at 355 nm (a CDOM 355) were calculated because this wavelength shows the maximum excitation for humic-like substances [66], i.e., those of terrestrial origin.Since the NCWAS is an area impacted by high river discharges we chose this absorption coefficient for the retrieval algorithms on the assumption that most of CDOM is of terrestrial origin.

Satellite Data and Processing
MODIS-Aqua data over the study areas in the NCWAS were acquired from CMEMS online catalogue [67].MODIS is an ocean colour sensor aboard the polar-orbiting Aqua (MODIS-A) and Terra (MODIS-T) satellite platforms with 36 spectral bands and three spatial resolutions of 250 m, 500 m and 1 km.For the Mediterranean Sea, R rs and diffuse attenuation coefficient of light at 490 nm is operationally produced and uploaded on the CMEMS catalogue by the Group for Satellite Oceanography (GOS-ISAC) of the Italian National Research Council, in Rome, for near real time data from MODIS-A and NPP-VIIRS sensors [1,68,69].
R rs (λ) is the fundamental quantity to be derived from ocean colour sensors and it is defined as the ratio of upwelling radiance and downwelling irradiance at any wavelength that can be expressed as the ratio of normalized water leaving radiance and the extra-terrestrial solar irradiance.The spectral remote-sensing reflectance R rs is defined as [70]: where the depth argument of "in air" indicates that R rs is calculated using the water-leaving radiance L w and E d in the air, just above the water surface.The R rs is a measure of how much of the downwelling radiance that is incident onto the water surface in any direction is eventually returned through the surface into a small solid angle ∆Ω centered on a particular direction (θ, φ).However, R rs is usually computed for nadir-viewing directions only.
In this study, nominal MODIS-A six R rs bands (412, 443, 488, 531, 547 and 667 nm) were used for CDOM retrieval, with a spatial resolution of 1.1 km and acquired in the same dates of the field campaign period 2013-2016 (temporal window between satellite overpass and the time of field sampling = ±12 h). Figure 2a shows an example of R rs MODIS imagery (R: 547 nm; G: 488 nm; B: 443 nm) resized on the Adriatic Sea in Italy centred on the Ancona coastal area (red box in Figure 2a).MODIS R rs imagery was acquired on 22 July 2015, i.e., on the same date of the eleven measurement stations performed in these coastal areas (white dots in Figure 2b; Ecosee/a 4 Cruise).In Figure 2c the graph shows the R rs spectra relative to the eleven coastal measurement stations.As expected, R rs spectra relative to stations 3, 10 and 11 show higher values in the blue band at 443 nm.This is also evident from the MODIS R rs RGB imagery in Figure 2 that shows both the CDOM gradient along the coastal area occurring in this date of acquisition (22 July 2015) and the higher values in the blue band at 443 nm for the stations more distant from the coast (i.e., 3, 10 and 11).

CDOM Models
Due to the lack of coincident in situ measurements and satellite observations, previous attempts to retrieve bio-information from satellite data in this complex optical region found better results on the retrieval of chlorophyll and suspended matter rather than of CDOM [71][72][73][74].Having a consistent number of in situ and satellite data, the main aim of this work is on the definition of empirical models useful to effectively retrieve aCDOM355 values using the optical satellite data acquired on the northwestern Adriatic coastal sea waters.Simple linear algorithms were chosen among several algorithms available in literature for CDOM retrieval, as they proved to be nearly effective with respect to the more complex models in CDOM retrieval [39].On the other hand, empirical algorithms are known to be valid only regionally as they are particularly sensitive to changes in the specific composition of water constituents when boundary conditions are changed [75].For the fitting and evaluation of our empirical models, the collected data were divided into two data sets: data collected from November 2013 to May 2015 (test data set, 93 observations) and data collected from July 2015 to July 2016 (validation data set, 79 observations).The initial data set was split according to two different time intervals in order to assess the prediction capability of the models.The spatial coverage of the two data sets was about the same, Figure 1.Even if the spatial coverage could be considered not very extensive, it must be kept in mind that this limitation was due to the necessity to have corresponding data from in situ and satellite measurements.The mathematical models were fitted on the test data set and the prediction performance evaluated on the validation data set.In all cases, a least squares analysis was carried out in order to identify the best relationship between the dependent variable, the CDOM absorption coefficient aCDOM355, and the independent variables.As independent variables, the following MODIS Rrs band ratios were considered: (a) Rrs(547)/Rrs( 412 The inverse band ratios (a) and (b), i.e., Rrs(412)/Rrs(547) and Rrs(488)/Rrs(547), were used to define empirical algorithms in coastal waters by many authors [23,37,39,76].However, the reference wavelength at 547 nm is affected by particles light scattering [77] that in coastal waters interferes with the CDOM retrieval.Since the NCWAS is strongly influenced by freshwaters/terrestrial input, a wavelength in the red or near infrared (667 nm) was added as numerator of the band ratios (c) and (d) in order to reduce the effects of detritus particles, as suggested by other authors [40,78,79].At first, a simple multilinear model was tested by ordinary least square (OLS) using both the dependent and

CDOM Models
Due to the lack of coincident in situ measurements and satellite observations, previous attempts to retrieve bio-information from satellite data in this complex optical region found better results on the retrieval of chlorophyll and suspended matter rather than of CDOM [71][72][73][74].Having a consistent number of in situ and satellite data, the main aim of this work is on the definition of empirical models useful to effectively retrieve a CDOM 355 values using the optical satellite data acquired on the northwestern Adriatic coastal sea waters.Simple linear algorithms were chosen among several algorithms available in literature for CDOM retrieval, as they proved to be nearly effective with respect to the more complex models in CDOM retrieval [39].On the other hand, empirical algorithms are known to be valid only regionally as they are particularly sensitive to changes in the specific composition of water constituents when boundary conditions are changed [75].For the fitting and evaluation of our empirical models, the collected data were divided into two data sets: data collected from November 2013 to May 2015 (test data set, 93 observations) and data collected from July 2015 to July 2016 (validation data set, 79 observations).The initial data set was split according to two different time intervals in order to assess the prediction capability of the models.The spatial coverage of the two data sets was about the same, Figure 1.Even if the spatial coverage could be considered not very extensive, it must be kept in mind that this limitation was due to the necessity to have corresponding data from in situ and satellite measurements.The mathematical models were fitted on the test data set and the prediction performance evaluated on the validation data set.In all cases, a least squares analysis was carried out in order to identify the best relationship between the dependent variable, the CDOM absorption coefficient a CDOM 355, and the independent variables.As independent variables, the following MODIS R rs band ratios were considered: (a) R rs (547)/R rs (412); (b) R rs (547)/R rs (488); (c) R rs (667)/R rs (412); (d) R rs (667)/R rs (488); (e) R rs (667)/R rs (443); and (f) R rs (531)/R rs (412).The inverse band ratios (a) and (b), i.e., R rs (412)/R rs (547) and R rs (488)/R rs (547), were used to define empirical algorithms in coastal waters by many authors [23,37,39,76].However, the reference wavelength at 547 nm is affected by particles light scattering [77] that in coastal waters interferes with the CDOM retrieval.Since the NCWAS is strongly influenced by freshwaters/terrestrial input, a wavelength in the red or near infrared (667 nm) was added as numerator of the band ratios (c) and (d) in order to reduce the effects of detritus particles, as suggested by other authors [40,78,79].At first, a simple multilinear model was tested by ordinary least square (OLS) using both the dependent and the independent variables transformed according to the log function [39].This first analysis provided a quite low value of the adjusted coefficient of determination (adjusted R 2 = 0.5671).Therefore, we started using a simple multilinear model with untransformed variables and, as shown in the Results Section, proceeded according to model check and diagnostics.Once obtained some final models, their prediction capability was assessed against the validation data set.
The accuracy of model prediction was evaluated according to the following statistical indicators: (a) the adjusted coefficient of determination (adjusted R 2 ) between the in situ validation values a CDOM 355 is 's and the calculated values a CDOM 355 c 's; (b) the mean and standard deviation of the absolute per cent difference (APD) (APD = 100 × |a CDOM 355 is − a CDOM 355 c |/a CDOM 355 is ); and (c) the root mean square error (RMSE) between a CDOM 355 is 's and a CDOM 355 c 's.All the statistical analyses were performed using the open source statistical software R ver.3.3.1 [80] at a significance level of 0.05.
Finally, the retrieved optimal model was used to construct CDOM time maps for this geographical area to further evaluate the applicability of the a CDOM 355 algorithm.These maps were then compared with the corresponding Po river discharges, which as aforementioned has a great influence in CDOM and nutrient circulation in the northwestern Adriatic Sea coastal areas.

Field Distribution of CDOM
The surface distribution of salinity and CDOM (a CDOM 355) are summarized in Table 2.Because of the complexity of the basin, the NCWSA was divided in three sub-areas: the northern part (Zone A), directly influenced from Po River runoff, and the central (Zone B) and southern parts (Zone C) both influenced by minor river discharges and indirectly by Po River runoff carried out along the northwestern Adriatic coastal areas by the WAC.River discharges along the Western Adriatic Sea coastal areas enrich seawater in CDOM, especially in Zone A (Figure 1).A strong linear relationship between a CDOM 355 and salinity confirms this assumption (R 2 = 0.80; Figure 3).Salinity is a useful parameter commonly used to characterize the main terrestrial origin of CDOM, i.e., riverine CDOM, along the NCWAS.
The NCWAS is characterized by a decrease of CDOM from north to south and from coast to off-shore coupled with a salinity increase as expected.Throughout the study period a wide range of CDOM and salinity (S) are observed (S ranges from 18 to 38.7 and a CDOM 355 from 0.08 to 2.18 m −1 ), as shown by the high standard deviations (Table 2).The high spatial and temporal variability of the data is however useful for satellite observation in order to define a properly CDOM retrieval algorithm.River discharges along the Western Adriatic Sea coastal areas enrich seawater in CDOM, especially in Zone A (Figure 1).A strong linear relationship between aCDOM355 and salinity confirms this assumption (R 2 = 0.80; Figure 3).Salinity is a useful parameter commonly used to characterize the main terrestrial origin of CDOM, i.e., riverine CDOM, along the NCWAS.

Simple OLS Linear Regression
Before starting the model tuning, a simple correlation analysis was carried out among the independent variables, i.e., the band ratios, to avoid the problem of collinearity.It turned out that the band ratio R rs (667)/R rs (443) was nearly perfectly correlated with the band ratio R rs (667)/R rs (488), Pearson correlation coefficient = 0.991, and the band ratio R rs (531)/R rs (412) was nearly perfectly correlated with the band ratio R rs (547)/R rs (412), Pearson correlation coefficient = 0.996.Since the aim of this work is the development of a good prediction algorithm for the dependent variable [81], the band ratios R rs (667)/R rs (443) and R rs (531)/R rs (412) were dropped in the following analysis.According to [78], the band ratio R rs (667)/R rs (488) was maintained because riverine CDOM still shows significant absorption at 488 nm, whereas chlorophyll does not.The band ratio R rs (547)/R rs (412) was maintained in order to compare our results to previous works [23,37].
At the beginning, a simple multilinear model with untransformed variables was chosen as starting algorithm: where the b i 's are the parameters to be estimated using satellite R rs band ratios and the in situ a CDOM 355 values (test data set, 93 observations collected from November 2013 to May 2015).The results of the OLS regression for this model on the test data set are briefly summarized in Table 3.The adjusted coefficient of determination R 2 of this model was 0.8502.From Table 3, only the intercept and the band ratio R rs (667)/R rs (488) can be considered significant for the model.Therefore, the following simplified algorithm (Model 1) was evaluated: The regression analysis of Model 1 provided the results presented in Table 4.The Model 1 fit is shown in Figure 4 along with the experimental data.The adjusted coefficient of determination R 2 was 0.8433.This high value of adjusted R 2 is comforting to pursue further the simple linear model.The diagnostic check of Model 1 was carried out according to the normal distribution of the residuals and the constancy of the variance (i.e., the homoscedasticy).The Shapiro-Wilk normality test on the model residuals gave a P value of 0.7093.However, as shown in Figure 5, the plot of the residuals against the fitted values displays a variance which increases with the increment of the fitted value (i.e., the heteroscedasticity).Heteroscedasticity was also confirmed by the studentized Breusch-Pagan test [82], which gave a p-value of 1.091 × 10 −6 .Besides heteroscedasticity, the residuals trend is slightly upward concave, Figure 5.This additional feature is discussed in Section 3.2.4.Even if the OLS estimates of the parameters and the coefficient of determination (R 2 ) are not biased by heteroscedasticity, this violation affects the confidence intervals of OLS parameters [82].Two different approaches can be adopted to overcome heteroscedasticity.The first approach is a dependent variable transformation according to the Cox-Box analysis [83].The second approach is to carry out a generalized least square (GLS) regression with a specific variance structure [84].Both approaches are adopted in the following.

Dependent Variable Transformation
The dependent variable transformation was carried out according to the Box-Cox transformation, i.e., y′ = y α , where y is the original dependent variable, y' the transformed variable and α the transformation exponent.The plot in Figure 6, as produced by the R package MASS [85], suggests a square root transformation [83], i.e., y′ = y 0.5 , with the following Model 2 to be tested: The diagnostic check of Model 1 was carried out according to the normal distribution of the residuals and the constancy of the variance (i.e., the homoscedasticy).The Shapiro-Wilk normality test on the model residuals gave a P value of 0.7093.However, as shown in Figure 5, the plot of the residuals against the fitted values displays a variance which increases with the increment of the fitted value (i.e., the heteroscedasticity).Heteroscedasticity was also confirmed by the studentized Breusch-Pagan test [82], which gave a p-value of 1.091 × 10 −6 .Besides heteroscedasticity, the residuals trend is slightly upward concave, Figure 5.This additional feature is discussed in Section 3.2.4.
Remote Sens. 2017, 9, 180 9 of 22 The diagnostic check of Model 1 was carried out according to the normal distribution of the residuals and the constancy of the variance (i.e., the homoscedasticy).The Shapiro-Wilk normality test on the model residuals gave a P value of 0.7093.However, as shown in Figure 5, the plot of the residuals against the fitted values displays a variance which increases with the increment of the fitted value (i.e., the heteroscedasticity).Heteroscedasticity was also confirmed by the studentized Breusch-Pagan test [82], which gave a p-value of 1.091 × 10 −6 .Besides heteroscedasticity, the residuals trend is slightly upward concave, Figure 5.This additional feature is discussed in Section 3.2.4.Even if the OLS estimates of the parameters and the coefficient of determination (R 2 ) are not biased by heteroscedasticity, this violation affects the confidence intervals of OLS parameters [82].Two different approaches can be adopted to overcome heteroscedasticity.The first approach is a dependent variable transformation according to the Cox-Box analysis [83].The second approach is to carry out a generalized least square (GLS) regression with a specific variance structure [84].Both approaches are adopted in the following.

Dependent Variable Transformation
The dependent variable transformation was carried out according to the Box-Cox transformation, i.e., y′ = y α , where y is the original dependent variable, y' the transformed variable and α the transformation exponent.The plot in Figure 6, as produced by the R package MASS [85], suggests a square root transformation [83], i.e., y′ = y 0.5 , with the following Model 2 to be tested: Even if the OLS estimates of the parameters and the coefficient of determination (R 2 ) are not biased by heteroscedasticity, this violation affects the confidence intervals of OLS parameters [82].Two different approaches can be adopted to overcome heteroscedasticity.The first approach is a dependent variable transformation according to the Cox-Box analysis [83].The second approach is to carry out a generalized least square (GLS) regression with a specific variance structure [84].Both approaches are adopted in the following.

Dependent Variable Transformation
The dependent variable transformation was carried out according to the Box-Cox transformation, i.e., y = y α , where y is the original dependent variable, y' the transformed variable and α the transformation exponent.The plot in Figure 6, as produced by the R package MASS [85], suggests a square root transformation [83], i.e., y = y 0.5 , with the following Model 2 to be tested:

Generalized Least Squares (GLS)
Another way of tackling the problem of heteroscedasticity is to perform a generalized least squares (GLS) regression analysis.In this case, the dependent and independent variables are untransformed, but the variance of the dependent variable is let to vary instead of being considered constant as in the OLS.Model 3 is, therefore, given by Model 3: aCDOM355 = b1 + b2 × Rrs(667)/Rrs(488) (6) Table 5 reports the regression results obtained for Model 2 on the test data set.The model fit is shown in Figure 7 along with the experimental data.The adjusted coefficient of determination was 0.8369.

Generalized Least Squares (GLS)
Another way of tackling the problem of heteroscedasticity is to perform a generalized least squares (GLS) regression analysis.In this case, the dependent and independent variables are untransformed, but the variance of the dependent variable is let to vary instead of being considered Both the plots of the residuals (not shown) and the studentized Breusch-Pagan test (p = 0.7571) indicate that the problem of heteroscedasticity was solved.

Generalized Least Squares (GLS)
Another way of tackling the problem of heteroscedasticity is to perform a generalized least squares (GLS) regression analysis.In this case, the dependent and independent variables are untransformed, but the variance of the dependent variable is let to vary instead of being considered constant as in the OLS.Model 3 is, therefore, given by Model 3: a CDOM 355 = b 1 + b 2 × R rs (667)/R rs (488) (6) with the fitting parameters b i 's to be estimated by GLS.The fitting results of this model, which were carried out by using the R package "nlme" with a variance structure of the type "varPower()" [86], are shown in Table 6. Figure 8 shows the experimental data and the model fit.The adjusted coefficient of determination was 0.6325.Comparing the values of intercept b 1 and slope b 2 of Model 1 (Table 4) and Model 2 (Table 6), we can see that they are statistically different.This is unexpected, since heteroscedasticity should not affect the estimate of the regression parameters.Evidently, heteroscedasticity is not the only issue that emerges from the analysis of the residuals plot.In the next section, this problem is analysed more deeply.with the fitting parameters bi's to be estimated by GLS.The fitting results of this model, which were carried out by using the R package "nlme" with a variance structure of the type "varPower()" [86], are shown in Table 6. Figure 8 shows the experimental data and the model fit.The adjusted coefficient of determination was 0.6325.Comparing the values of intercept b1 and slope b2 of Model 1 (Table 4) and Model 2 (Table 6), we can see that they are statistically different.This is unexpected, since heteroscedasticity should not affect the estimate of the regression parameters.Evidently, heteroscedasticity is not the only issue that emerges from the analysis of the residuals plot.In the next section, this problem is analysed more deeply.

Curvature of the Residuals
As shown in Figure 5, a slight curved trend is present in the residuals of Model 1.This curvature usually implies the presence of higher order terms in the model [87].In this section, besides the possibility of introducing terms of higher order, also the interaction between independent variables was considered.The detailed steps of this kind of analysis are described in Crawley [88].According to this procedure, preliminarily to the regression analysis, two explorative investigations should be carried out: (1) a regression tree to see if the model has a complex structure, i.e., there are interaction terms between the independent variables; and (2) a generalized additive model (GAM) analysis to see if high order terms are necessary.Once ascertained the general model, a repeated OLS regression is performed cancelling, each time, the non-significant terms starting from the highest interaction terms.When all the surviving terms are significant, the regression procedure stops.However, if heteroscedasticity is present in the final model diagnostic, the regression procedure must be repeated

Curvature of the Residuals
As shown in Figure 5, a slight curved trend is present in the residuals of Model 1.This curvature usually implies the presence of higher order terms in the model [87].In this section, besides the possibility of introducing terms of higher order, also the interaction between independent variables was considered.The detailed steps of this kind of analysis are described in Crawley [88].According to this procedure, preliminarily to the regression analysis, two explorative investigations should be carried out: (1) a regression tree to see if the model has a complex structure, i.e., there are interaction terms between the independent variables; and (2) a generalized additive model (GAM) analysis to see if high order terms are necessary.Once ascertained the general model, a repeated OLS regression is performed cancelling, each time, the non-significant terms starting from the highest interaction terms.When all the surviving terms are significant, the regression procedure stops.However, if heteroscedasticity is present in the final model diagnostic, the regression procedure must be repeated from the beginning after the transformation of the dependent variable.In our case, the regression tree and the GAM analysis both indicated that we had to start with a model comprising all the main effects, i.e., all the independent band ratios, all the interaction terms and a quadratic term for the band ratio R rs (667)/R rs (488).At the end of the regression steps, the studentized Breusch-Pagan test indicated that the dependent variable had to be transformed, also in this case as y = y 0.5 .Repeating all the regression steps, the final model was the one with the dependent variable transformed as y = y 0.5 and with the following lengthy form: The adjusted coefficient of determination R 2 was 0.9873.The results of the regression analysis carried out by OLS are summarized in Table 7.

Model Validation
In order to validate the four models tested so far, the calculated values a CDOM 355 c of each model were compared to the in situ values of the validation data a CDOM 355 is (validation data set, 79 observations collected from 2015 to 2016).The calculated a CDOM 355 c were obtained by inserting the appropriate values of the band ratios of the validation data set into the equations of the four fitted models.The accuracy of the four models was then assessed by comparing the calculated values against the in situ values of the validation data set according to following statistical indicators: (a) the adjusted coefficient of determination (adjusted R 2 ) between a CDOM 355 is and a CDOM 355 c ; (b) the mean and standard deviation of the APD (APD = 100 × |a CDOM 355 is − a CDOM 355 c |/a CDOM 355 is ); and (c) the RMSE between a CDOM 355 is and a CDOM 355 c .Table 8 summarizes the results of the validation.For comparison, we have considered three additional empirical models developed by other authors.When a band of the additional models was not coincident with any of the bands considered for this work, the closest band among those at hand was chosen.The second column of Table 8 indicates a very good correlation between the calculated and the in situ values for the first three models, with adjusted R 2 values in the range of 0.832-0.843.The adjusted R 2 value for Model 4 was far below in spite of its highest fitting of the test data.The good agreement for the first three models must be truly appreciated considering that the validation data were collected more than one year after the test data.According to the APD values, Model 2 has the best prediction accuracy, while Model 4 the lowest one.The best value, 31 for Model 2, is comparable with APD's found by other authors [39,71] even if Mannino et al. obtained values as low as 18-20 [37].On the absolute scale, the prediction accuracy was almost comparable for the first three models, as indicated by the RMSE values.With the highest RMSE value, the low prediction capability of Model 4 was confirmed.The performances of the three additional models were quite low and comparable to Model 4, in some cases even worse.

Application of Satellite-Derived CDOM
In order to verify the ability of Model 2 (i.e., the one with the best prediction accuracy) to capture the dynamic range of CDOM for the test areas, maps of satellite a CDOM 355 were produced for the studied region.The selected days to produce the maps were chosen on the basis of daily flow rate of the Po River [89] and the availability of cloud free MODIS-A R rs imagery in order to show the distribution of riverine CDOM along the NCWAS.The processed MODIS-A satellite R rs images, as shown as examples in Figure 9, indicate that the proposed Model 2 algorithm can be confidently applied to generate consistent MODIS-A maps of a CDOM 355.The satellite retrieved CDOM distribution fits quite well based on the field measurements and the processes that affect the distribution of biochemical properties along the NCWAS.
Furthermore, the field measurements as well as optical satellite data demonstrate a decrease in a CDOM 355 from North (Po River mouth) towards South and from nearshore to offshore.Increasing and decreasing of CDOM due to differences in Po River flow rate is a quite evident process (as expected) along the Western coast of the Adriatic Sea.

Discussion
According to the above model analysis, the band ratio R rs (667)/R rs (488) was the most significant independent variable for the retrieval of a CDOM 355 in the study area of NCWAS, which is characterized by high river discharges.The results confirm that the performance of empirical algorithms in a complex basin, i.e., a basin strongly influenced by freshwater/terrestrial input, can be significantly improved by selecting at least one band with a wavelength greater than 600 nm [34,40,78].The additional spectral bands in the red or near infrared could be also helpful in better accounting for detritus particles because they are less sensible to the constituents in the seawater column so that the R rs (667) can be assumed to be mainly invariant as not influenced by detritus particles.On the sea surface at the wavelength of 667 nm the total light absorption is highest for seawater rather than for other constituents in the seawater column so that the R rs (667) can be assumed to be constant [34].Moreover, riverine CDOM shows a higher absorption at 488 nm than the chlorophyll [40].On the other hand, the R rs (547) nm is mostly driven by particles light scattering that in coastal water could reduce the uncertainty of CDOM retrieval when used as a reference point in band ratios [77].Furthermore, the NCWAS is well described by many authors to be impacted by particles coming from many rivers [60][61][62]90].Therefore, we applied the ratio R rs (667)/R rs (488) that is more sensitive to changes in CDOM rather than in chlorophyll because, in low salinity waters as in the NCWAS, most of the light attenuation is controlled by CDOM [14,91].The high river discharges on the NCWAS should therefore imply the use of a band ratio normally used in retrieving algorithms for river-plume waters.In this work, only simple linear models were considered.However, as shown in the Section 3, they proved to be quite effective for modelling the available data.Of the final four algorithms, only Model 2 and Model 4 fulfil the fundamental requirements for a correct linear regression.Considering the prediction capability, from our analysis Model 2 can be indicated as the best algorithm in order to retrieve a CDOM 355 in NCWAS starting from the band ratio R rs (667)/R rs (488).The predictions of this model in data collected more than one year after the model fitting showed a very good comparison with the in situ measured data with an adjusted R 2 value of 0.8322, Table 8.Just for comparison, similar validation analyses for coastal sea waters gave R 2 values of: 0.37 in the same geographical area for a CDOM 412 [71], 0.61 for a CDOM 443 in the western Canada coastal waters [92], 0.62 for a CDOM 412 in the Gulf of Mexico [91], 0.70 for a CDOM 442 in the North Sea and Western English Channel [32], 0.83 a CDOM 420 in a Finnish lake [93] and 0.96 for a CDOM 355, a CDOM 412 and a CDOM 443 in coastal waters in U.S. Middle Atlantic Bight [37].To better illustrate the agreement between calculated and measured values, in Figure 10 a graphical comparison between the in situ validation data and the Model 2 calculated data is shown.Figure 10 indicates that the prediction capability of the model seems to be rather good.It is however fair to add that the model performances could be overestimated due to the size of the validation data set, which is not very large, and the limited range of a CDOM 355 values, which roughly spans only one order of magnitude.
Plots of satellite-derived distributions of a CDOM 355 are spatially and temporally consistent with changes in Po River discharges (Figure 9) and the NCWAS.Since terrestrial DOM is the main source of CDOM to estuaries and coastal seas, a CDOM 355 values should generally decrease from the estuary to off-shore.The CDOM maps of Figure 9, moreover, show that the minimum CDOM concentration appears in 2015 only after a long period of reduced Po River flows (daily flow less than 500 m 3 •s −1 ), nevertheless determining the distribution of the sediments along the western Adriatic coast.The increase in surface a CDOM 355 corresponds to high river discharge.
As an example, Figures 9b and 11 show how an exceptional Po River discharge, happened in autumn 2014 (mean discharge 10 November-10 December 2014: 5028.2 m 3 •s −1 ) spreads high concentrations of riverine CDOM over the NCWAS compared with other periods, reaching the Istrian coast [44].In particular, Figure 11 shows the CDOM map obtained by applying Model 2 to the R rs Landsat 8 (30 m/pixel) bands closest to the R rs (667) and R rs (488) MODIS bands.In this case, we exceptionally tested the Model 2 equation on the Landsat R rs imagery, which was besides the severe cloud cover affecting the study area the only available for the day of the maximum Po River discharge between 2014 and 2016 (i.e., 19 November 2014; see plot of Figure 9e).Landsat retrieved CDOM map, because of the higher GSD, moreover, correctly depicts the detrital contribution of all the minor rivers' discharges occurring in the North and West Adriatic coasts, as well as the Venice Lagoon sub-basins.However, it is important to remark that MODIS R rs are recursively generated as standard product [68,69], while Landsat R rs were retrieved by applying the following commonly used R rs equation using as input Landsat radiance data.
where L tot (λ) = L 0 (λ) − L↑(λ) is the total spectral radiance signal above the water surface; L surf (λ) is the contribution to the radiance due to the reflection of the water surface; L w (λ) is the water leaving spectral radiance; and E d (λ) = E(λ)cosθ z τ ↓ (λ) is the incident irradiance above the water surface.Once the surface reflection is removed, the R rs can be considered Lambertian for short variation of the viewing angle around the zenith.In order to estimate and then remove the signal reflected by the water surface (L surf ) we used the Hydrolight 4.3 software that is a radiative transfer numerical model that computes radiance distributions and derived quantities for water bodies [94].
surface at the wavelength of 667 nm the total light absorption is highest for seawater rather than for other constituents in the seawater column so that the Rrs(667) can be assumed to be constant [34].Moreover, riverine CDOM shows a higher absorption at 488 nm than the chlorophyll [40].On the other hand, the Rrs(547) nm is mostly driven by particles light scattering that in coastal water could reduce the uncertainty of CDOM retrieval when used as a reference point in band ratios [77].Furthermore, the NCWAS is well described by many authors to be impacted by particles coming from many rivers [60][61][62]90].Therefore, we applied the ratio Rrs(667)/Rrs(488) that is more sensitive to changes in CDOM rather than in chlorophyll because, in low salinity waters as in the NCWAS, most of the light attenuation is controlled by CDOM [14,91].The high river discharges on the NCWAS should therefore imply the use of a band ratio normally used in retrieving algorithms for river-plume waters.In this work, only simple linear models were considered.However, as shown in the Section 3, they proved to be quite effective for modelling the available data.Of the final four algorithms, only Model 2 and Model 4 fulfil the fundamental requirements for a correct linear regression.Considering the prediction capability, from our analysis Model 2 can be indicated as the best algorithm in order to retrieve aCDOM355 in NCWAS starting from the band ratio Rrs(667)/Rrs(488).The predictions of this model in data collected more than one year after the model fitting showed a very good comparison with the in situ measured data with an adjusted R 2 value of 0.8322, Table 8.Just for comparison, similar validation analyses for coastal sea waters gave R 2 values of: 0.37 in the same geographical area for aCDOM412 [71], 0.61 for aCDOM443 in the western Canada coastal waters [92], 0.62 for aCDOM412 in the Gulf of Mexico [91], 0.70 for aCDOM442 in the North Sea and Western English Channel [32], 0.83 aCDOM420 in a Finnish lake [93] and 0.96 for aCDOM355, aCDOM412 and aCDOM443 in coastal waters in U.S. Middle Atlantic Bight [37].To better illustrate the agreement between calculated and measured values, in Figure 10 a graphical comparison between the in situ validation data and the Model 2 calculated data is shown.Figure 10 indicates that the prediction capability of the model seems to be rather good.It is however fair to add that the model performances could be overestimated due to the size of the validation data set, which is not very large, and the limited range of aCDOM355 values, which roughly spans only one order of magnitude.Plots of satellite-derived distributions of aCDOM355 are spatially and temporally consistent with changes in Po River discharges (Figure 9) and the NCWAS.Since terrestrial DOM is the main source of CDOM to estuaries and coastal seas, aCDOM355 values should generally decrease from the estuary to off-shore.The CDOM maps of Figure 9, moreover, show that the minimum CDOM concentration appears in 2015 only after a long period of reduced Po River flows (daily flow less than 500 m 3 •s −1 ), nevertheless determining the distribution of the sediments along the western Adriatic coast.The increase in surface aCDOM355 corresponds to high river discharge.
As an example, Figure 9b and Figure 11 show how an exceptional Po River discharge, happened in autumn 2014 (mean discharge 10 November-10 December 2014: 5028.2 m 3 •s −1 ) spreads high concentrations of riverine CDOM over the NCWAS compared with other periods, reaching the Istrian coast [44].In particular, Figure 11 shows the CDOM map obtained by applying Model 2 to the Rrs Landsat 8 (30 m/pixel) bands closest to the Rrs(667) and Rrs(488) MODIS bands.In this case, we exceptionally tested the Model 2 equation on the Landsat Rrs imagery, which was besides the severe cloud cover affecting the study area the only available for the day of the maximum Po River discharge between 2014 and 2016 (i.e., 19 November 2014; see plot of Figure 9e).Landsat retrieved CDOM map, because of the higher GSD, moreover, correctly depicts the detrital contribution of all the minor rivers' discharges occurring in the North and West Adriatic coasts, as well as the Venice Lagoon subbasins.However, it is important to remark that MODIS Rrs are recursively generated as standard product [68,69], while Landsat Rrs were retrieved by applying the following commonly used Rrs equation using as input Landsat radiance data.
Rrs(λ) = Ltot(λ) − Lsurf(λ)/Ed(λ) = Lw(λ)/Ed(λ) (8) where Ltot(λ) = L0(λ) − L↑(λ) is the total spectral radiance signal above the water surface; Lsurf(λ) is the contribution to the radiance due to the reflection of the water surface; Lw(λ) is the water leaving spectral radiance; and Ed(λ) = E(λ)cosθzτ↓(λ) is the incident irradiance above the water surface.Once the surface reflection is removed, the Rrs can be considered Lambertian for short variation of the viewing angle around the zenith.In order to estimate and then remove the signal reflected by the water surface (Lsurf) we used the Hydrolight 4.3 software that is a radiative transfer numerical model that computes radiance distributions and derived quantities for water bodies [94].Landsat atmospheric correction and R rs retrieval are anyhow out of the scope of this study and the use of Landsat data is not an important step for consistency validation of the CDOM algorithm here proposed and validated only for MODIS-A R rs on the NCWAS coastal areas.In addition, the central wavelengths of the two bands applied for the a CDOM algorithm (667 nm and 488 nm) from Landsat-8 and MODIS are close, while their bandwidths and SNRs (signal to noise ratio) are different.
The decrease in surface a CDOM 355 towards Southern Adriatic Sea and off-shore during spring 2014 and summer 2015 is consistent with photobleaching due to high solar radiation (both periods; [10,20]) as well as low river discharge (summer period).

Conclusions and Future Work
Based on the in situ data collected from several cruises in the North-Central Western Adriatic Sea (Mediterranean Sea) and on the corresponding MODIS ocean colour data, empirical models for the estimation of Coloured Dissolved Organic Matter (CDOM) using sensor data were tested.Data were split into two data sets: a test data set for model fitting (data collected from 2013 to 2015) and a validation data set (data collected from 2015 to 2016) for model assessment.The model set up was based on some statistical indicators: adjusted R 2 , absolute percentage difference (APD) and root mean square error (RMSE).
Following statistical check and diagnostic, four final simple linear models were obtained.Three of them gave very good results in terms of comparison between the calculated and the in situ measured values of the validation data set (adjusted R 2 in the range 0.83-0.84).On the base of APD and RMSE, the best model resulted to be a simple linear model (Model 2) with the dependent variable transformed according to y = y 0.5 and the band ratio R rs (667)/R rs (488) as the only independent variable.The use of the R rs (667) instead of the R rs (547) as a reference point in band ratio seems to reduce the uncertainty of CDOM retrieval in the NCWAS characterized by high light particles scattering.However, in this view, future work will be dedicated to demonstrate that retrievals are effective even for high particulate concentrations, e.g., by showing their efficacy even for high R rs (547) conditions and for particularly high river outflows.
The validation methods and processing of MODIS-A data show that Model 2 provides a good retrieval for a CDOM 355 and it is consistent with the spatial and temporal distribution of freshwater discharges along the Western Adriatic Sea.
The spatial and temporal limitations associated with in situ and remote sensing data has introduced some uncertainties.Anyway, the retrieved algorithm (Model 2) is very promising for CDOM mapping in this complex basin and should be considered when applying ocean colour remote sensing data to quantify other ocean constituents (e.g., chlorophyll).
MODIS-A data provision started on 2002 and its quality in terms of bands and Signal to Noise Ratio is going down, therefore, future work will be focused on developing new CDOM algorithms for the recently launched multispectral optical satellite sensors (e.g., Landsat 8 and Sentinel 2 and 3 with a GSD of 30 m) including also VIIRS observations and next generation hyperspectral sensor (e.g., EnMAP and PRISMA with a GSD of 30 m and the VNIR VENµS Superspectral Camera with a GSD of 10 m) and their validation strategy using novel in situ measurements.

Figure 1 .
Figure 1.Study area.In the upper pictures, the blue dotted box indicates Zone A, the green dotted box Zone B and the yellow dotted box Zone C (see text).In the top left picture, stations pertaining to the test data set are indicated by blue dots, in the top right picture the red dots indicate stations pertaining to the validation data set.Many stations were sampled more than once.In the lower picture, the red box indicates the location of the study area in the Mediterranean basin.

Figure 1 .
Figure 1.Study area.In the upper pictures, the blue dotted box indicates Zone A, the green dotted box Zone B and the yellow dotted box Zone C (see text).In the top left picture, stations pertaining to the test data set are indicated by blue dots, in the top right picture the red dots indicate stations pertaining to the validation data set.Many stations were sampled more than once.In the lower picture, the red box indicates the location of the study area in the Mediterranean basin.
Remote Sens. 2017, 9, 180 6 of 22 coastal area occurring in this date of acquisition (22 July 2015) and the higher values in the blue band at 443 nm for the stations more distant from the coast (i.e., 3, 10 and 11).

Figure 2 .
Figure 2. (A) Overview of MODIS Rrs imagery (R: 547 nm; G: 488 nm; B: 443 nm) with the Rrs spectra relative to the oceanographic cruises stations acquired in the Ancona coastal area (Italy) on 22 July 2015 shown in inset (B) as coloured dots (with the relative station number).Eastings and Northings are in geographic coordinates (latitude/longitude). (C) Rrs spectra of the eleven measurement stations shown in (B) and represented with the same colour and number of the dots.

Figure 2 .
Figure 2. (A) Overview of MODIS R rs imagery (R: 547 nm; G: 488 nm; B: 443 nm) with the R rs spectra relative to the oceanographic cruises stations acquired in the Ancona coastal area (Italy) on 22 July 2015 shown in inset (B) as coloured dots (with the relative station number).Eastings and Northings are in geographic coordinates (latitude/longitude). (C) R rs spectra of the eleven measurement stations shown in (B) and represented with the same colour and number of the dots.

Figure 3 .
Figure 3. Relationship between aCDOM355 and salinity (S) from in situ observations.The symbols represents the three areas of the Western Adriatic Sea, see text.The solid line is a regression line.Figure 3. Relationship between a CDOM 355 and salinity (S) from in situ observations.The symbols represents the three areas of the Western Adriatic Sea, see text.The solid line is a regression line.

Figure 3 .
Figure 3. Relationship between aCDOM355 and salinity (S) from in situ observations.The symbols represents the three areas of the Western Adriatic Sea, see text.The solid line is a regression line.Figure 3. Relationship between a CDOM 355 and salinity (S) from in situ observations.The symbols represents the three areas of the Western Adriatic Sea, see text.The solid line is a regression line.

Figure 4 .
Figure 4. Plot of the in situ aCDOM355 values versus the in situ band ratio Rrs(667)/Rrs(488) of the test data set.The solid and dashed lines are the OLS fitted Model 1 and the 95% confidence interval, respectively.

Figure 5 .
Figure 5. Plot of the residuals of Model 1 versus fitted values.Note the spread increase of the residual values with the increase of the fitted values.

Figure 4 .
Figure 4. Plot of the in situ a CDOM 355 values versus the in situ band ratio R rs (667)/R rs (488) of the test data set.The solid and dashed lines are the OLS fitted Model 1 and the 95% confidence interval, respectively.

Figure 4 .
Figure 4. Plot of the in situ aCDOM355 values versus the in situ band ratio Rrs(667)/Rrs(488) of the test data set.The solid and dashed lines are the OLS fitted Model 1 and the 95% confidence interval, respectively.

Figure 5 .
Figure 5. Plot of the residuals of Model 1 versus fitted values.Note the spread increase of the residual values with the increase of the fitted values.

Figure 5 .
Figure 5. Plot of the residuals of Model 1 versus fitted values.Note the spread increase of the residual values with the increase of the fitted values.

Figure 7 .
Figure 7. Plot of the in situ aCDOM355 0.5 values versus the in situ band ratio Rrs(667)/Rrs(488) of the test data set.The solid and dashed lines are the OLS fitted Model 2 and the 95% confidence interval, respectively.

Figure 6 .
Figure 6.Plot of the log-likelihood function as a function of λ for the variable transformation y = y α .

Model 2 :
aCDOM355 0.5 = b1 + b2 × Rrs(667)/Rrs(488) by OLS (5) Table5reports the regression results obtained for Model 2 on the test data set.The model fit is shown in Figure7along with the experimental data.The adjusted coefficient of determination was 0.8369.Both the plots of the residuals (not shown) and the studentized Breusch-Pagan test (p = 0.7571) indicate that the problem of heteroscedasticity was solved.

Figure 7 .
Figure 7. Plot of the in situ aCDOM355 0.5 values versus the in situ band ratio Rrs(667)/Rrs(488) of the test data set.The solid and dashed lines are the OLS fitted Model 2 and the 95% confidence interval, respectively.

Figure 7 .
Figure 7. Plot of the in situ a CDOM 355 0.5 values versus the in situ band ratio R rs (667)/R rs (488) of the test data set.The solid and dashed lines are the OLS fitted Model 2 and the 95% confidence interval, respectively.

Figure 8 .
Figure 8. Plot of the in situ aCDOM355 values versus the in situ band ratio Rrs(667)/Rrs(488) of the test data set.The solid and dashed lines are the GLS fitted Model 3 and the 95% confidence interval, respectively.

Figure 8 .
Figure 8. Plot of the in situ a CDOM 355 values versus the in situ band ratio R rs (667)/R rs (488) of the test data set.The solid and dashed lines are the GLS fitted Model 3 and the 95% confidence interval, respectively.

Figure 9 .
Figure 9. Satellite-derived aCDOM355 maps obtained using MODIS-A Rrs computed using the Model 2 equation (aCDOM355 0.5 = b1 + b2 × Rrs(667)/Rrs(488)) for the days: (A) 20 May 2014; (B) 11 December 2014; (C) 30 March 2015; and (D) 6 August 2015.The calculated maps matches pretty well the flow rates of the Po River (E).CDOM maps are overlaid (only for visualization purposes) on ASTER GDEM World Wide Elevation data (1 arc-second of spatial resolution) where dark green colour depicts land and the blue-violet colour the sea (in our case, no-data values due to the clouds).

Figure 9 .
Figure 9. Satellite-derived a CDOM 355 maps obtained using MODIS-A R rs computed using the Model 2 equation (a CDOM 355 0.5 = b 1 + b 2 × R rs (667)/R rs (488)) for the days: (A) 20 May 2014; (B) 11 December 2014; (C) 30 March 2015; and (D) 6 August 2015.The calculated maps matches pretty well the flow rates of the Po River (E).CDOM maps are overlaid (only for visualization purposes) on ASTER GDEM World Wide Elevation data (1 arc-second of spatial resolution) where dark green colour depicts land and the blue-violet colour the sea (in our case, no-data values due to the clouds).

Figure 10 .
Figure 10.Plot of the in situ aCDOM355 values versus the aCDOM355 values calculated according to Model 2. All data are from the validation data set.The solid line corresponds to the one-to-one relationship and the dashed lines are the 95% prediction intervals of Model 2.

Figure 10 .
Figure 10.Plot of the in situ a CDOM 355 values versus the a CDOM 355 values calculated according to Model 2. All data are from the validation data set.The solid line corresponds to the one-to-one relationship and the dashed lines are the 95% prediction intervals of Model 2.

Figure 11 .
Figure 11.CDOM (m −1 ) map retrieved by applying Model 2 equation to the Landsat 8 Rrs acquired on the NCWAS on 19 November 2016.CDOM map is overlaid (only for visualization purposes) on ASTER GDEM World Wide Elelvation data (1 arc-second of spatial resolution) where dark green colour depicts land and the light grey colour the sea (in our case no-data values due to the clouds).

Figure 11 .
Figure 11.CDOM (m −1 ) map retrieved by applying Model 2 equation to the Landsat 8 R rs acquired on the NCWAS on 19 November 2016.CDOM map is overlaid (only for visualization purposes) on ASTER GDEM World Wide Elelvation data (1 arc-second of spatial resolution) where dark green colour depicts land and the light grey colour the sea (in our case no-data values due to the clouds).

Table 1 .
List of field data from cruises for the absorption coefficient of Coloured Dissolved Organic Matter (a CDOM ) and salinity (S).

Table 2 .
Mean and standard deviation (s.d.) of in situ a CDOM 355 and salinity (S) measurements as a function of the area, see Figure1.The last column shows the number of observations (n).

Table 3 .
Results of the OLS regression of a CDOM 355 algorithm (3) on the test data set.

Table 4 .
Results of the OLS regression of Model 1 on the test data set.

Table 5 .
Results of the OLS regression of Model 2 on the test data set.

Parameter (Model 2) Parameter Estimate Standard Error t-Value p
* Significant at p = 0.05.

Table 6 .
Results of the GLS regression Model 3 on the test data.

Table 6 .
Results of the GLS regression Model 3 on the test data.
Parameter (Model 3) Parameter Estimate Standard Error t-Value* Significant at p = 0.05.

Table 7 .
Results of the OLS regression of Model 4 on the test data set.

Table 8 .
Comparison results of the prediction ability of the four models on the validation data set.Three additional models taken from bibliography are also included.