1. Introduction

Environmental monitoring and climate change studies require quality controlled ocean color products [1–3]. The fulfillment of such a demand often relies on the availability of field measurements of apparent optical properties (AOPs) for the vicarious calibration of space-born sensors [4] and the assessment of the related remote sensing radiometric data (e.g., [5, 6]). In situ AOPs are also essential to develop bio-optical algorithms for the generation of high-level products (e.g., pigment indices [7, 8]). The need of highly accurate in situ measurements is thus the rationale for continuous developments in measurement protocols and data reduction methods, mostly when considering applications in coastal regions challenged by the complexity of seawater bio-optical properties [9].

Processing codes for the determination of radiometric data products from in situ measurements apply schemes whose validity is sometimes restricted to ideal cases and whose implementation is often matter of subjective interpretations. This is documented by an inter-comparison of processing codes which exhibited a range of different results from the same in-water radiometric data despite computational methods were identical [10]. Recognizing that theoretical investigations can assist in understanding, quantifying, and possibly minimizing uncertainties affecting derived products, this study focuses on computational methods for the reduction of in-water radiometric profile data. Considered quantities are the upwelling radiance L_u and the downward irradiance E_d, henceforth generically denoted ℜ. Matter of investigation is the regression of radiometric measurements performed at different depths to derive the subsurface value ℜ₀ and the diffuse attenuation coefficient K_ℜ[11].

Current operational processors [12, 13] determine ℜ₀ and K_ℜ through the linear regression of log-transformed ℜ(z) versus depth z[11]. Solutions adopted for lessening perturbations induced by wave focusing and defocusing of light in radiometric measurements [14–19] include removing most perturbed samples, binning data in regular depth intervals and performing an incremental regression in successive layers [10, 15, 20]. Least-squares fits have also been applied to derive K_ℜ after interpolating log[ℜ(z)] with Hermitian cubic polynomials [21]. Since the log-transformation might bias results [22], this study evaluates nonlinear data reductions to compute ℜ₀ and K_ℜ directly from ℜ(z). Linear and nonlinear solutions are compared and discussed considering both virtual optical profiles generated through Monte Carlo (MC) simulations [23, 24] and actual radiometric profiles from field measurements [11, 13]. The aim is to investigate uncertainties of regression data products in the presence of light field perturbations due to surface waves. The removal of light focusing and defocusing effects in radiometric profile data and the assessment of related regression products is behind the objectives of this study.

The manuscript is organized as follows. Regression methods are described in Section 2. Simulated and in situ radiometric data are presented in Section 3. Results of Section 4 are further discussed in Section 5 through a set of case studies addressing: 1) the vertical distribution of perturbing effects; 2) changes in the distribution of the traveling direction of photons; and 3) the number of radiometric samples per unit depth. Summary and final remarks are in Section 6.

2. Methods

2.1. The linear method LN

The dependence of ℜ on z is expressed by:

2.2. The nonlinear method NL

The average of log-transformed values tends to be smaller than the logarithm of their average because of the inequality between the arithmetic and the geometric means [25]. This can bias regression results computed from log[ℜ(z_i)]. Corrections for bias minimization (e.g., [22]) are of difficult application because perturbations due to light focusing and defocusing largely vary in the water column. This work hence presents a nonlinear NL approach [26] to compute ℜ₀ and K_ℜ by directly minimizing

The minimization of SSE^NL requires a nonlinear solution because the partial derivatives of Eq. 6 do not lead to a linear system. Two numerical optimization schemes have been evaluated in this study: the Levenberg-Marquardt (LM) [27] and the Trust-Region (TR) [26, 28] algorithms. These algorithms produced equivalent results for most of the in situ and simulated radiometric profile data. However, in a few cases LM generated SSE^NL values much larger with respect to TR. The TR was then selected in virtue of an expected better reliability. The use of a weighted nonlinear scheme is also discussed in Section 5 in view of comprehensively addressing potentials of NL methods.

2.2.1. Trust-Region algorithm

The TR algorithm starts from an initial solution guess [ℜ̂₀, K̂_ℜ], here set to the LN result. A function q (e.g., second order Taylor expansion) is then used to approximate SSE^NL in a neighborhood Π of [ℜ̂₀, K̂_ℜ] called trust region. The approximation function q is so that its minimum [ℜ̃₀, K̃_ℜ] can be found easier than the minimum of SSE^NL.

The TR algorithm relies on the following iterative steps to satisfy convergence criteria: 1) if SSE^NL([ℜ̃₀, K̃_ℜ]) < SSE^NL([ℜ̂₀, K̂_ℜ]) then [ℜ̂₀, K̂_ℜ]=[ℜ̃₀, K̃_ℜ] and q is updated, otherwise the starting point remains the same and Π is shrunk; then 2) a new minimum of q is computed in Π. Additional details on the TR algorithm can be found in literature [26, 28]. ${\Re }_0^{\text{NL}}$ and $K_{\Re }^{\text{NL}}$ henceforth indicate the NL regression solution.

3. Data

3.1. Monte Carlo simulations

Radiometric data have been generated with the Monte Carlo (MC) code for Ocean Color Simulations MOX [23] accounting for a number of environmental conditions, including: 1) the sea surface wave length l and height h; 2) the sun-zenith angle θ_s; and 3) IOPs defined by the beam attenuation and absorption coefficients (c and a respectively). For all simulations, the above-water irradiance is E_s =1 in units of W · m⁻² · nm⁻¹. The seawater refractive index is set to 1.34 and the wavelength dependence is implicitly defined through the values of IOPs. The complete list of MOX simulation parameters is given in Table 1, whereas Fig. 1 presents a schematic of photon tracing.

Table 1. MOX simulation parameters. Inputs in brackets indicate sets of values applied in different simulations. The coefficients of the Fournier-Forand volume scattering function [29, 30] are: slope of the Junge distribution m =3.5835 and refraction index n =1.34.

View Table | View all tables in this article

 

Fig. 1 Schematics of MOX simulation domain and photon tracing (the depicted path can be referred equally to a photon starting from the sun or from the sky). The radiometric sensor of a virtual free-fall profiler is represented by the “collecting bins” located at the nodes of the 2D grid. Simulations presented in this study were performed with resolution and collecting bin size of 1 cm, over a domain 50 m deep and 10 m width.

Download Full Size | PDF

The MOX code has been optimized for speed by applying high-performance computing solutions on supercomputers [24, 31–33] to trace the number of photons (e.g., up to 5 × 10¹⁰ for the present work) required to produce results with a precision analogous to that of field measurements [23].

3.1.1. Simulation case-naming

A specific naming convention is applied to identify MOX simulation cases (Table 2). Taking the r30c0d6a0d50 case-name as an example: 1) the first letter corresponds to the sea-surface model ( r for a single harmonic component and q for two harmonic components); 2) the following two-digit number indicates the sun elevation (i.e., 30 means θ_s = 30°); and 3) seawater attenuation and absorption coefficients are at the end (i.e., c0d6 indicates c = 0.6 m⁻¹ and a0d50 indicates a = 0.50 m⁻¹).

Table 2. MOX case-names with corresponding simulation parameters. Taking the r30c0d6a0d50 case-name as an example: 1) the first letter corresponds to the sea-surface model ( r for a single harmonic component and q for two harmonic components); 2) the following two-digit number indicates the sun elevation (i.e., 30 means θ_s = 30°); and 3) seawater attenuation and absorption coefficients are at the end (i.e., c0d6 indicates c = 0.6 m⁻¹ and a0d50 indicates a = 0.50 m⁻¹).

View Table | View all tables in this article

3.1.2. Virtual radiometric profiles

Virtual profiles have been generated for a deployment speed v_d of 0.25ms⁻¹ and a sampling frequency ν_p of 6 Hz, assuming instantaneous measurements of depth and radiometric values. The linear wave theory [34–36] has been applied to model the response of the pressure transducer and compute an adjusted vertical position of each virtual data point.

The comparison between NL and LN results is based on the following data sets: 1) one hundred single-casts, each made by a single virtual profile; 2) one hundred multi-casts, each formed by five profiles randomly sampled from the single-casts pool; and 3) one all-cast, containing all 100 single-casts. The reference extrapolation layer to derive regression parameters is between 0.25 and 5 m (e.g., Panel 2(b)), as suggested by real conditions encountered in coastal regions. Additional extrapolation intervals considered in the study are [0.25, 15] and [5, 15] m.

3.2. In-situ measurements

With particular reference to optically complex waters, the determination of highly accurate in-water radiometric products often requires profiling near the surface to minimize the perturbing effects of non-homogeneities in the IOP vertical distribution, and to account for the rapid decrease of light signal with depth in highly attenuating waters. An additional requirement is the capability of producing a number of measurements per unit depth, not significantly affected by the tilt of the profiling system, to lessen the effects of wave perturbations [16]. As a consequence, the precision of derived subsurface radiometric products is a function of the sampling depth-interval and of the depth resolution defined by the system acquisition rate and deployment speed.

In situ radiometric profiles of E_d and L_u used in this study were collected during an oceanographic cruise performed in the Black Sea in 2012. The system for in-water optical profiling was equipped with OCR-507 (Satlantic Inc., Halifax) radiance and irradiance sensors with spectral channels of 10 nm bandwidth centered at the nominal wavelengths 412, 443, 490, 510, 555, 665, and 683 nm. In-water profiles were performed with 6 Hz acquisition rate and ∼ 0.3 m s⁻¹ deployment speed. Expected uncertainties in derived remote sensing reflectance R_RS are ∼4 – 5% in the 412–555 nm spectral range and ∼ 6% at 665 nm. Radiometric products and their determination are discussed in [13, 16]. Measurements suitable for the determination of subsurface values were collected in the near surface layer by applying the multi-cast technique to increase the number of samples per unit depth [16, 23]. Specifically, data from multiple casts performed within a time frame of 10 minutes were grouped to form a single radiometric profile with an increased depth resolution.

4. Results

4.1. Statistical figures

The percent difference between LN and NL regression results, with ℑ generically indicating the subsurface value or the diffuse attenuation coefficient, is computed as:

Statistical figures for the comparison of LN and NL results are the mean difference μ_ℑ and standard deviation of the differences σ_ℑ, respectively:

Data fitting is also evaluated through the Mean Square Error (MSE) between radiometric samples and regression values

4.2. Simulation results

MOX simulations displayed in Figs. 2–3 illustrate radiometric fields and optical profiles in the case of sea-surface described by a single harmonic component. Results for two harmonic components are in Figs. 4–5. Statistical figures in Table 3(a) indicate that the values of $E_{\text{d}0}^{\text{LN}}$ ( $K_{E_{\text{d}}}^{\text{LN}}$) tend to be lower than those of $E_{\text{d}0}^{\text{NL}}$ ( $K_{E_{\text{d}}}^{\text{NL}}$). This is explained by the log-transformation for the LN data regression in conjunction with the fact that light focusing and defocusing effects generally decrease with depth (Figs. 2 and 5). Differences between $E_{\text{d}0}^{\text{LN}}$ and $E_{\text{d}0}^{\text{NL}}$ appear larger in the regression layers [0.25, 15] and [5, 15] m with respect to [0.25, 5] m. δ_Ed0 tends also to increase with the sun elevation. The differences between $K_{E_{\text{d}}}^{\text{LN}}$ and $K_{E_{\text{d}}}^{\text{NL}}$ show instead a limited dependence on both depth intervals and sun zenith angles. LN results larger than those obtained from the NL regression are discussed in Section 5.

 

Fig. 2 MOX simulations illustrating E_d radiometric fields and optical profiles in the case of a single harmonic component at the sea-surface (l = 5 m and h = 0.5 m). Results for sun-zenith angles θ_s = 30° and θ_s = 60° are in the top and bottom row panels, respectively.

Download Full Size | PDF

 

Fig. 3 As in Fig 2, but with two harmonic components at the sea-surface (l₁ = 5 m, h₁ = 0.5 m and l₂ = 0.5 m, h₂ = 0.05 m).

Download Full Size | PDF

 

Fig. 4 As in Fig 2, but for L_u.

Download Full Size | PDF

 

Fig. 5 As in Fig 4, but with two harmonics at the sea-surface.

Download Full Size | PDF

Table 3. LN vs. NL regression results expressed by the mean difference μ_ℑ (Eq. 8) in the case of simulated data. Standard deviation values (σ_ℑ, Eq. 9) are within brackets.

View Table | View all tables in this article

Examples of δ_Ed0 and δ_KEd (Eq. 7) frequency distributions are displayed in Fig. 6. Column panels correspond to different MOX simulations. A single wave harmonic is used to model the sea-surface. Results derived from the [0.25, 5] m layer indicate that LN tends to overestimate with respect to NL when θ_s = 30° (first two row panels). An opposite tendency is seen for θ_s = 60° (bottom row panels). Single-cast and multi-cast results represented by white and black bars, respectively, confirm that an increased number of samples per unit depth permits to improve the precision of regression results with respect to the single-cast case [16]. The agreement between LN and NL data products is better for L_u (Table 3(b)) than for E_d. The reason is that L_u is less affected by light focusing and defocusing effects (cfr. Fig 2 and 4), hence $L_{\text{u}0}^{\text{LN}}$ is also less subject to biases induced by log-transformation. The regression layer exhibiting lower μ_Lu (i.e., usually less than 1%) is [0.25, 5] m. The frequency distributions of differences presented in Fig. 7 confirm the importance of the multi-cast scheme also for the regression of L_u data.

 

Fig. 6 Frequency distribution of δ_Ed0 and δ_KEd values derived from MOX simulations in the case of a single wave component at the sea-surface. The comparison between NL and LN results is based on the following data partitioning: 1) one hundred single-casts, each made by single virtual profile; 2) one hundred multi-casts, each formed by five individual profiles randomly sampled from the single-casts pool; and 3) one all-cast, containing all 100 single-casts.

Download Full Size | PDF

 

Fig. 7 As in Fig 6, but for L_u.

Download Full Size | PDF

4.3. In situ measurement results

In situ radiometric profiles applied in this study are presented in Fig. 8 and 9. Panel insets indicate the above-water irradiance E_s corresponding to in-water samples. Each multi-cast is identified by the name of the corresponding field station. The depth intervals for the regression of E_d and L_u data are [0.25, 5] m and [0.25, 2] m, respectively (with the exception of the station k0790 for which the interval [0.3, 3] m is used for both radiometric quantities).

 

Fig. 8 In-situ optical profiles.

Download Full Size | PDF

 

Fig. 9 Continued from Fig 8.

Download Full Size | PDF

The comparison of LN and NL regression results is detailed in Table 4(a) and 4(b) for E_d and L_u, respectively. Rows indicate statistical values for data pertaining to a specific multi-cast set. Columns correspond to different wavelengths.

Table 4. LN vs. NL regression results expressed by the mean difference μ_ℑ (Eq. 8). Analysis based on in situ radiometric profiles.

View Table | View all tables in this article

The analysis of in situ measurements exhibits $E_{\text{d}0}^{\text{LN}}$ lower than $E_{\text{d}0}^{\text{NL}}$. An underestimate of $K_{E_{\text{d}}}^{\text{LN}}$ is also seen with respect to $K_{E_{\text{d}}}^{\text{NL}}$. A substantial equivalence of regression results is observed for $L_{\text{u}0}^{\text{LN}}$ and $L_{\text{u}0}^{\text{NL}}$, as well as for $K_{L_{\text{u}}}^{\text{LN}}$ and $K_{L_{\text{u}}}^{\text{NL}}$ (i.e., generally within 1%).

5. Discussion

5.1. Vertical distribution of perturbing effects

Both virtual and in situ radiometric profiles indicate that LN tend to produce E_d0 and K_Ed values smaller than those computed with the NL method. Negative μ_Ed0 but null μ_KEd would be expected if perturbations induced by light focusing and defocusing effects were multiplicative factors sampled from the same distribution within the entire regression layer (i.e., homoscedastic upon log-transformation). Negative μ_KEd values obtainable from radiometric profiles in clear waters can be explained by a reduction of perturbations with depth. However, a peculiar case is represented by E_d values from the r30c0d6a0d50 MOX simulation (Table 2). The corresponding irradiance field and an example of virtual multi-cast profile are displayed in Fig. 2 (i.e., Panel 2(a) and 2(b), respectively). The sea surface is modeled by a single harmonic component and the sun elevation is 30°. E_d samples selected for this analysis are detailed in Fig. 10, where the same data are presented on linear and log-scale in the left and right column panels, respectively. Regression results for the [0.25, 5], [0.25, 15] and [5, 15] m layer are from the top to the bottom row panels. Each regression layer is highlighted in gray. The green color and the ⊕ symbol indicate the LN regression line and $E_{\text{d}0}^{\text{LN}}$. The corresponding NL results are presented in red and with the ⊗ symbol. Insets in the log-scale panels report regression results (i.e., E_d0 and K_Ed), while those in the linear scale panels show data product differences (Eq. 7) and MSE values (Eq. 10).

 

Fig. 10 E_d case study based on the r30c0d6a0d50 MOX simulation. The same data are presented on linear and log scale in the left and right column panels, respectively. Regression results for the [0.25, 5], [0.25, 15] and [5, 15] m layer are from the top to the bottom row panels.

Download Full Size | PDF

Perturbations due to light focusing and defocusing are larger at the bottom of the [0.25, 5] m regression layer (see Panel 10(a) and 10(b)) and this leads to $K_{E_{\text{d}}}^{\text{LN}}>K_{E_{\text{d}}}^{\text{NL}}$. The value of E_d0 is also slightly larger than that of $E_{\text{d}0}^{\text{NL}}$. Panels 10(c) and 10(d) show results for the [0.25, 15] m regression layer. In this case, perturbations are larger in the upper part of the regression layer, and statistical figures of the Panels 10(c) inset indicate that $E_{\text{d}0}^{\text{LN}}<E_{\text{d}0}^{\text{NL}}$. The comparison of Panel 10(a) with 10(c) also shows that NL produces almost identical results in the [0.25, 5] m and the [0.25, 15] m layers (i.e., $E_{\text{d}0}^{\text{NL}}=0.994$ for the former and $E_{\text{d}0}^{\text{NL}}=0.995$ for the latter). Hence, the NL method performs well in fitting data of the [0.25, 5] m layer significantly affected by wave perturbations. A larger variation is instead observed for the $E_{\text{d}0}^{\text{LN}}$ values. Results from the [5, 15] m layer are shown in Panels 10(e) and 10(f). Regression values from the NL method better match those obtained in the other two layers when compared to results from the LN scheme. Additionally, MSE^NL < MSE^LN values highlight the better fitting performance of the NL method. It is finally highlighted that results with $K_{E_{\text{d}}}^{\text{LN}}>K_{E_{\text{d}}}^{\text{NL}}$ are not documented by the in situ data included in this work (Table 4). Nevertheless, the case formerly discussed represents a realistic condition although not typical.

5.2. Distribution of the photon traveling direction

This second analysis concerns L_u profile data. Results are presented in Fig. 11, where the layout of panels is analogous to that of Fig. 10. The selected MOX simulation is r60c1d0a0d20 with sea surface described by a single harmonic component, θ_s = 60°, c = 1.0 and a = 0.20 m⁻¹ (see Fig. 4 and Table 2). When compared to E_d, the lower perturbations due to the effects of light focusing and defocusing on L_u increase the convergence between the LN and NL solution—cfr. Panel 11(a) and Panel 10(a) for regression results determined in the [0.25, 5] m layer. The convergence between LN and NL results decreases in the [5, 15] m layer and even more in the [0.25, 15] m layer. Indeed, LN and NL regression lines of Panel 11(a) indicate a slope change in L_u data, with larger K_Lu values above ∼ 8 m depth. Note that IOPs are constant and there is no significant radiance contribution due to the sea-bottom reflectance (Figs. 4(k) and 4(l)).

 

Fig. 11 As in Fig. 10, but considering L_u data from the r60c1d0a0d20 MOX simulation.

Download Full Size | PDF

The diffuse attenuation would also be constant if the iso-depth distribution of the photon direction

 

Fig. 12 Average photon traveling direction. Results for θ_s =30° and θ_s =60° in Panel 12(a) and 12(b), respectively (c = 1.0 and a = 0.20 m⁻¹ in both cases). Details for spots labeled in the right panel are displayed in Fig. 13.

Download Full Size | PDF

 

Fig. 13 Examples of probability density (adimensional) of photon direction at position selected in Fig. 12. Both x and z are in units of m. The nadir is at θ =180°.

Download Full Size | PDF

Panel 11(c) shows that the LN regression line tends to be evenly driven by samples above and below the ∼ 8 m depth, here considered as a proxy to identify the level from which the diffuse coefficient becomes stationary. Small residuals between regression and sampled values become much larger when considering log-transformed quantities. For instance, if 10⁻⁴ is an hypothetical sample and 10⁻⁶ is the corresponding regression value, then: ${\text{SSE}}_i^{\text{LN}}={\left[-4+6\right]}^2$ and ${\text{SSE}}_i^{\text{NL}}={\left[{10}^{-4}-{10}^{-6}\right]}^2$, see Eq. 4 and 6, respectively. This residual is about eight orders of magnitude larger for SSE^LN than for SSE^NL. The NL line can hence overlap better with samples closer to the surface, as documented by the fact that $K_{L_{\text{u}}}^{\text{LN}}$ is 5.5% lower than $K_{L_{\text{u}}}^{\text{NL}}$, and $L_{\text{u}0}^{\text{LN}}$ underestimates $L_{\text{u}0}^{\text{NL}}$ of 8.0% (Panel 11(d)). In comply with results of Section 5.1, this analysis indicates that the agreement between regression results obtained in the [0.25, 15] and [0.25, 5] m layers is better for the NL than for the LN method. These findings also illustrate how differences between LN and NL regression results may depend on factors contradicting the assumption of validity of the exponential decay law (see Eq. 1 and 2) because of a peculiar combination of IOP values, illumination conditions and perturbations induced by light focusing and defocusing.

5.3. Number of samples per unit depth

This last discussion concerns the effect of the profile data density N_d expressing the number of samples per unit depth, upon the reproducibility of LN and NL regression results [16, 23]. This analysis also includes results obtained by applying a weighted nonlinear regression scheme WN. In summary, this scheme uses the residuals between NL regression and sampled values to define a weighted SSE and recompute regression parameters. By this mean, WN results should in principle be less dependent on data samples affected by larger perturbations (details in Appendix A).

The reference term for this analysis is the NL subsurface value ${\widehat{\Re }}_0^{\text{NL}}$ derived from the samples of the all-cast profile (Section 3.1.2) comprised in the [0.25, 5]m layer. The corresponding sub-surface value computed considering only the j-th single-cast and the N_d data density is denoted as ${\left[{\Re }_0^{\text{X}}\left(N_{\text{d}}\right)\right]}_j$, where X indicates the regression methods. In analogy with Eq. 7–9, the percent difference ${\widehat{\delta }}_j^{\text{X}}\left(N_{\text{d}}\right)$, the mean difference μ̂^X(N_d) and the standard deviation of differences σ̂^X(N_d) are:

Top panels of Fig. 14 report E_d results from MOX simulations with different sun elevations (θ_s = 30° and 60° on the left and right, respectively). Bottom panels show results for different sea-surface geometries (one and two harmonics on the left and right, respectively). The same layout is applied for L_u data of Fig. 15. LN, NL and WN results are represented by green, red and blue lines, respectively. Overall, simulation cases include 4 IOP configurations, 2 sun zenith values, and 2 sea-surface geometries. Each panel hence contains 8 curves of the same color (although some lines overlap).

 

Fig. 14 E_d0 variability as a function of the number of samples per unit depth. Top panels report results for θ_s = 30° and 60° on the left and right, respectively. Bottom panels show results for a sea-surface with one and two harmonics on the left and right, respectively.

Download Full Size | PDF

 

Fig. 15 As in Fig 14, but for L_u.

Download Full Size | PDF

Panels of Fig. 14 indicate that a larger sun elevation, as well as two harmonics at the sea surface, tend to increase the N_d value required to allow for a reproducibility of subsurface estimates below 5%. This is explained by an increase of the vertical gradient of perturbations due to light focusing in the [0.25, 5] m layer (see Fig. 2 and 3). It is also noted that in these cases LN results tend to have a better reproducibility, but also a larger bias (see Table 3). As an example, the q60c0d6a0d50 simulation case where σ̂^LN is lower than σ̂^NL corresponds to a bias δ_Ed0 ∼−15% (results are analogous when considering WN instead of NL). Fig. 14 and 15 also indicate that the complex structure of perturbations due to light focusing and defocusing limits the capability of the WN method to substantially increase the precision of subsurface values with respect to NL.

The documented small bias between $L_{\text{u}0}^{\text{LN}}$ and $L_{\text{u}0}^{\text{NL}}$ values computed with samples in the [0.25, 5] m layer (Section 4) conforms with results of Fig. 15 indicating the reduced number of data points needed to satisfy a 5% threshold with respect to the E_d case. This is in line with previous findings about depth resolution requirements for optical profiling in coastal waters [16].

6. Summary and conclusions

This study has addressed the reduction of in-water radiometric profiles with the objective of investigating solutions minimizing uncertainties of regression results like subsurface radiance and irradiance values (E_d0 than L_u0) and diffuse attenuation coefficients (K_Ed and K_Lu). Analyses have been conducted using radiometric profile data generated through Monte Carlo simulations and field measurements. A nonlinear NL approach has been investigated as an alternative to the standard linear LN method.

The study demonstrates the effectiveness of Monte Carlo radiative transfer codes in supporting investigations requiring accurate modeling of the physical and geometrical properties of realistic sea surfaces This work is then an additional example of the importance of virtual environments built on the basis of high-performance computing solutions, to assist the analysis of actual data and processing methods.

Results indicate that the LN method, relying on log-transformed data, tends to underestimate regression results with respect to NL which operates directly on the input radiometric values. The log-transformation is thus identified as a source of bias in data products. Observed differences between $L_{\text{u}0}^{\text{LN}}$ and $L_{\text{u}0}^{\text{NL}}$ are of the order of 1 – 2%, that is well below the target uncertainty for data products from in situ measurements (i.e., 5%). Instead, difference between $E_{\text{d}0}^{\text{LN}}$ and $E_{\text{d}0}^{\text{NL}}$ can easily exceed 5%. Difference between LN and NL regression results are much larger for E_d than L_u because the bias due to the log-transformation depends on the amplitude and the vertical distribution of perturbations induced by light focusing and defocusing. This supports the use of the NL method for the reduction of downward irradiance profile data.

Nonlinear methods such as the Trust Region algorithm, cannot ensure that the NL solution correspond to a global minimum of the error function. The NL method should then be applied accounting for this constraint, and it might be worth recomputing NL regression results using slightly different starting points to verify the convergence of the minimization scheme. It is also suggested to check that the MSE value is lower for the NL than for the LN method.

Study results were obtained from the analysis of simulated radiometric profile data characterized by deep water and vertically homogeneous IOPs. It is likely that the inclusion of bottom effects and optical stratifications in simulated data would further increase the complexity of the problem by imposing limitations to the width of the extrapolation layer. This may then lead to increased data products uncertainties due to a lower number of regression samples. It is however expected the overall conclusions about the LN and NL approaches are consistent with to those discussed here. It is finally highlighted that the precision of radiometric data products derived from the extrapolation of in-water profile data mainly depends on the instrument specifications and the applied measurement protocol (e.g, sampling frequency and deployment speed). The main requirement ensuring repeatability of regression results is a number of radiometric samples per unit depth sufficient to statistically represent in-water radiometric perturbations. This work hence further remarks the importance of applying the multi-cast measurement scheme.