Machine learning for estimating phytoplankton size structure from satellite ocean color imagery in optically complex Pacific Arctic waters

In response to recent advances in satellite ocean color remote sensing, we have developed a chlorophyll a size distribution (CSD) model using machine learning (ML) approaches for optically complex Pacific Arctic waters. Previous CSD models have used principal component analysis (PCA) to retrieve spectral features from satellite-estimated phytoplankton absorption coefficient (a_ph(λ)) by assuming a strong correlation between the spectral features and phytoplankton size structure determined from the exponent of the CSD (η). A weakness of such approach is that it relies on satellite retrievals of a_ph(λ), which can be highly uncertain due to the optical effects of water constituents other than phytoplankton. In this study, we have developed a method based on ML to use remote sensing reflectance (R_rs(λ)) for directly retrieving η, thus avoiding uncertainties due to the inversion of a_ph(λ) from R_rs(λ). Results show superior performance of the ML-based CSD models compared to the PCA-based model utilizing both R_rs(λ) and a_ph(λ) as predictors of η. For direct R_rs(λ)-based retrievals, a CSD model based on multivariable linear regression produced the best performance among all models considered. Nevertheless, models using in-situ a_ph(λ) yielded better accuracy, reflecting a closer optical linkage between η and a_ph(λ) than between η and R_rs(λ). Our choice of an R_rs(λ)-based model for satellite application is therefore practical, motivated by the limitations and uncertainty of a_ph(λ) inversions in optically complex waters. Another key finding is that more complex ML approaches do not always produce more effective models than standard linear regression. Indeed, multivariable linear regression outperformed other ML approaches for retrieving η directly from R_rs(λ), whereas support vector machine performed the best among diverse ML approaches in the case of a_ph(λ). Overall, this study found benefits in using R_rs(λ) with ML to improve the retrieval accuracy of η for Pacific Arctic waters.

1 Introduction

Satellite remote sensing is a cost-effective tool that can provide observations across a range of temporal and spatial scales. One of the primary parameters retrieved from ocean color satellite data is the mass concentration of chlorophyll a (Chl a; see Table 1 for symbols, definitions, and units), the primary pigment associated with photosynthesis and a key indicator of phytoplankton biomass. Satellite-derived Chl a observations have revolutionized our understanding of climate systems, marine ecosystems, and biogeochemical processes (McClain, 2009). However, Chl a alone does not provide a full description of the fundamental ecosystem functions of phytoplankton, such as nutrient uptake and cycling, energy transfer through marine food webs, deep-ocean carbon export, and gas exchange with the atmosphere (Mouw et al., 2017).

Due to the significance of phytoplankton community composition in ocean biogeochemical processes, continuous research and innovation in satellite ocean color techniques have extended our capabilities from routinely estimating Chl a concentration to retrieving phytoplankton functional types (PFTs) (Gordon et al., 1980; Mouw et al., 2017). PFTs are conceptual groupings of phytoplankton species that have similar biogeochemical functions (e.g., nitrogen fixers, calcifiers, dimethylsulfide producers, and silicifiers) and other characteristics such as cell size (pico-, nano-, and micro-phytoplankton). PFTs are often defined based on phytoplankton size class (PSC), phytoplankton taxonomic composition (PTC), or particle size distribution (PSD), and the choice of partitioning depends on the question at hand (Mouw et al., 2017), with no universally accepted standard (Reynolds et al., 2002). In particular, PSC serves as a useful index of the trophic state, carbon export efficiency, and productivity (Hood et al., 2006; Le Quéré et al., 2005) and, therefore, comprises the majority of PFT research.

A wide range of satellite-based methods for global estimations of PFTs have been developed to date (IOCCG, 2014). Mouw et al. (2017) provide a “user guide” for applying remote sensing techniques to monitor PFTs, explaining details of various PFT algorithms and their associated uncertainties and discussing the advantages and disadvantages of different approaches. Satellite estimation of PFTs generally exploit spectral features in remote sensing reflectance (R_rs(λ)), absorption coefficient of phytoplankton (a_ph(λ)), and/or backscattering coefficient of particles (b_bp(λ)) caused by variations in PFT composition (Fujiwara et al., 2011; Kostadinov et al., 2010; Li et al., 2013; Roy et al., 2017). The ocean color variables used in these spectral approaches are grouped into two categories: apparent optical properties (AOPs, e.g., R_rs(λ)) and inherent optical properties (IOPs, e.g., a_ph(λ)). Remotely sensed IOPs are derived from spectral inversion of R_rs(λ) (Mobley, 1994), thereby introducing additional uncertainties for IOP-based methods compared to R_rs(λ)-based methods.

For global estimation of PSC, Waga et al. (2017) developed a Chl a size distribution (CSD) model that retrieves the synoptic size structure of the phytoplankton community by determining the exponent of CSD (CSD slope; η). As opposed to other PSC approaches, η represents the size structure of the phytoplankton community with a single value; thus, the output of the approach can be easily incorporated into ocean biogeochemical models. Akin to the PSD (Kostadinov et al., 2010; Roy et al., 2017), the arbitrariness of the arrangement of the size range is another advantage of this approach, where other methods generally adopt a fixed target group or size class (e.g., <2, 2–20, and >20 µm for pico-, nano-, and micro-phytoplankton, respectively). More specifically, once η is determined, fractional contributions of phytoplankton biomass at diverse size ranges can be estimated from η. Moreover, there is flexibility in computing η with different combinations of size-fractionated Chl a, generating a comparable variable across datasets that often comprise various size ranges of size-fractionated Chl a data.

The spectral features of a_ph(λ) can reveal specific information regarding variations in the composition and size structure of phytoplankton assemblage (Bricaud and Morel, 1986a). For example, how pigments are distributed within a phytoplankton cell affects the magnitude of a_ph(λ), while pigment composition influences the spectral shape of a_ph(λ). Waga et al. (2017) applied principal component analysis (PCA) to normalized a_ph(λ) spectra derived from in situ measurements at seven wavelengths (412, 443, 469, 488, 531, 547, and 555 nm) that are consistent with spectral bands of the Moderate Resolution Imaging Spectroradiometer (MODIS). This method assumes that PCA captures spectral features of a_ph(λ) as a simpler set of principal component (PC) scores while still maintaining significant patterns and trends. The relationship between η and the resulting PC scores was then quantified by ordinary least squares regression, enabling η to be estimated from satellite derivations of a_ph(λ) (Waga et al., 2017). In order to investigate spatiotemporal variations in the size structure of phytoplankton communities and its impacts on the marine ecosystems in the Pacific Arctic, the CSD model was subsequently optimized for the Pacific Arctic based on a regional in situ dataset (Waga et al., 2019a). However, in Arctic coastal waters, phytoplankton absorption is typically low (only 16 % of non-water absorption at 443 nm) relative to colored dissolved organic matter (CDOM) and non-algal particles (NAP) and, as a result, IOP inversion algorithms for estimating a_ph(λ) are characterized by high uncertainty (Matsuoka et al., 2007). Therefore, direct approaches to estimate η utilizing R_rs(λ) may be advantageous in Arctic coastal environments, even though R_rs(λ) itself is not solely influenced by phytoplankton.

The present study develops the CSD model for the Pacific Arctic utilizing diverse supervised machine learning (ML) approaches, ranging from simple linear regression to convoluted methods such as neural networks (Chen et al., 2015, 2018; Li et al., 2020, 2023; Waga et al., 2022), support vector machines (Deng et al., 2019; Selvaraju et al., 2021; Su et al., 2015), Gaussian processes (Pasolli et al., 2010), and ensemble methods (Bao et al., 2023; Qi et al., 2022; Qiao et al., 2022; Zhang et al., 2023). A main advantage of ML is the ability to parameterize general relationships from training data without predefined or explicit equations (Marzban, 2009). To date, a variety of ML models have been used for retrieval of various ocean parameters, including the diffuse attenuation coefficient (Chen et al., 2015), particle backscattering coefficient (Sauzède et al., 2016), Chl a concentration (Chen et al., 2021; Hu et al., 2021; Kolluru and Tiwari, 2022; Mukonza and Chiang, 2022; Syariz et al., 2020), and reconstructions of ocean color data (Chen et al., 2019; Fasnacht et al., 2022; Krasnopolsky et al., 2016). The current study aims to (1) parameterize CSD models for the Pacific Arctic using spectral features of R_rs(λ) and a_ph(λ), (2) assess satellite algorithm performance using an in situ dataset, and (3) compare newly developed models with the previously developed PCA-based CSD model. The updated CSD model provides accurate estimates of spatiotemporal variations in PSC in the Pacific Arctic, providing key information on how recent environmental changes are affecting the foundation of marine food webs in a changing Arctic.

Table 1 Definitions and units of all symbols used in the text, figures, and equations.

2 Material and methods

An updated CSD model is proposed in this study to enable reasonable estimation of spatiotemporal variations in PSC for optically complex Pacific Arctic waters. See Sects. S1–S4 in the Supplement for complete materials and methods.

2.1 In situ data

Multiple research cruises were conducted in the Pacific Arctic during the summer months from 2007 to 2021 (Table 2). A total of 177 open ocean and coastal sampling locations were visited in the sub-Arctic Bering Sea and the west Beaufort Sea, including the Stefansson Sound near Prudhoe Bay along the northern coast of Alaska (Fig. 1). A companion map, color-coded by cruise year, is provided in Fig. S1 in the Supplement. At each station, spectral radiometric measurements were made during daylight hours, and water samples were collected for a_ph(λ) and size-fractionated Chl a (hereafter referred to as a_{ph_obs}(λ) and Chl a_{size_obs}, respectively).

Figure 1 Sampling locations of in situ data used in this study. Colors of each plot indicate the exponent of chlorophyll a (Chl a) size distribution (CSD slope; η_obs), whereas background color represent the bathymetry.

Table 2 Details of cruises, number of samples (N) obtained during each cruise, and filter pore sizes used to collect size fractionated chlorophyll a samples. Note that the cruise period indicates the date span of in situ data collected.

2.1.1 Phytoplankton pigments

Chl a_{size_obs} was determined using a 10 AU fluorometer (Turner Designs), except for ten samples from the 2021 cruise in Prudhoe Bay (PB21), for which Chl a_{size_obs} was determined using high performance liquid chromatography (HPLC). HPLC analysis provides the concentration of not only Chl a but also other major phytoplankton pigments (i.e., fucoxanthin, peridinin, 19^′-hexanoyloxyfucoxanthin, 19^′-butanoylofucoxanthin, alloxanthin, chlorophyll b, neoxanthin, prasinoxanthin, violaxanthin, lutein, and zeaxanthin). At each station in all the cruises, both fractionated and unfractionated (i.e., without filtration using filters of different pore sizes for size fractionation) samples were collected. Unfractionated HPLC samples were collected at each station in all the cruises.

2.1.2 Absorption coefficient

Particles in surface seawater samples were collected on a GF/F filter until the filter had sufficient coloration to measure a_{ph_obs}(λ). The absorption coefficient of particles (a_{p_obs}(λ)) on the filter was measured in the spectral range from 300 to 850 nm at 1 nm intervals using an MPS-2400 (Shimadzu Corporation), MPS-2450 (Shimadzu Corporation) or Cary 100 (Agilent Technologies) spectrophotometer. The quantitative filter technique (QFT) was used to determine a_{ph_obs}(λ) for samples measured with the MPS-2400 and MPS-2450 instruments (i.e., all cruises but PB21), following the procedure described by Mitchell (1990), whereas a_{ph_obs}(λ) for the PB21 samples was determined with GF/F filters placed inside a 15 cm integrating sphere connected to the Cary 100 (IOCCG, 2018). Following the measurement for a_{p_obs}(λ), the absorption coefficient of NAP (a_{NAP_obs}(λ)) was measured after soaking the filter in 95 % methanol or sodium hypochlorite, and a_{ph_obs}(λ) was finally obtained by subtracting a_{NAP_obs} (λ) from a_{p_obs} (λ). The absorption coefficient of CDOM (a_{CDOM_obs}(λ)) at wavelengths from 250 to 750 nm at 1 nm intervals was measured using the same spectrophotometers as for the particulate absorption measurements, with the exception of the PB21 samples, which were analyzed using a Cary 300 (Agilent Technologies) spectrophotometer with 5 cm quartz cuvettes.

2.1.3 Remote sensing reflectance

In situ spectral radiance and irradiance measurements were acquired using a PRR-800/810 (Biospherical Instruments), C-OPS (Biospherical Instruments), or HyperPro (Satlantic) spectroradiometer. Each spectroradiometer has different spectral resolutions and ranges: the PRR-800/810 and C-OPS collected at 17 (380 to 765 nm) and 19 wavelengths (320 to 875 nm), respectively, whereas the HyperPro acquired data between 400 and 800 nm at approximately 3 nm intervals. Remote sensing reflectance (R_{rs_obs}(λ)) was calculated as the ratio of the water-leaving radiance (L_w(λ)) to the above-water downward spectral irradiance (E_s(λ)):

$$\begin{array}{}\text{(1)}& R_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)=L_{\mathrm{w}}\left(\mathit{\lambda }\right)/E_{\mathrm{s}}\left(\mathit{\lambda }\right).\end{array}$$

R_{rs_obs}(λ) was resampled at ten MODIS bands in the visible range (i.e., 412, 443, 469, 488, 531, 547, 555, 645, 667, and 678 nm) from the original wavelengths of each instrument using spline interpolation (Wang et al., 2015). Finally, a modified version of the Quasi-Analytical Algorithm (QAA; Lee et al., 2002) for the Pacific Arctic (Fujiwara et al., 2016) was used to estimate a_ph(λ) (a_{ph_QAA}(λ)) from in situ R_rs(λ) (R_{rs_obs}(λ)) and satellite R_rs(λ) (R_{rs_sat}(λ)). Here, a_{ph_QAA}(λ) estimated from R_{rs_obs}(λ) and R_{rs_sat}(λ) is denoted as a_{ph_QAAobs}(λ) and a_{ph_QAAsat}(λ), respectively. To avoid the retrieval of negative a_{ph_QAA}(λ), the modified version of QAA uses an optimized spectral slope of the absorption coefficient of combined CDOM and NAP (S_dg) obtained by reconstructing the S_dg based on a dataset collected in the Pacific Arctic (Fujiwara et al., 2016). The a_{ph_QAAobs}(λ) was used to validate the performance of the modified version of the QAA by comparing it with a_{ph_obs}(λ).

2.1.4 Pigment-based identification of phytoplankton taxonomic composition

An open-source R software package, phytoclass (version 1.0.0), was used to determine the Chl a biomass of different phytoplankton groups from their accessory pigments (Hayward et al., 2023). The phytoclass package is a Chl a taxonomic partitioning software package similar to the widely used CHEMTAX software (Mackey et al., 1996). However, phytoclass has been shown to be more accurate and does not rely on initial assumptions of pigment to Chl a ratios for each phytoplankton group (Hayward et al., 2023). Eight target taxonomic groups (diatoms, chrysophytes, dinoflagellates, prymnesiophytes, chlorophytes, prasinophytes, cryptophytes, and cyanobacteria) and 11 marker pigments for each taxonomic group (peridinin, 19^′-butanoyloxyfucoxanthin, fucoxanthin, 19^′-hexanoyloxyfucoxanthin, neoxanthin, prasinoxanthin, violaxanthin, alloxanthin, lutein, zeaxanthin, and chlorophyll b) were selected following Zhuang et al. (2016), as these groupings have been used previously for CHEMTAX analysis in the Chukchi Sea shelf region.

2.2 Satellite data

The MODIS sensor onboard NASA's Aqua satellite (MODIS-A), operational since 2002, provides the longest time series among all currently operational ocean color sensors, which is an attractive advantage for decadal-scale monitoring and retrospective analyses. Level-3 standard mapped images of 4 km spatial resolution monthly climatological R_{rs_sat}(λ) at ten bands in the visible range (i.e., 412, 443, 469, 488, 531, 547, 555, 645, 667, and 678 nm) and daytime sea surface temperature (SST) derived by MODIS-A (version R2022.0) were downloaded from NASA's Ocean Color website. The R_{rs_sat}(λ) data were then used to compute a_{ph_QAAsat}(λ) by using the modified QAA algorithm (Fujiwara et al., 2016).

2.3 Chlorophyll a size distribution model

The exponent of the CSD (η), representing the size structure of phytoplankton communities, was determined following the method of Waga et al. (2017). Assuming the CSD follows a Junge-type power law distribution, the total Chl a (Chl a_total) and Chl a_size in a size range from D₁ to D₂ can be expressed as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{Chl}{a}_{\mathrm{total}}=\underset{D_{\min }}{\overset{D_{\max }}{\int }}\mathrm{Chl}{a}_{\mathrm{0}}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\eta }}\mathrm{d}D,\text{(3)}& {\displaystyle }\mathrm{Chl}{a}_{\mathrm{size}}=\underset{D_{\mathrm{1}}}{\overset{D_{\mathrm{2}}}{\int }}\mathrm{Chl}{a}_{\mathrm{0}}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\eta }}\mathrm{d}D,\end{array}$$

where Chl a₀ is the Chl a at a reference diameter D₀ (here, 0.7 µm). In this study, D_min and D_max were defined as 0.7 and 200 µm, respectively. η was derived as the slope of the linear regression in log-space computations between the inverse log-transformed median diameters (from D₁ to D₂), and Chl a_size normalized by the bin width. An advantage of the CSD model is its robustness when using different sets of Chl a_size to retrieve η (Waga et al., 2017).

A large η indicates a greater contribution of smaller-sized phytoplankton, whereas a small η suggests that larger-sized phytoplankton dominate. The fraction of Chl a_size can be derived using η as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{F}_{\mathrm{size}}& ={\displaystyle \frac{\mathrm{Chl}{a}_{\mathrm{size}}}{\mathrm{Chl}{a}_{\mathrm{total}}}}={\displaystyle \frac{\underset{D_{\mathrm{1}}}{\overset{D_{\mathrm{2}}}{\int }}\mathrm{Chl}{a}_{\mathrm{0}}{\left(\frac{D}{D_{\mathrm{0}}}\right)}^{-\mathit{\eta }}\mathrm{d}D}{\underset{D_{\min }}{\overset{D_{\max }}{\int }}\mathrm{Chl}{a}_{\mathrm{0}}{\left(\frac{D}{D_{\mathrm{0}}}\right)}^{-\mathit{\eta }}\mathrm{d}D}}\\ & ={\displaystyle \frac{D_{\mathrm{2}}^{\mathrm{1}-\mathit{\eta }}-D_{\mathrm{1}}^{\mathrm{1}-\mathit{\eta }}}{{\mathrm{200}}^{\mathrm{1}-\mathit{\eta }}-{\mathrm{0.7}}^{\mathrm{1}-\mathit{\eta }}}}\end{array}\end{array}$$

In this study, the size ranges for pico-, nano-, and micro-phytoplankton were defined as 0.7–2, 2–20, and 20–200 µm, respectively. To estimate the fraction of Chl a within the size ranges for pico- (F_pico), nano- (F_nano), and micro-phytoplankton (F_micro), D₁ and D₂ in Eq. (4) were set as the lower and upper limits of each size range. For clarification purposes, the size fractions determined from in situ Chl a_size observations are denoted as F_{size_obs}, whereas those estimated through a CSD model with Eq. (4) using η were represented as F_{size_MDL}.

2.4 Model development

The CSD model was trained using 70 % of the entire dataset (i.e., training subset), randomly determined using the MATLAB randsample function (R2025b), while the remaining 30 % was used for final validation (i.e., validation subset). The details of model development based on the PCA and supervised ML approaches are described in Sect. 2.4.1 and 2.4.2, respectively.

2.4.1 PCA approach

The previous version of the CSD model for the Pacific Arctic (Waga et al., 2019a) used the spectral shape of a_ph(λ) to estimate η. To capture the spectral features of a_ph(λ), PCA was applied to normalized a_{ph_obs}(λ) (${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$) at ten MODIS-A bands. The formula for ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ is:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)& =[a_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)\\ & -\mathrm{mean}\left(a_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)\right)]/\mathrm{SD}(a_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right))\end{array}\end{array}$$

where mean(a_{ph_obs}(λ)) and SD(a_{ph_obs}(λ)) are the arithmetic mean and standard deviation calculated across the ten wavelengths for each individual a_{ph_obs}(λ) spectrum, respectively. The input values for the PCA comprise a matrix (m×N) composed of ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ values, where m and N are the number of the wavelengths and number of samples, respectively. Assuming the resulting PC scores correlate with η, η was estimated as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\eta }={\left[{\mathit{\beta }}_{\mathrm{0}}+\exp \sum _{i=\mathrm{1}}^k{\mathit{\beta }}_i{S}_i\right]}^{-\mathrm{1}}\text{(7)}& {\displaystyle }{S}_i=\sum _{j=\mathrm{1}}^m{w}_{i,j}{\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left({\mathit{\lambda }}_j\right),\end{array}$$

where S_i and w_i,j are the ith PC score and the loading factors for ith PC at wavelength j. In addition, m and k represent the number of wavelengths and the number of PCs (k=4 in this study). The model parameters β₀ and β_i are the regression coefficients between η and PC scores.

By substituting for the calculation of S_i in Eq. (6), we obtained new equations as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathit{\eta }={\left[{\mathit{\beta }}_{\mathrm{0}}+\exp \sum _{j=\mathrm{1}}^m{C}_j{\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left({\mathit{\lambda }}_j\right)\right]}^{-\mathrm{1}}\text{(9)}& {\displaystyle }{C}_j=\sum _{i=\mathrm{1}}^k{\mathit{\beta }}_i{w}_{i,j},\end{array}$$

where β₀ and C_j are the final model parameters. Once the model parameters were determined based on a_{ph_obs}(λ), the same coefficients were used in the case of a_{ph_QAAobs}(λ) and a_{ph_QAAsat}(λ) to produce estimates of η. For the R_rs-based models, normalized R_{rs_obs}(λ) (${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$) was calculated in the same manner as Eq. (5), and η was determined by employing ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ in Eqs. (6)–(9) in place of ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. Note that η determined by Chl a_{size_obs}, estimated through the CSD model using in situ measurements ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ or ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and satellite products (${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{QAAsat}}\left(\mathit{\lambda }\right)$ or ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{sat}}\left(\mathit{\lambda }\right))$ are denoted as η_obs, and η_MDLobs and η_MDLsat, respectively.

2.4.2 Supervised ML approach

In addition to the PCA approach used in prior work (Waga et al., 2017, 2019a, b, 2021b), CSD models were trained with various ML approaches. Since we know both the input (i.e., ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ or ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$) and corresponding output (i.e., η_obs) values, supervised ML was used to train CSD models. To this end, we leveraged the Regression Learner App in the MATLAB Statistics and Machine Learning toolbox, a user-friendly resource that enables simple data exploration, feature selection, specification of validation schemes, model training, and model evaluation. This application includes commonly used regression methods, e.g., linear regression models, regression trees, Gaussian process regression models, support vector machines, kernel approximation models, ensembles of regression trees, and neural network regression models.

To avoid the possibility of missing certain representative samples and/or overfitting the models, repeated five-fold cross-validation (ten repeats) was carried out by randomly dividing the training subset into five equally sized sets (or five-folds). Evaluation of the trained models was performed five times, each time excluding one-fold from the training subset and using it for validation. Each observation in the training subset was assigned to an individual group and stayed in that group for the duration of the procedure so that each observation was allowed to be used one time for testing and four times for training the model. Finally, the performance of the trained models was determined as the average of the performance metrics from the five iterations.

The MATLAB Regression Learner App returns three other statistical metrics besides the coefficient of determination (r²): the root mean square error (RMSE), mean squared error (MSE), and mean absolute error (MAE) between the observed and predicted values, defined as:

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathrm{RMSE}=\sqrt{\sum _{n=\mathrm{1}}^N{\left(X_n-Y_n\right)}^{\mathrm{2}}/N}\text{(11)}& {\displaystyle }\mathrm{MSE}=\sum _{n=\mathrm{1}}^N{\left(X_n-Y_n\right)}^{\mathrm{2}}/N\text{(12)}& {\displaystyle }\mathrm{MAE}=\sum _{n=\mathrm{1}}^N\left|X_n-Y_n\right|/N\end{array}$$

where X_n and Y_n represent the nth observed and predicted values, respectively. Once CSD models based on each ML method were finalized, the best ML-based CSD model for each predictor (i.e., ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$) was determined based on the four aforementioned statistical metrics. Once the best-performing models for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ among diverse regression methods were determined, they were used for final validation and further analysis.

2.5 Model validation metrics

The performance of the resulting PCA-based CSD models and the best-performing ML-based CSD models were compared using the validation subset. Bias is a key metric for the performance assessment of satellite products (Seegers et al., 2018), defined as:

$$\begin{array}{}\text{(13)}& \mathrm{Bias}={\mathrm{10}}^{\left({\sum }_{n=\mathrm{1}}^N\left(X_n-Y_n\right)/N\right)}\end{array}$$

Following recommended validation procedures for satellite ocean color algorithms (Seegers et al., 2018), the performance of the CSD models, as well as the modified QAA, was evaluated based on MAE (Eq. 12) and bias.

3 Results

3.1 Phytoplankton size structure and taxonomic composition

The measured exponent of CSD (η_obs) values ranged from 0 to 2.24 with corresponding Chl a_{total_obs} values of 18.84 and 0.05 mg m⁻³, respectively (Table 3). Figure 2 depicts the Chl a_{total_obs} and η_obs values with regard to the relative contributions of F_{size_obs}. High Chl a_{total_obs} was characterized by communities having a predominant contribution of F_{micro_obs} and correspondingly lower contributions of both F_{pico_obs} and F_{nano_obs}. A similar but opposite pattern was found in η_obs, with small η_obs values clearly associated with largeF_{micro_obs}. This opposite pattern resulted from the fact that small η_obs values represent significant contributions of F_{micro_obs} essentially associated with high Chl a_{total_obs}. In addition, F_{micro_obs} and F_{pico_obs} ranged between 0.01–0.94 and 0.00–0.80, respectively, suggesting that our dataset covered a wide range of PSCs in the Pacific Arctic. According to Eq. (4), the smallest η_obs corresponded to 0.9, 0.09, and 0.01 of F_{size_MDL} for micro-, nano-, and pico-phytoplankton, whereas the largest η_obs corresponded to 0.02, 0.26, and 0.73, respectively.

Table 3 Summary statistics of primary variables used in this study. Note that these variables were determined by in situ observations. Abbreviation: Chl a, chlorophyll a; η, exponent of Chl a size distribution (CSD); F_size, fractional contribution of micro-, nano-, and pico-plankton; a_ph(443) phytoplankton absorption coefficient at 443 nm; a_NAP(443), absorption coefficient of non-algal particles (NAP) at 443 nm; a_CDOM(443), absorption coefficient of colored dissolved organic matter (CDOM) at 443 nm; and R_rs(443), remote sensing reflectance at 443 nm.

Figure 2 Ternary diagrams depicting phytoplankton size composition. Each diagram illustrates fractional contribution of micro- (F_{micro_obs}), nano- (F_{micro_obs}), and picophytoplankton (F_{micro_obs}) to total phytoplankton biomass, colored with (a) total Chl a (Chl a_{total_obs}) and (b) η_obs, respectively.

Figure 3 illustrates the biomass and fractional contribution to total Chl a of phytoplankton taxa determined by phytoclass, with respect to η_obs. The pigment ratios used in this study are detailed in Table S1 in the Supplement. Diatoms dominated in terms of both biomass and fractional contribution for small η_obs values and gradually decreased as the η_obs value increased (p<0.01). A similar but opposite pattern was observed for prymnesiophytes, indicating a gradual increase in the fractional contribution with increasing η_obs values (p<0.01). Interestingly, diatoms and prymnesiophytes were the only taxa that dominated the phytoplankton communities, while other taxa remained only minor contributors across the η_obs range. More specifically, prasinophytes and cryptophytes showed slight increases in their fractional contribution up to >0.30 at η_obs values ranging from 0.70–2.00, while their Chl a biomass in all cases remained less than 0.20 mg m⁻³. Other taxa showed negligible variations in biomass, whereas their fractional contributions fluctuated in response to reduced Chl a for the entire phytoplankton community but was statistically insignificant (p≥0.01).

Figure 3 Variations in major phytoplankton groups with reference to CSD slope. (a) Biomass and (b) fractional contribution of each phytoplankton taxa to total phytoplankton biomass (Chl a) determined by phytoclass, with respect to η_obs value. Plots and vertical bars denote the average and standard deviations of each value within the respective η_obs bins. Abbreviations: Chryso, chrysophytes; Dino, dinoflagellates; Prym, prymnesiophytes; Chloro, chlorophytes; Pras, prasinophytes; Crypto, cryptophytes; Cyano, cyanobacteria.

3.2 Phytoplankton absorption and remote sensing reflectance spectra

Since the Pacific Arctic is characterized as optically complex, i.e., the contributions of different water constituents (phytoplankton, NAP, and CDOM) are highly variable, the fractional contribution of each constituent to the total absorption by seawater (a_{total_obs}(λ)) was investigated using in situ data (Fig. S4). The ratio of a_{ph_obs}(λ) to a_{total_obs}(λ) was typically <0.30, even at wavelengths of maximum pigment absorption (i.e., 443, 469, and 488 nm) and weak pure water absorption (a_w(λ)), whereas a_{CDOM_obs}(λ) comprised 0.66±0.15 (mean ± SD) of a_{total_obs}(412). At longer wavelengths (i.e., 645, 667, and 678 nm), a_w(λ) contributed significantly to total absorption, with average values of >0.95. Overall, phytoplankton was the dominant constituent to a_{total_obs}(443) for only 30 of the 177 samples, suggesting that estimations of a_ph(λ) from R_rs(λ) using the QAA algorithm are likely to have large uncertainties for the majority of samples due to the significant contributions to absorption by other water constituents.

Figure 4 shows spectral variations in R_{rs_obs}(λ), a_{ph_obs}(λ), ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ at ten MODIS-A bands, with respect to η_obs. Larger spectral variations in R_{rs_obs}(λ), with a distinct peak at green wavelengths (i.e., 531, 547, and 555 nm), were found for smaller η_obs values, whereas larger η_obs values corresponded to relatively flat spectral shapes, with only small peaks at shorter wavelengths (i.e., 469 and 488 nm). a_{ph_obs}(λ) also showed similar differences in spectral shape and magnitude with η_obs values, except with peaks at blue wavelengths. In contrast, ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ emphasize only spectral shape by normalizing the range of variability in R_{rs_obs}(λ) and a_{ph_obs}(λ) (Fig. 6c, d). Regarding ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, sharper peaks at blue wavelengths (i.e., 412, 443, and 469 nm) with the maximum value at 443 nm were observed for large η_obs. Moreover, ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ increased more prominently with increasing wavelength from its minimum near 550 nm at smaller η_obs, whereas larger η_obs corresponded to less pronounced increases in ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ over this spectral range. Overall, the spectral features of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ exhibited clear variations associated with η_obs values, with ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ exhibiting larger variations associated with η_obs across the wide range of wavelengths compared to ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ also exhibited larger spectral variations, but differences associated with η_obs were smaller in magnitude. The performance of the modified QAA for MODIS-A bands, determined by comparing a_{ph_QAAobs}(λ) with a_{ph_obs}(λ), is shown in Table S2. According to the validation results, a_ph(λ) values at longer wavelengths (645, 667, and 678 nm) exhibited poor QAA estimation accuracy and were removed from the model development based on PCA and ML approaches. It is noteworthy that the MAE for these wavelengths represents between 25 % and 30 % of the pure water values (Pope and Fry, 1997). While this might appear large in an absolute sense, the red portion of the spectrum contains limited phytoplankton taxonomic information outside of the chlorophyll absorption band at 678 nm (Huot et al., 2005).

Figure 4 Spectral variations in key optical properties. Spectral variations in (a) remote sensing reflectance (R_{rs_obs}(λ)), (b) a_{ph_obs}(λ), (c) normalized R_{rs_obs}(λ) (${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right))$, and (d) normalized a_{ph_obs}(λ) (${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right))$ with respect to η_obs. Vertical bars represent the standard deviations at each wavelength for each η_obs range.

3.3 CSD model development: PCA approach

The spectral features of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ captured by PCA were used to develop the CSD model. Variations in the loading factors, which describe how much each variable contributes to a particular principal component at ten wavelengths (i.e., 412, 443, 469, 488, 531, 547, 555, 645, 667, and 678 nm) and seven MODIS-A bands (i.e., 412, 443, 469, 488, 531, 547, and 555 nm) for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, respectively, are shown in Figs. S2 and S3.

The spectral features captured by PCA demonstrate optical signatures of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. The regression coefficients β₀ and β_i of the logistic-type function (Eqs. 8 and 9) were therefore determined by least squares regression between the first four PC scores of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ or ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and η_obs. The resulting regression coefficients were then used to compute the model parameter C_j (Eq. 8). Here, PCA and subsequent procedures for β_i and C_j retrievals were conducted separately for two sample groups exhibiting either a_ph(412)≥a_ph(469) or a_ph(412)<a_ph(469) regarding ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, whereas the procedures for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ were performed on the entire dataset (unpartitioned) for model training. The partitioning of the model parameters for ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ was based on the trial-and-error approach (Waga et al., 2017) because a single combination of regression coefficients cannot capture the entire variations in the spectral shape of ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ in response to changing η_obs. The partitioning sequence aimed to avoid underestimation that was observed for higher η_obs (Waga et al., 2017). Since no specific pattern in η_obs estimation was identified for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, this study did not exploit the portioning approach for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. The resulting model parameters are summarized in Table S3. The resulting PCA-based CSD models for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ were hereafter denoted as $\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ and $\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$, respectively.

3.4 CSD model development: supervised ML approach

Additional CSD models were developed using a supervised ML approach through MATLAB's Regression Learner App, setting ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ or ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ as input and η_obs as output. Performance statistics for the top five and bottom five models are presented in Table 4. Comprehensive results for the 28 models appear in Tables S4 (${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right))$ and S5 (${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right))$. The best model for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ was a linear regression with linear preset, whereas that for ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ was a support vector machine (SVM) with medium Gaussian preset. These models achieved the best performance on the majority of four statistical metrics (i.e., RMSE, MSE, r², and MAE) relative to the other candidates and were thus selected as the ML-based CSD models for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$; hereafter, $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ and $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$, respectively. The model parameters for $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ is reported in Table S6.

Upon statistical evaluation, we found random patterns in relationships between model performance and regression methods. For example, the linear regression with linear interaction preset showed the second worst performance while the standard linear preset showed the best performance among all 28 models tested with ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ as input. The SVM showed the best (medium Gaussian preset) and worst (cubic preset) performance for ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. The models trained with the neural network method tended to show poor estimation accuracy for both ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. Overall, the performance of the CSD models developed by the supervised ML approach varied largely among the regression methods used in the training process, indicating that care should be taken when choosing a regression method for model development.

Table 4 Training results of the top five and bottom five CSD models based on diverse machine learning approaches (i.e., model type and preset). The four statistical metrics, including the root mean square error (RMSE), mean squared error (MSE), coefficient of determination (r²), and mean absolute error (MAE), are given as mean ± SD derived from ten repeats of five-fold cross-validation.

3.5 CSD model validation

Validation results of the four CSD models, i.e., $\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$, $\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$, $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$, and $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$ are shown in Fig. 5, with respect to the fractional contribution of a_{ph_obs}(443) to a_{total_obs}(443). The ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$-based models performed relatively well for both PCA and ML approaches, whereas, the PCA-based ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ model underestimated η_obs, with the range of estimated values (∼0.4–1.3) much lower than the measured range (∼0.2–2.2). In addition, the ML-based models ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ and $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$) showed better performance compared to the PCA-based models ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ and $\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$). Overall, the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$ performed the best among the four CSD models developed in this study. However, satellite retrieval of a_ph(λ) in optically complex waters amplifies uncertainty in retrieving η_MDLsat for the CSD models exploiting a_ph(λ). Validation results of the CSD model using the a_{ph_QAAobs}(λ), estimated from R_{rs_obs}(λ) through the modified QAA, showed diminished performance, especially for the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$ (Fig. 5f). In this sense, the best-performing model for applications with R_{rs_sat}(λ) is $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$. The $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ yielded statistical measures of 0.21 and 1.16 for MAE and bias, respectively. Out of the 53 samples in the validation dataset, estimates for 44 samples (i.e., 83 %) were within ±35 % of the in situ measured values. The associated average and median percent errors with respect to in situ values were 28.0 % and 16.2 %, respectively.

Figure 5 Validation results of the developed models. Comparison between measured (η_obs) and model-estimated η values (η_MDLobs) with respect to the fractional contribution of a_{ph_obs}(443) to a_{total_obs}(443). Upper (a–c) and lower panels (d–f) show CSD models developed by PCA- and ML-based approaches, respectively. Panels (c) and (f) show the results of the same CSD models in panels (b) and (e) but using ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{QAAobs}}\left(\mathit{\lambda }\right)$, whereas panels (b) and (c) use ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ determined from in situ observations. MAE denotes the median absolute error.

3.6 CSD slope distribution in the Pacific Arctic

Seasonal variations in climatological η_MDLsat distribution derived by the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ from R_{rs_sat}(λ) in the Pacific Arctic are shown in Fig. 6. The η_MDLsat values were persistently low in the western side of the Bering Strait, whereas those on the eastern side were generally high throughout the season. Such west-east contrast was also found on the Bering Sea shelf, with low η_MDLsat values in the west and high η_MDLsat values in the east. These spatial dynamics in the η_MDLsat would likely reflect current patterns in the Pacific Arctic. Indeed, SST shows coincident patterns with such spatial variations in η_MDLsat values (Fig. S5), with relatively higher water temperatures tending to contain higher η_MDLsat as well. The climatological mean η_MDLsat in the Pacific Arctic decreased from 1.88 to 1.52 from July to September (Fig. 7), suggesting an overall shift from smaller to larger phytoplankton communities over the season. More specifically, the fractional contribution of micro-phytoplankton (pico-phytoplankton) to total phytoplankton biomass changed from 0.04 to 0.13 (0.61 to 0.44) between July and September.

Figure 6 Monthly climatology of η_MDLsat values in (a) June, (b) July, (c) August, and (d) September in the Pacific Arctic for 2002–2022 (derived from R_{rs_sat}(λ) using the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$). White areas indicate no valid retrievals due to cloud and/or sea-ice cover.

Figure 7 Histograms of monthly climatology of η_MDLsat values in (a) June, (b) July, (c) August, and (d) September in the Pacific Arctic for 2002–2022.

4 Discussion

4.1 Taxonomic composition and size structure of phytoplankton community

Numerous studies have reported that the size structure of phytoplankton communities has strong linkages with the taxonomic composition (Finkel et al., 2010). Diatoms and dinoflagellates are generally classified as micro-phytoplankton; prymnesiophytes, chrysophytes, chlorophytes, and cryptophytes are classified as nano-phytoplankton; and prasinophytes and cyanobacteria are grouped into pico-phytoplankton. According to pigment-based taxonomic identification, diatoms and prymnesiophytes were the main phytoplankton taxa contributing to variations in the size structure of the phytoplankton communities (Fig. 3b). More specifically, a higher fractional contribution of diatoms was associated with smaller η_obs values, suggesting a large-sized phytoplankton-dominated condition. In contrast, a higher fractional contribution of prymnesiophytes resulted in larger η_obs values, indicating a small-sized phytoplankton dominated condition. Overall, shifts in the relative fractions of micro- and nano-size classes drove the change in η_obs, while pico-plankton had less impact.

4.2 Responses of optical signatures to phytoplankton size structure

The absolute concentration of phytoplankton pigments in seawater typically affects first-order variability in the magnitude of R_rs(λ), with secondary impacts on R_rs(λ) spectral shape associated with diversity in dissolved and particulate properties, such as phytoplankton community composition (Ciotti et al., 2002). Therefore, spectral variations in the magnitude-normalized ${\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)$ can be reasonably assumed to coincide with changes in the size structure of the phytoplankton community. Indeed, the spectral shape of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ showed a transition of the peak wavelength from green to blue with increasing η_obs values (Fig. 4b). Likewise, the magnitude of a_ph(λ) is related to pigment composition and concentration, whereas size information is contained in the shape of the absorption spectrum due to pigment packaging within cells (Bricaud and Morel, 1986b). For example, we found a sharp absorption peak in ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ around 443 nm that appeared to be positively correlated with CSD slope (Fig. 4d). Overall, our study demonstrated strong influences of the size structure of phytoplankton communities on ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, as reported in previous studies (Mouw et al., 2017). Although we found clear linkages in the spectral shape of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ with η_obs, it is important to note that ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ is influenced not solely by phytoplankton but also by CDOM and NAP. Since Chl a is generally uncorrelated with CDOM and NAP in coastal waters, the combined impact of absorption and scattering by all water constituents on water-leaving radiance likely accounts for the somewhat poorer performance of the remote sensing-based models compared to the in situ pigment absorption based ML models.

4.3 Comparison of PCA- and ML-based approaches

The output from PCA consisted of two terms: loading factors and PC scores. Loading factors define the rotations of the axes. PC scores are linearly uncorrelated variables that represent the positions of samples in the new rotated axes, and each is the linear combination of original spectra with corresponding loading factors (Wang et al., 2015). The PCA-based approach adopted here assumes that PC scores are correlated with η values, yet this assumption would not have been necessarily valid in this study. Indeed, the PCA-based CSD model showed a degraded performance compared to that of the ML-based model particularly for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, suggesting that the PCA could have added uncertainties in the retrieval of η. In fact, simple and direct linear regression resulted in better performance of the CSD model in the case of utilizing ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. In addition, the first two PC modes explained about 95 % of spectral variations in ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. This fact suggests that the other two PC modes (i.e., PC modes 3 and 4) contribute little to explaining the entire spectral variation but may have added uncertainties, especially considering the relatively small dataset used in the current study. Note that the PCA-based approach is a dimensionality reduction method often used to reduce the dimensionality of large data sets by transforming a large set of variables into a smaller one that still contains most of the information in the large set (Corte-Real, 2020). In the case of hyperspectral data, the input variables can easily be hundreds of wavelengths, which imposes a significant computational cost. The PCA can aggregate important spectral features into PC scores and may prove beneficial for developing robust remote sensing algorithms based on hyperspectral data. However, because we are using multispectral data with limited number of predictor variables, this potential benefit of PCA is not realized in our study.

While the conventional least square regression has been used for decades in the development of satellite ocean color algorithms (Fujiwara et al., 2011; Hirata et al., 2011; O'Reilly et al., 1998; Waga et al., 2017, 2019a), more complex ML methods are increasingly being applied and many studies have reported their capability for improved ocean color product retrievals (Chen et al., 2019; Hu et al., 2021, 2018). The least square regression is a statistical method that fits a pre-defined equation to specific data. Due to is relative simplicity, it cannot fully extract hidden patterns in data and/or elicit a deep characterization of intricate relationships between a number of interdependent variables (Martens and Dardenne, 1998). However, the ML approach of learning relationships between the input values and the corresponding output values without predefined or explicated equations requires an extensive dataset that covers complex behaviors in the data and a wide range of environmental conditions (Marzban, 2009). Once trained, ML approaches are powerful tools for the fast and efficient processing of large datasets, such as geospatial satellite data (Paul and Huntemann, 2021; Waga et al., 2022).

One of the key findings of this study is that more complex ML approaches (e.g., support vector machine, ensemble, and neural network) do not always produce more effective models than simple ML approaches (e.g., standard linear regression) (Table 4). While more complex models generally perform better than simpler ones (Makridakis et al., 2022), a complicated or flexible model will pose challenges for interpretation and can end up overfitting random effects (i.e., noises) that are unique to the dataset used for training. If these random effects are not present in new data to which the model is applied, then the model can produce incorrect results when it uses relationships developed based on random phenomena in the training dataset. Thus, the limited size of our dataset (i.e., only 177 samples) likely contributed to the poor performance of the complex ML models. Nonetheless, the CSD model trained with a support vector machine was selected as the best model for ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. This indicates that the poor performance of complex ML approaches for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ may also be associated with other regression-related factors (e.g., number of features, classifier hyper-parameter optimization, and number of cross-validation folds) rather than simply the number of samples used for training (Vabalas et al., 2019). One potential explanation for the better performance with the simple linear regression approach for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ is that variance in ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ for each η_obs range was larger compared to that of ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ (Fig. 4c). Complex ML approaches applied to ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ likely introduced errors related to the variance in the relationship between the spectral features and η_obs, whereas a simple ML approach captured only predominant features with lesser effects of the variance. Finally, we wish to also express that the type of batch approach employed by MATLAB's Machine Learning App is useful for identifying what type of model might perform well for the problem at hand, however it should not be taken as canon as more complex ML approaches often require careful customization and model design.

4.4 Methodological uncertainties and limitations

A major challenge of the ML approach, with some exceptions, such as linear regression, is that it is difficult or impossible to derive a mechanistic understanding of the model-predicted relationship between the input and output values (Ray, 2019). For this reason, the ML approaches are sometimes called “black boxes.” This lack of transparency can be problematic in interpreting the results generated by the model (Vollmer et al., 2020; Wachter et al., 2017). While ML approaches have been employed in numerous fields besides satellite remote sensing, they have not adequately addressed the issue of causality, which is essential to support wider dissemination and acceptance of the proposed models (Hall et al., 2022). What can be said at this point is that the selection of an ML approach carries with it trade-offs between accuracy and interpretability. Establishing procedures for interpreting how ML models learn and arrive at answers is crucial to not only selecting the appropriate model approach but also for improving reliability and building confidence in the selected approach.

The superior in-situ performance of a_ph(λ)-based models reflects a stronger physical coupling between η and a_ph(λ) (Fig. 5). Our preference for the R_rs(λ)-based model is operational, as it avoids uncertainties due to the inversion of a_ph(λ) from R_rs(λ) in optically complex waters and yields reliable retrievals for satellite applications; it should not be taken as evidence that η is more fundamentally linked to R_rs(λ) than to a_ph(λ). A further explanation for why the a_ph(λ)-based model performed better than the R_rs(λ)-based model pertains to measurement uncertainty related to the temporal and spatial scales of the input observations. Field data for the a_ph(λ)-based model, including pigments and absorption, were derived from analyses of well-mixed water drawn from relatively small sample volumes of few liters, resulting in high confidence that type and concentration of material analyzed for absorption was similar to the material extracted for pigments. By comparison, in situ measurements of radiometry for the computation of R_rs(λ) were measured away from the ship to avoid effects on the light field, at times that were often offset from water sampling by tens of minutes, and represented signals integrated across thousands of liters of near-surface ocean water. Therefore, uncertainty regarding sample similarity was far greater for R_{rs_obs}(λ) than for a_{ph_obs}(λ).

Our outcome metric was η, computed from within-sample size fractions rather than absolute Chl a. Prior work (Waga et al., 2017) showed that η is insensitive to reasonable choices of pore-size boundaries: percent differences in the resulting η_obs were under 5 % across three different typical Chl a_size cutoffs (i.e., >20, 2–20, and < 2 µm; >10, 2–10, and <2 µm; and >20, 5–20/<5 µm). Nevertheless, we acknowledge a small residual uncertainty for cruises that used different filters, which could add noise in heterogeneous conditions. To assess any such effect, we conducted a sensitivity check that removes cruises with differing pore-size splits (i.e., 2007, 2009, 2010) and compared model ranking and error metrics on the reduced subset. These results are summarized in Table S7, which suggests consistent findings across the entire dataset (Table 4).

Absolute Chl a can differ across analytical methods (Wang et al., 2025), yet our modeling targets a dimensionless outcome (i.e., η) computed from within-sample size fractions rather than absolute concentrations. This proportion-based normalization places fluorometer and HPLC observations on a common scale and helps mitigate method-specific bias in total Chl a. The HPLC-based Chl a_size subset in our compilation is small, which limits our ability to estimate a stable cross-method offset in η or to perform a rigorous calibration. Looking ahead, a targeted cross-calibration, paired fluorometer- and HPLC-based Chl a_size measurements collected contemporaneously across key water masses, would better quantify any residual method dependence in the retrieval of Chl a_size and further strengthen future assessments.

Overall, our dataset is heterogeneous in time, space, and methods, which introduces non-exchangeability among samples and elevates the risk of biased validation. We used a standard repeated five-fold cross-validation and an external 30 % subset to validate the performance of the developed models, but these procedures do not fully control for grouping by cruise, pore-size scheme, analytical approach, or region. As a result, cross-validated skill may be optimistic if folds inadvertently mix samples that are more similar to each other than to the broader population, and the external split may still reflect historical or regional structure (Stock, 2022; Stock and Subramaniam, 2022). Our purpose here is model ranking rather than precise absolute skill; nevertheless, the uncertainty associated with non-stratified resampling should be borne in mind when interpreting differences among approaches. A more conservative assessment is to partition the data into discrete “blocks” according to certain criteria, which enables the creation of independent training and validation folds using stratified blocking (e.g., temporal and spatial blocks) (Zhang et al., 2023). Such cross-validation strategies are preferable for heterogeneous datasets and are recommended for future work and community benchmarks.

4.5 Performance of CSD model in optically complex Pacific Arctic waters

Considering the estimation error associated with the semi-analytical IOP inversion algorithm (i.e., the modified QAA), the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$ contains large uncertainties in the retrieval of η (Fig. 5). This is primarily because the poor performance of the modified QAA in optically complex waters hampered the a_ph(λ) retrieval (Table S2), and estimation errors were propagated to the ${\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)$-based CSD model for application to satellite data. In other words, the performance of the ${\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)$-based CSD model could be improved if a more accurate IOP inversion algorithm were to be established for optically complex waters. Moreover, hyperspectral satellite sensors, such as the NASA Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission's primary sensor, the Ocean Color Instrument (OCI), and the planned Surface Biology and Geology (SBG) and Geostationary Littoral Imaging Radiometer (GLIMR) sensors will have the capability to capture more detailed spectral features of a_ph(λ) (Dierssen et al., 2023; Werdell et al., 2018), which will greatly benefit satellite-based monitoring of phytoplankton communities (Isada et al., 2015).

Considering that the accuracy goal for satellite-derived Chl a is defined as within ±35 % of the true value (Hooker and McClain, 2000), and a variety of ocean color products, such as primary productivity (Behrenfeld and Falkowski, 1997), utilize Chl a as one of the input parameters, we conclude that the CSD model developed in this study performs sufficiently well in the Pacific Arctic, presuming adequate correction for atmospheric effects in the satellite data. Since this region receives a large amount of freshwater containing CDOM and NAP delivered from rivers (Matsuoka et al., 2007), it was expected that the performance of the CSD model relying on ${\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)$ would be influenced by CDOM and NAP, which often dominate the optical properties of seawaters in this region (Chaves et al., 2015; Mustapha et al., 2012; Wang and Cota, 2003). However, the validation results suggest that the $\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$ performed with consistent accuracy regardless of the fractional contribution of a_{ph_pbs}(λ) to a_{total_obs}(λ) at 443 nm (Fig. 5).

4.6 Distribution of CSD slope in the Pacific Arctic

The Pacific Arctic, with a large continental shelf extending from the northern Bering Sea to the southern Chukchi Sea and northwards, has been characterized by a tight pelagic-benthic coupling (Grebmeier et al., 1988, 1989; Grebmeier and McRoy, 1989), with up to 70 % of primary production ultimately reaching the seafloor (Walsh et al., 1989). The seasonal cycle of sea-ice formation and melting provides suitable conditions for phytoplankton growth (Stabeno et al., 2010), with large spring diatom blooms occurring at the marginal ice edge and under the ice (Laney and Sosik, 2014; Waga et al., 2021a). The northern Bering and Chukchi Seas are reported to have the highest sinking particulate organic carbon fluxes (0.8–2.5 g C m⁻² d⁻¹) within the world ocean, and the particles collected by moored sediment traps consist of aggregates composed of diatoms exclusively (O'Daly et al., 2020). On the continental shelves in the Pacific Arctic, much of the organic carbon produced in the euphotic layer is directly transported to the seafloor with little or no grazing by zooplankton (Campbell et al., 2009). This strong pelagic-benthic coupling has maintained areas of persistently high benthic biomass, also called benthic hotspots (Grebmeier et al., 2015a), which serve as important foraging areas for upper trophic level benthivores, such as bearded seals, walrus, gray whales, and diving seabirds (Grebmeier, 2006). These hotspots are supported by influxes of organic carbon introduced by vertical transport from the overlying water column and lateral advection (Grebmeier et al., 2015b). Regarding the vertical transport of organic carbon, Waga et al. (2019a) reported that the size structure of phytoplankton communities has a significant relationship with Chl a concentration in the underlying seafloor sediments, suggesting a connection between phytoplankton cell size and benthic macrofaunal biomass in this region.

We found clear spatial variation in the distribution of η_MDLsat in the Pacific Arctic (Fig. 6). For example, on the Bering Sea shelf, the Siberian coast exhibited smaller η_MDLsat values, whereas larger values were found along the Alaskan coast. Throughout the seasons, there were west-east gradients showing smaller and larger η_MDLsat values on the Siberian and Alaskan sides of the Bering Strait, respectively. Since a small CSD slope represents a greater proportion of larger-sized phytoplankton, this result indicates larger-sized phytoplankton typically dominated along the Siberian coast, and smaller-sized phytoplankton dominated along the Alaskan coast. In the Pacific Arctic, three major water masses prevail: i.e., the Alaskan Coastal Water, Anadyr Water, and Bering Shelf Water (Coachman et al., 1976; Danielson et al., 2017). The Alaskan Coastal Water is identified with relatively high temperatures and low salinity due to freshwater input flows along the western coast of Alaska out to the Beaufort Sea (Coachman et al., 1976). The Anadyr Water, which flows along the eastern coast of Siberia, has low temperatures and high salinity, and supplies large amounts of nutrients to the Bering Sea and Bering Strait (Coachman et al., 1976). The Bering Shelf Water flows between Anadyr Water and Alaskan Coastal Water on the Bering Sea shelf and forms as these two water masses mix as they pass through the Bering Strait (Grebmeier et al., 1988). In addition to these general current patterns, satellite images of SST (Fig. S5) show distinct signatures of cold-water outcroppings in the western side of the Bering Strait, particularly in July and August. Such signatures were associated with friction between the current and the sea floor (Kawaguchi et al., 2020) and accompanied by upward nutrient flux to the surface from the nutrient-rich bottom layer of Anadyr Water (Nishioka et al., 2021), resulting in smaller η_MDLsat around the Bering Strait. These water mass distributions matched the spatial pattern in the η_MDLsat in the Pacific Arctic, suggesting a tight relationship between nutrient availability and phytoplankton cell size (Ko et al., 2020; Suzuki et al., 2021).

The η_MDLsat values in the Pacific Arctic showed clear seasonal changes from June to September (Fig. 7). According to previous studies in this region (Waga et al., 2021a; Waga and Hirawake, 2020), ice-associated spring blooms mature primarily within 20 d after sea-ice retreat and then decay gradually until fall blooms occur. Although the timing and presence/absence of spring and fall blooms largely depend on sea-ice conditions and other factors such as wind forcing (Fujiwara et al., 2018; Nishino et al., 2015), June and July are generally characterized as the post-bloom period and August and September are the typical fall bloom period. Such onset and decay of phytoplankton blooms are strongly linked to the size composition of phytoplankton communities in the Pacific Arctic (Waga and Hirawake, 2020), as shown in seasonal variations in η_MDLsat values.

5 Conclusions

This study developed a CSD model in optically complex Pacific Arctic waters by employing machine learning methods, which exploit hidden, complex relationships between optical signatures and phytoplankton size composition. Considering the large uncertainties in the inversion of a_ph(λ) from satellite-derived R_rs(λ), we used R_rs(λ) directly as a model input instead of a_ph(λ), though a_ph(λ) is more directly related to the size composition of phytoplankton communities. Neglecting the estimation errors produced from IOP inversion and considering only remotely sensed radiances and phytoplankton absorption spectra from water samples, the best-performing model among the four CSD models examined in this study was the ML-based model with normalized a_ph(λ) spectra used as input ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{SVM}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$), followed by the ML-based model with R_rs(λ) ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{LR}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$), the PCA-based model with a_ph(λ) ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{a}}_{\mathrm{ph}}\left(\mathit{\lambda }\right)}$), and finally the PCA-based model with R_rs(λ) ($\mathrm{CSD}{\mathrm{model}}_{\mathrm{PCA}-{\widehat{R}}_{\mathrm{rs}}\left(\mathit{\lambda }\right)}$). Within our dataset, the PCA-based CSD model showed a degraded performance compared to that of the ML-based model for both ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. Although the PCA-based approach assumes that PC scores are correlated with η values, this assumption would not have been necessarily valid, particularly for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. In addition, this study utilized the first four PC modes as representative for spectral features of ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$. The first two PC modes explained about 95 % of spectral variations in ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$ and ${\widehat{a}}_{\mathrm{ph}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, whereas the latter two modes contributed little to explaining the entire spectral variation but may have added uncertainties associated with the PCA step. Another key finding is that more complex ML approaches do not always produce more effective models than standard linear regression. Indeed, simple linear regression outperformed other ML approaches for ${\widehat{R}}_{\mathrm{rs}\mathrm{\_}\mathrm{obs}}\left(\mathit{\lambda }\right)$, whereas the CSD model developed with support vector machine was selected as the best for a_ph(λ). Overall, we found benefits in using ML tools to modify and improve the retrieval accuracy of the previously developed CSD model in the Pacific Arctic. Future innovations in machine learning, satellite (and airborne) ocean color sensor capabilities, and IOP algorithms can further contribute to robust, synoptic remote sensing monitoring of phytoplankton size structure in optically complex waters, such as the Arctic Ocean, where rapid change is altering the dynamics of phytoplankton with cascading effects on higher trophic levels, ecosystem functioning, and marine resources.