## 2 Methods
### Multivariate Optimization Procedure  
For a vector of concentrations \(C\) and a given hydro-optical model the residual between measured and reconstructed reflectance at \(j\)-th wavelength can be computed as follows:  
\[g_{j}=\left(S_{j}-T_{j}\right)/\text{ }T_{j} \tag{4}\]  
The multidimensional least-square solution using all wavelengths is found by minimizing the sum of squared residuals by simultaneous variation of concentrations in the vector \(C\):  
\[f\left(\textbf{{C}}\right)=\sum_{j}g_{j}^{2}\left(\textbf{{C}}\right) \tag{5}\]  
The absolute minimum of \(f\)_(C)_ can be found with the Levenberg-Marquardt finite difference algorithm [4; 5], which assures a rapid convergence of the iterative procedure. The following iterative expression is used to this goal:\[\mathbf{C}_{k+1}=\mathbf{C}_{k}+\zeta_{k}\left(F_{k}^{\,\prime}F_{k}+\mu_{k}D_{k}\right)^ {-1}F_{k}^{\,\prime}\!\left(1-\frac{T(\mathbf{C}_{k})}{\left\{S_{k}\right\}}\right) \tag{6}\]  
where \(k=\) iteration step, \(D_{k}=\)_diag_(\(F_{k}^{\,\prime}F_{k}\)) = diagonal matrix consisting of the elements \(F_{k}^{\,\prime}F_{k}\), with \(F(\mathbf{C})=\left|\frac{\partial T_{j}}{\partial\mathbf{C}_{j}}\right|=\) matrix composed of \(n\) x \(m\) elements (\(n=\) number of wavelengths, at which the measurements of \(S_{j}\) have been conducted, \(m=\) dimension of the concentration vector \(\mathbf{C}\)), \(F^{\,\prime}(\mathbf{C})=\) transposed matrix, \(\mu_{k}=\) direction of minimization, \(\zeta_{k}=\) length of the minimization step. The criterion for the termination of the iterative procedure is \(f(\mathbf{C})\leq 10^{-5}\). The actual length, \(\zeta_{k}\) can be chosen easily and effectively through applying the method of reduction, which consists in the following: at a certain iteration step \(k\)  a deliberately high value of \(\zeta_{k}\) is chosen (e.g. 50). From eq. 3 a new concentration vector \(\mathbf{C}_{k+1}\) is determined, for which the inequality \(f(\mathbf{C}_{k+1})<f(\mathbf{C}_{k})\) is tested. If it is not satisfied, \(\zeta_{k}\) is reduced by a factor of two, and a new round of search is initiated. This iterative procedure is continued until the above inequality is satisfied, i.e. until the newly determined concentration vector proves to be closer to the minimum of the function of residuals \(f(\mathbf{C})\) as compared to the concentration vector established in the previous iteration step. In the course of the iterative descent to the global minimum of the function of residuals, the above procedure could end up in a _local_ minimum. If the determined concentration vector, \(\mathbf{C}_{k}\)  happens to be on the downhill and the length of the iteration step \(\zeta\)  is sufficiently high, then the next concentration vector \(\mathbf{C}_{k+1}\) might occur on the downhill in the direction to a deeper minimum. In this case the process of search will move to a new point provided the inequality \(f(\mathbf{C}_{k+1})<f(\mathbf{C}_{k})\) is satisfied. Thus, such a search with a variable length \(\zeta\) can, under favourable conditions, evade a local _minor_ minimum and end either in a _deeper_ or even in the _global_ minimum. However, equally possible is the situation that the concentration vector \(\mathbf{C}_{k+1}\) is on the _uphill_ side of a certain local minimum of \(f\)_(\(\mathbf{C}\))_. In this case the search ends in this local minimum, which corresponds to an inadequate solution of the inverse problem. Naturally, the probability for the latter option increases with measurement errors. To avoid the above difficulty, it is advisable to initiate the iteration procedure for an array of initial guess values, \(\mathbf{C}_{\theta}\). Afterwards the deepest minimum of \(f\) (\(\mathbf{C}\)) is selected. The number \(\mathbf{N}\) of the initial vectors \(\mathbf{C}_{\theta}\) should not be excessively high because the computation time for the inverse problem solution increases proportionally with \(\mathbf{N}\).