## 2 Data and Methods
### Pigment Retrieval from \(R_{\text{rs}}(\lambda)\)
#### 2.5.1 \(a_{\text{Gau}}(\lambda)\) from \(R_{\text{rs}}(\lambda)\)
In the determination of Gaussian parameters (\(\sigma_{i}\) and \(\lambda_{i}\)) for \(a_{\text{ph}}(\lambda)\) in the global scale, we also tested various combinations of parameters using data published in the literature [63; 64; 65]. The existing parameters were not successful at obtaining satisfactory results for every data range due to various reasons, including the fact that the initial datasets used to obtain the parameters had a small dynamic range unable to cover varied conditions such as coastal regions and non-bloom natural oceanic waters. Thus, a refinement of the parameters \(\sigma_{i}\) and \(\lambda_{i}\) was conducted using the \(a_{\text{ph}}(\lambda)\) dataset obtained from SeaBASS to improve the overall performance of the MuPI model for global oceans. The non-linear least square fitting procedure in MATLAB was used to solve Equation (1). A set of refined Gaussian parameters for \(\sigma_{i}\) and \(\lambda_{i}\) were obtained and are presented in Table 2. For oceanic waters, as the absorption coefficient from water molecules contributes \(>\)80% of the total absorption coefficient for wavelengths \(>\)550 nm, it is difficult to obtain accurate \(a_{\text{Gau}}(\lambda)\) by directly inverting \(R_{\text{rs}}(\lambda)\) in the longer wavelengths. On the other hand, since the Gaussian peaks at 434 and 492 nm cover the main absorption features of the different pigments, the two \(a_{\text{Gau}}(\lambda)\) at 434 and 492 nm were chosen as the independent variables in this effort. Following [PERSON] et al. [73], the implementation of this model used two Gaussian peak heights [\(a_{\text{Gau}}(\lambda_{1})\) and \(a_{\text{Gau}}(\lambda_{2})\)] to reconstruct \(a_{\text{ph}}(\lambda)\), in which empirical relationships as shown in Table 2 were used. These relationships between \(a_{\text{Gau}}(\lambda)\) were obtained by regression analysis with the purpose of reducing the unknowns in the \(R_{\text{rs}}(\lambda)\) inversion procedure [63; 73]. With this design, there will be five unknowns [\(a_{\text{Gau}}(\lambda_{1})\), \(a_{\text{Gau}}(\lambda_{2})\), \(b_{\text{bp}}(\lambda_{0})\), \(a_{\text{dg}}(\lambda_{0})\) and \(S\)] to be retrieved from a \(R_{\text{rs}}(\lambda)\) spectrum, which is obtained by a minimization of the cost function (Equation (8)):  
\[\delta=\frac{\sqrt{\frac{1}{N_{\lambda}}\sum\limits_{i=1}^{N_{\lambda}}\left( \hat{R}_{\text{rs}}(\lambda_{i})-R_{\text{rs}}(\lambda_{i})\right)^{2}}}{\frac {1}{N_{\lambda}}\sum\limits_{i=1}^{N_{\lambda}}R_{\text{rs}}(\lambda_{i})} \tag{8}\]  
with \(N_{\lambda}\) as the wavelength number, \(R_{\text{rs}}(\lambda)\) as the measured, and \(\hat{R}_{\text{rs}}(\lambda)\) the modeled spectrum, respectively. Basically, \(\delta\) value provides a measure of the relative difference between the input and output \(R_{\text{rs}}\) spectra. The generalized reduced gradient (GRG) nonlinear optimization procedure [79] was used to solve Equation (8). The statistical indices used to estimate the agreement between the two values (\(\hat{R}\) and \(R\)) was the unbiased absolute percentage difference (UAPD), defined as Equation (9) and root mean square error (RMSE, Equation (10)) with \(N\) as the number of samples. \[\text{UAPD}=\frac{\left|\hat{R}-R\right|}{0.5\left(\hat{R}+R\right)}\times 1 00\% \tag{9}\]  
\[R_{\text{rs}}(\lambda)=0.52{\sum_{i=1}^{2}g_{i}\bigg{[}\frac{b_{\text{b}}( \lambda)}{a(\lambda)+b_{\text{b}}(\lambda)}\bigg{]}^{i}}/\Bigg{\{}1-1.7{\sum_{i=1 }^{2}g_{i}\bigg{[}\frac{b_{\text{b}}(\lambda)}{a(\lambda)+b_{\text{b}}(\lambda) }\bigg{]}^{i}\Bigg{\}} \tag{2}\]0