## Abstract

This study addresses the regression of in-water radiometric profile data with the objective of investigating solutions to minimize uncertainties of derived products like subsurface radiance and irradiance (*L*_{u0} and *E*_{d0}) and diffuse attenuation coefficients. Analyses are conducted using radiometric profiles generated through Monte Carlo simulations and field measurements. A nonlinear NL approach is presented as an alternative to the standard linear method LN. Results indicate that the LN method, relying on log-transformed data, tends to underestimate regression results with respect to NL operating on non-transformed data. The log-transformation is thus identified as the source of biases in data products. Observed differences between LN and NL regression results for *L*_{u0} are of the order of 1–2%, that is well below the target uncertainty for data products from *in situ* measurements (i.e., 5%). For *E*_{d0}, instead, differences can easily exceed 5% as a result of more pronounced light focusing and defocusing effects due to wave perturbations. This work also remarks the importance of applying the multi-cast measurement scheme as a mean to increase the precision of data products.

© 2013 OSA

## 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 ℜ_{0} and the diffuse attenuation coefficient *K*_{ℜ}[11].

Current operational processors [12, 13] determine ℜ_{0} 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 ℜ_{0} 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:

_{0}indicates the value just below the sea-surface and

*K*

_{ℜ}> 0 for

*z*positive downward is the diffuse attenuation coefficient. The exponential law is used to approximate Eq. 1 in a water layer were

*K*

_{ℜ}(

*z*) does not vary appreciably. After taking the logarithm, Eq. 2 becomes The linear LN method relies on Eq. 3 to compute regression parameters from radiometric values ℜ(

*z*) at different depths (

_{i}*i*= 1,..., N) by minimizing the sum-of-squares SSE error

_{0}and

*K*

_{ℜ}are obtained solving the linear system

*α*=log(ℜ

_{0}). Subsurface value and diffuse attenuation coefficient computed with the linear method are henceforth denoted by ${\Re}_{0}^{\text{LN}}$ and ${K}_{\Re}^{\text{LN}}$.

#### 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 ℜ

_{0}and

*K*

_{ℜ}by directly minimizing

*z*) values—cfr. Eq. 4 and 6.

_{i}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 [ℜ̂_{0}, *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 [ℜ̂_{0}, *K*̂_{ℜ}] called *trust region*. The approximation function *q* is so that its minimum [ℜ̃_{0}, *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}([ℜ̃_{0}, *K*̃_{ℜ}]) < SSE^{NL}([ℜ̂_{0}, *K*̂_{ℜ}]) then [ℜ̂_{0}, *K*̂_{ℜ}]=[ℜ̃_{0}, *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

^{−2}· nm

^{−1}. 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.

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^{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

^{−1}and

`a0d50`indicates

*a*= 0.50 m

^{−1}).

### 3.1.2. Virtual radiometric profiles

Virtual profiles have been generated for a deployment speed *v _{d}* of 0.25ms

^{−1}and a sampling frequency

*ν*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.

_{p}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^{−1} 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:

*j*-th radiometric cast within the set of virtual or

*in situ*samples. ℑ

^{NL}is chosen as the reference in Eq. 7 acknowledging that NL regression results are not biased by the log-transformation of radiometric data.

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

_{0}and

*K*

_{ℜ}are the results of the LN or the NL regression. The comparison between MSE

^{LN}and MSE

^{NL}values allows for verifying the performance of the TR algorithm, which might be affected by the presence of local minima. SSE values could be used for this purpose as well. However, by removing the dependence on the number of samples, MSE offers the additional capability of cross-checking results derived for different environmental cases or regression layers.

#### 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.

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

*θ*= 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

_{s}*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.

#### 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).

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.

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).

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

^{−1}(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)).

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

*p*(

*x|z*) is equal to 1/

*x*

_{M}being

*x*

_{M}the width of the simulation domain). A change of

*p*

_{dir}(

*θ|z*) can be induced by: 1) variations of the Fresnel refraction angle at elements defining the sea-air interface, which implies that photons refracted close to nadir tend to penetrate deeper in the water column; and 2) photon scattering, whose effect depends on the volume scattering phase function. Hence, the depth at which the diffuse attenuation becomes asymptotically constant is also a function of the beam attenuation and absorption coefficients. The heterogeneity of the distribution of the photon traveling direction is documented in Fig. 12(a) and 12(b) for

*θ*= 30° and 60°, respectively. Examples of

_{s}*p*

_{dir}(

*θ|x*,

*z*) at four spots selected in Fig. 12(b) are displayed in Fig. 13 to illustrate the cause of changes with depth of the diffuse attenuation coefficient for the considered case. A more detailed description of the joint effects of sea-surface geometry and the multiple photon scattering is out of the scope of the present study.

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^{−4} is an hypothetical sample and 10^{−6} 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:

_{cast}= 100 is the number of single casts.

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).

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.

## A. Weighted nonlinear regression WN

In the weighted least squares regression WN, an additional scale factor (i.e., the weight *w*) is included in the fitting process with respect to the standard nonlinear scheme. The sum of squared residuals is given by

*y*and

_{i}*ŷ*are the observed and the fitted response values, respectively. Fluctuations of the radiometric field decrease with depth, and the same reduction can be imposed to the

_{i}*w*weights. In fact, the weight factor should be inversely proportional to the squared deviation (variance) between average and observed radiometric values at each depth. The WN method hence allows for determining ℜ

_{i}_{0}and

*K*

_{ℜ}by giving higher “importance” to deeper data where light focusing effects tend to reduce. The variance of

*y*tends to the mean squared error of measured signal when depth grows and wave focusing effects becomes negligible. Hence, the WN method can also be extended to take into account the precision of radiometric measurements affecting profile data.

_{i}The following procedure computes the weighting factor as a function of depth:

- A first fitted response value
*ŷ*_{1}=_{i}*ŷ*(0) ·*e*^{−γ1·zi}is derived by applying the nonlinear minimization technique for summed square of residuals (*y*−_{i}*ŷ*)_{i}^{2}. - Approximation of square of residuals (
*y*−_{i}*ŷ*(0) ·*e*^{−γ1·zi})^{2}and application of the nonlinear minimization for $S={\sum}_{i=1}^{N}{\left[{\left({y}_{i}-\widehat{y}\left(0\right)\cdot {e}^{-{\gamma}_{1}\cdot {z}_{i}}\right)}^{2}-{\widehat{\epsilon}}^{2}\left({z}_{i}\right)\right]}^{2}$.

The nonlinear data regression can be further simplified when making the hypothesis that perturbations follow an exponential decay. In this case there are two unknown parameters in the expression for the fitted response value *ŷ _{i}* =

*ŷ*(0) ·

*e*

^{−}

*. The first one,*

^{γ·zi}*ŷ*(0) yields a linear dependence, the second one,

*γ*, correspond to the linear exponential component. Assuming that the nonlinear parameter

*γ*is known, the value of

*ŷ*(0) can be expressed as

*ŷ*is a nonlinear function of the single parameter

_{i}*γ*(the overline indicates averaged values). Hence, the problem of nonlinear minimization in the case of exponential law approximation can be reduced to a nonlinear optimization with only one unknown.

## Acknowledgments

This study has been partially supported by the European Space Agency within the framework of the MERIS Validation Activities under contract n. 12595/09/I-OL with the “Faculdade de Ciencias e Tecnologia” of the “Universidade Nova de Lisboa”, and contract n. 21653/08/I-OL with the Institute for Environment and Sustainability of the Joint Research Centre. The work also benefited of support from the NATO Science for Peace and Security Program, project n. 982678. Access to the Milipeia cluster (University of Coimbra, Portugal) was granted through the project “Large-scale parallel Monte Carlo simulations for Ocean Colour applications.”

## References and links

**1. **W. B. Philip and C. D. Scott, “Modelling regional responses by marine pelagic ecosystems to global climate change,” Geophys. Res. Lett. **29**, 1806 (2002).

**2. **J. L. Sarmiento, R. Slater, R. Barber, L. Bopp, S. C. Doney, A. C. Hirst, J. Kleypas, R. Matear, U. Mikolajewicz, P. Monfray, V. Soldatov, S. A. Spall, and R. Stouffer, “Response of ocean ecosystems to climate warming,” Global Biogeochem. Cycles **18**, 23 (2004) [CrossRef] .

**3. **M. J. Behrenfeld, R. T. ÓMalley, D. A. Siegel, C. R. McClain, J. L. Sarmiento, G. C. Feldman, A. J. Milligan, P. G. Falkowski, R. M. Letelier, and E. S. Boss, “Climate-driven trends in contemporary ocean productivity,” Nature **444**, 752–755 (2006) [CrossRef] [PubMed] .

**4. **B. A. Franz, S. W. Bailey, P. J. Werdell, and C. R. McClain, “Sensor-independent approach to the vicarious calibration of satellite ocean color radiometry,” Appl. Optics **46**, 5068–5082 (2007) [CrossRef] .

**5. **T. Kajiyama, D. D’Alimonte, and G. Zibordi, “Match-up analysis of MERIS radiometric data in the Northern Adriatic Sea,” IEEE Geosci. Remote Sens. Lett. (2013). Accepted for publication [CrossRef] .

**6. **G. Zibordi, F. Mélin, J.-F. Berthon, and E. Canuti, “Assessment of MERIS ocean color data products for European seas,” Ocean Sci. Discuss. **10**, 219–259 (2013) [CrossRef] .

**7. **T. Kajiyama, D. D’Alimonte, and G. Zibordi, “Regional algorithms for European seas: a case study based on MERIS data,” IEEE Geosci. Remote Sens. Lett. **10**, 283–287 (2013) [CrossRef] .

**8. **D. D’Alimonte, G. Zibordi, J.-F. Berthon, E. Canuti, and T. Kajiyama, “Performance and applicability of bio-optical algorithms in different European seas,” Remote Sens. Environ. **124**, 402–412 (2012) [CrossRef] .

**9. **S. Sathyendranath, “Remote sensing of ocean colour in coastal, and other optically-complex waters,” International Ocean-Colour Coordinating Group, IOCCG Report NUMBER 3 (2000).

**10. **S. Hooker, G. Zibordi, J.-F. Berthon, D. D’Alimonte, S. Maritorena, S. Mclean, and J. Sildam, *Results of the Second SeaWiFS Data AnalysisRound Robin, (DARR-00)*(NASAGSFC, Greenbelt, MD, USA, 2001), vol. 15 of *SeaWiFS Technical Report SERIES,* chap. 1, pp. 4–45.

**11. **G. Zibordi and K. Voss, *Field Radiometry and Ocean Colour Remote Sensing*(Springer, 2010), chap. 18, pp. 307–334.

**12. **P. J. Werdell and S. W. Bailey, “An improved in-situ bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. **98**, 122–140 (2005) [CrossRef] .

**13. **G. Zibordi, J.-F. Berthon, F. Mélin, and D. D’Alimonte, “Cross-site consistent in situ measurements for satellite ocean color applications: the BiOMaP radiometric dataset,” Remote Sens. Environ. **115**, 2104–2115 (2011) [CrossRef] .

**14. **J. R. Zaneveld, E. Boss, and P. Hwang, “The influence of coherent waves on the remotely sensed reflectance,” Opt. Express **9**, 260–266 (2001) [CrossRef] [PubMed] .

**15. **J. R. V. Zaneveld, E. Boss, and A. Barnard, “Influence of surface waves on measured and modeled irradiance profiles,” Appl. Optics **40**, 1442–1449 (2001) [CrossRef] .

**16. **G. Zibordi, D. D’Alimonte, and J.-F. Berthon, “An evaluation of depth resolution requirements for optical profiling in coastal waters,” J. of Atm. and Ocean. Tech. **21**, 1059–1073 (2004) [CrossRef] .

**17. **Y. You, D. Stramski, M. Darecki, and G. W. Kattawar, “Modeling of wave-induced irradiance fluctuations at near-surface depths in the ocean: a comparison with measurements,” Appl. Optics **49**, 1041–1053 (2010) [CrossRef] .

**18. **M. Hieronymi and A. Macke, “On the influence of wind and waves on underwater irradiance fluctuations,” Ocean Sci. **8**, 455–471 (2012) [CrossRef] .

**19. **M. Hieronymi, A. Macke, and O. Zielinski, “Modeling of wave-induced irradiance variability in the upper ocean mixed layer,” Ocean Sci. **8**, 103–120 (2012) [CrossRef] .

**20. **J. L. Muller and R. W. Austin, *Ocean Optics Protocols SeaWiFS for Validation, Revision 1*(NASA GSFC, Greenbelt, MD, USA, 1995), vol. 25 of *SeaWiFS Technical Report SERIES,* chap. 6, pp. 48–59.

**21. **D. A. Siegel, *Results of the SeaWiFS Data Analysis Round-Robin, July 1994 (DARR-94)*(NASA GSFC, Greenbelt, MD, USA, 1995), vol. 26 of *SeaWiFS Technical Report SERIES,* chap. 3, pp. 44–48.

**22. **J. J. Beauchamp and J. S. Olson, “Corrections for bias in regression estimates after logarithmic transformation,” Ecology **54**, 1403–1407 (1973) [CrossRef] .

**23. **D. D’Alimonte, G. Zibordi, T. Kajiyama, and J. C. Cunha, “Monte Carlo code for high spatial resolution ocean color simulations,” Appl. Optics **49**, 4936–4950 (2010) [CrossRef] .

**24. **T. Kajiyama, D. D’Alimonte, J. Cunha, and G. Zibordi, “*High-performance ocean color Monte Carlo simulation in the Geo-info project*,” in “Parallel Processing and Applied Mathematics,”, vol. 6068 of *Lecture Notes in Computer Science,*R. Wyrzykowski, J. Dongarra, K. Karczewski, and J. Wasniewski, eds. (Springer, 2010), vol. 6068 of Lecture Notes in Computer Science, pp. 370–379.

**25. **D. Schattschneider, “Proof without words: The arithmetic mean-geometric mean inequality,” Math. Mag. **59**, 11 (1986) [CrossRef] .

**26. **G. A. F. Seber and C. J. Wild, *Nonlinear regression*, Wiley series in probability and statistics (J. Wiley & Sons, 2003).

**27. **P. E. Gill, W. Murray, and M. H. Wright, *The Levenberg-Marquardt Method*(Academic Press, 1981), chap. 4.7.3, pp. 136–137.

**28. **Y. Yuan, “*A Review of Trust Region Algorithms for Optimization*,” in “ICIAM 99,” (Oxford University, 2000), pp. 271–282.

**29. **G. R. Fournier and J. L. Forand, “Analytic phase function for ocean water,” in “Ocean Optics XII,” (1994), no. 2558 in SPIE, pp. 194–201.

**30. **G. R. Fournier and M. Jonasz, “Computer-based underwater imaging analysis,” P. Soc. Photo-opt. Ins. **3761**, 62–70 (1999).

**31. **T. Kajiyama, D. D’Alimonte, and J. C. Cunha, “Performance prediction of ocean color Monte Carlo simulations using multi-layer perceptron neural networks,” (2011), vol. 4, pp. 2186–2195. Proceedings of the International Conference on Computational Science, ICCS2011.

**32. **T. Kajiyama, D. DAlimonte, and J. C. Cunha, “A high-performance computing framework for large-scale ocean color Monte Carlo simulations,” Concurrency Computat.: Pract. Exper pp. 1–22 (2011). Submitted for publications.

**33. **T. Kajiyama, D. D’Alimonte, and J. C. Cunha, “A statistical approach to performance tuning of Monte Carlo ocean color simulations,” in “Parallel and Distributed Computing, Applications and Technologies 2012,” (Beijing, China, 2012).

**34. **J. Escher and T. Schlurmann, “On the recovery of the free surface from the pressure within periodic traveling water waves,” J. Nonlinear Math. Phys. **15**, 50–57 (2008) [CrossRef] .

**35. **N. L. Jones and S. G. Monismith, “Measuring short-period wind waves in a tidally forced environment with a subsurface pressure gauge,” Limnol. Oceanogr.: Methods **5**, 317–327 (2007) [CrossRef] .

**36. **C. Tsai, M. Huang, F. Young, Y. Lin, and H. Li, “On the recovery of surface wave by pressure transfer function,” Ocean Eng. **32**, 1247–1259 (2005) [CrossRef] .