Influence of timing of sea ice retreat on phytoplankton size during marginal ice zone bloom period on the Chukchi and Bering shelves

The size structure and biomass of a phytoplankton community during the spring bloom period can affect the energy use of higher-trophic-level organisms through the predator–prey body size relationships. The timing of the sea ice retreat (TSR) also plays a crucial role in the seasonally ice-covered marine ecosystem, because it is tightly coupled with the timing of the spring bloom. Thus, it is important to monitor the temporal and spatial distributions of a phytoplankton community size structure. Prior to this study, an ocean colour algorithm was developed to derive phytoplankton size index FL, which is defined as the ratio of chlorophyll a (chl a) derived from cells larger than 5 μm to the total chl a, using satellite remote sensing for the Chukchi and Bering shelves. Using this method, we analysed the pixel-bypixel relationships between FL during the marginal ice zone (MIZ) bloom period and TSR over the period of 1998–2013. The influences of the TSR on the sea surface temperature (SST) and changes in ocean heat content (1OHC) during the MIZ bloom period were also investigated. A significant negative relationship between FL and the TSR was widely found in the shelf region during the MIZ bloom season. However, we found a significant positive (negative) relationship between the SST (1OHC) and TSR. Specifically, an earlier sea ice retreat was associated with the dominance of larger phytoplankton during a colder and weakly stratified MIZ bloom season, suggesting that the duration of the nitrate supply, which is important for the growth of large-sized phytoplankton in this region (i.e. diatoms), can change according to the TSR. In addition, under-ice phytoplankton blooms are likely to occur in years with late ice retreat, because sufficient light for phytoplankton growth can pass through the ice and penetrate into the water columns as a result of an increase in solar radiation toward the summer solstice. Moreover, we found that both the length of the ice-free season and the annual median of FL positively correlated with the annual net primary production (APP). Thus, both the phytoplankton community composition and growing season are important for the APP in the study area. Our findings showed a quantitative relationship between the interannual variability of FL, the TSR, and the APP, which suggested that satellite remote sensing of the phytoplankton community size structure is suitable to document the impact of a recent rapid sea ice loss on the ecosystem of the study region.


Introduction
The wide continental shelf from the northern Bering Sea to the southern Chukchi Sea is known as one of the most biologically productive regions in the world ocean (Springer and McRoy, 1993).Nutrient-rich water forms during wintertime as a consequence of strong vertical mixing, which supports an extremely high primary productivity from spring to autumn (Springer and McRoy, 1993).The high primary production in the region is not completely consumed by the grazers in the water column because of the low grazing pressure, fast sinking rate of large chain-forming diatoms, and shallow bottom depths (Grebmeier et al., 1988(Grebmeier et al., , 2006a)).As a result, this region supports a large amount of benthic biomass.Subsequently, top predators such as walrus (Odobenus rosmarus), grey whale (Eschrichtius robustus), and various sea birds also gather to feed on the rich benthic ecosystem (Sheffield et al., 2001;Moore et al., 2003;Feder et al., 2005).Thus, the marine ecosystems in this region have been characterized by rather simple food webs and efficient energy transport to higher trophic levels through tight pelagic-benthic coupling (Grebmeier and Dunton, 2000).However, recent sea ice decline and ocean warming in the northern Bering Sea have triggered a northward shift of this tightly coupled pelagicbenthic ecosystem (Grebmeier et al., 2006b).Ocean warming can modify the ecosystem structure from a pelagic-benthic to a pelagic-pelagic type that is generally found in the southern region (Piepenburg, 2005).Marine ecosystems in the northern Bering and Chukchi seas are now facing extreme climate change.Therefore, it is important to improve our comprehension of how these ecosystems respond to environmental change.Several studies have already reported changes in the composition and distribution of mammals, sea birds, and benthic species in the Chukchi Sea (Grebmeier, 2012 and references there in), although studies on the changes in the phytoplankton taxa are relatively scarce.Since even small changes in the lower trophic levels can greatly affect highertrophic-level organisms through the efficient and short energy transfer pathway of the arctic marine ecosystems (Grebmeier et al., 2010), it is essential to clarify how primary producers (i.e.phytoplankton) respond to recent environmental changes in terms of their size structure and productivity.
Satellite ocean colour remote sensing has enabled us to monitor the spatial and temporal dynamics of marine phytoplankton in diverse ecosystems.Several studies have attempted to measure not only the chlorophyll a (chl a) and primary productivity over the euphotic zone (PP eu ), but also the phytoplankton functional types (PFTs) from space (e.g.Alvain et al., 2005;Uitz et al., 2006;Kostadinov et al., 2009;Brewin et al., 2010;Mouw and Yoder, 2010;Hirata et al., 2011).The satellite remote sensing of the PFTs is a powerful way to assess biological contributions to biogeochemical cycles and ecosystem variability over large spatial scales and continuous temporal scales.For example, Takao et al. (2012) found a significant relationship between the interannual variability of the phytoplankton community composition and the primary production in the Indian sector of the Southern Ocean.Hirata et al. (2013) showed the potential of the longterm monitoring of ecological provinces using the satellitederived PFTs and the PFTs derived from a marine ecosystem model at the global scale.Thereby, the satellite remote sensing of the PFTs is expected to contribute to an understanding of the phytoplankton dynamics in the Chukchi and Bering seas.However, the Chukchi and Bering seas are known to be optically complex waters, and their optical properties are quite different from global pelagic waters.Large amounts of absorption from coloured dissolved organic matter (CDOM) and highly packaged phytoplankton absorption cause large estimation errors in ocean colour products (Matsuoka et al., 2007;Naik et al., 2013).To overcome this problem, the use of regionally optimized ocean colour algorithms instead of globally developed algorithms is required to accurately monitor the chl a, PFT, and PP eu .Cota et al. (2004) successfully modified a chl a algorithm to estimate the chl a for highly packaged, high-latitude phytoplankton from space.Regionally optimized PFT algorithms for the Chukchi and Bering seas were proposed by Fujiwara et al. (2011) to estimate the percentage of the contribution of > 5 µm algal cells to the total chl a biomass as an index for the phytoplankton community size structure (hereafter F L ).An optimal algorithm to derive PP eu was also developed by Hirawake et al. (2012) for the western Arctic seas.The application of these regional ocean colour algorithms to long-time-series satellite data is expected to contribute to our comprehension of phytoplankton community changes.
In this study, we focused on the timing of the sea ice retreat (TSR) and the phytoplankton size composition in the marginal ice zone (MIZ), which is the area where ice melt has just recently occurred.It has been shown that the TSR can alter the phytoplankton community composition in the northern Chukchi Sea in late summer (Fujiwara et al., 2014).However, little is known about the relationship between the TSR and the phytoplankton size composition in the shelf region of the Bering and Chukchi seas, and during the spring phytoplankton bloom period.Spring phytoplankton blooms generally occur in the MIZ in the Arctic Ocean (Kahru et al., 2011;Perrette et al., 2011).Because zooplankton strongly depend on the timing and magnitude of the spring bloom for their growth and production (e.g.Hunt et al., 2002Hunt et al., , 2011;;Søreide et al., 2010), the influence of the TSR on the phytoplankton size composition during the MIZ bloom period is crucial to evaluate the bottom-up effects of the primary production on the food web.In this study, we aimed to clarify the relationships between the temporal and spatial distributions of the phytoplankton size composition and the marine environment during the MIZ bloom periods using satellite remote sensing.Second, we aimed to assess how the phytoplankton size contributes to the annual net primary production (APP).

Satellite data
Satellite-derived level-2 and level-3 standard mapped images of the daily ocean colour data set were obtained from GSFC/DAAC NASA (1 km and 9 km resolutions, respec-  and 412, 443, 488, 555, and 667 nm for Aqua-MODIS), euphotic depth (Z eu ) derived from the quasi-analytical algorithm (QAA) proposed by Lee et al. (2007), and photosynthetic available radiation (PAR) were downloaded.The reprocessing versions for the products were 2010.0 and 2013.1 for SeaWiFS and Aqua-MODIS, respectively.The R rs (λ) values from SeaWiFS were converted to R rs (λ) values for MODIS using conversion factors that removed the biases between them (see Appendix A).The definitions of the symbols are listed in Table 1.The satellite-observed sea ice concentration (SIC) and sea surface temperature (SST) were used to assess the environmental conditions in the region.
The SIC derived by the Special Sensor Microwave Imager-Defense Meteorological Satellite Program (SSMI-DMSP) using NASA Team algorithm 2 (Markus and Cavalieri, 2000) was downloaded from the National Snow and Ice Data Center (NSIDC).The temporal and spatial resolutions were daily and 25 km, respectively.Following Perrette et al. (2011), nearest-neighbour interpolation was performed on the SIC data set to convert the 25 km spatial resolution into a 9 km resolution.Thus, the data were on the same scale as the ocean colour images.The AVHRR-derived SST was downloaded from the Physical Oceanography Data Active and Archive Center (PODAAC) for 1998-2012, and the Aqua-MODISderived SST was obtained from GSFC/DAAC NASA for 2003-2013 to cover the ocean colour time-series.

Calculation of heat flux
The development of thermal stratification is a major controlling factor of the nutrient supply to the upper layer from spring to summer in both the Bering and Chukchi seas (e.g.Hill et al., 2005;Coyle et al., 2008;Sambrotto et al., 2008).
To understand how thermal stratification develops in correspondence to the yearly change in the sea ice retreat timing in the region, the changes in the ocean heat content ( OHC) during the MIZ bloom season (see 2.3) were used as a proxy of the surface-mixed-layer-depth development.
OHC was calculated following the equation of Mizobata and Shimada (2012): where SHTFL, LHTFL, NSWRS, and NLWRS are the sensible heat flux, latent heat flux, net shortwave radiation, and net longwave radiation, respectively.All of the products used in Eq. (1) were obtained from the NCEP/NCAR daily reanalysis data set (Kalnay et al., 1996).These Gaussian gridded reanalysis products were regridded into polar-stereographic projections and then interpolated into the 9 km resolution using the nearest-neighbour method to achieve the same resolution as the satellite products.

Data processing
The chl a concentration was computed using the Arctic-OC4L algorithm proposed by Cota et al. (2004), which was optimized for the optical properties of phytoplankton in the Arctic Ocean: where R is the maximum band ratio defined as log(R rs ( 443 band because there is no observation band around 510 nm in MODIS.We found there was no statistically significant difference in the chl a values with and without the inclusion of R rs (510) (data not shown).
An index for the phytoplankton community size composition (F L ) was defined as the ratio of the chl a attributed to cells larger than 5 µm (chl a >5 µm ) to the total chl a (chl a total ): (3) The satellite F L value was estimated using the phytoplankton size derivation model (SDM) proposed in Fujiwara et al. (2011), which was optimized for the phytoplankton communities and optical properties in the Bering and Chukchi seas.The SDM requires the ratio of the phytoplankton absorption coefficient (a ph (λ)) and spectral slope of the particle backscattering coefficient (γ ) to compute F L .a ph (λ) was calculated using QAA version 5 (QAA v-5; Lee et al., 2009), but the spectral slope of the absorption coefficient of gelbstoff+detritus (S dg ) was modified for the region (Appendix B) to avoid the retrieval of a negative a ph (λ).γ was empirically quantified using the R rs (λ) ratio (Fujiwara et al. 2011).Then, F L could be derived using Eq. ( 4): • 100(%), ( 4) where X 0 = 3.175, X 1 = −0.570and X 2 = −0.565,and λ 1 = 488 and λ 2 = 555.The PP eu was derived using the absorption-based productivity model (ABPM) proposed in Hirawake et al. (2011Hirawake et al. ( , 2012)).Although the ABPM was originally based on a vertically generalized productivity model (VGPM;Behrenfeld and Falkowski, 1997), it retrieves optimal values of chl anormalized productivity (P B opt ) from a ph (λ) instead of from SST and chl a.The use of a ph (λ) is suitable to discuss the effect of ocean warming on PP eu , because the PP eu value derived by ABPM is independent of the temperature (Hirawake et al., 2011).The ABPM also effectively reduces the optical effect of high-CDOM water in the study region via its a ph (λ) calculation steps using QAA (Hirawake et al., 2012).The ABPM was used to determine PP eu in the following equations: and, where E 0 and DL indicate the daily surface photosynthetically available radiation and day length, respectively.DL was computed from the date and latitude following the method of Brock (1981).
Since ocean colour products can only be retrieved from cloud-and ice-free pixels in the satellite imagery, we computed the 9-day moving average for a ph (443) and Z eu to increase the retrievals of valid PP eu values.We also applied nearest-neighbour spatial interpolation to the daily satellite images.The TSR was calculated using the daily SIC data for each year.In this study, non-ice-covered pixels were defined using a SIC < 10 % following Pabi et al. (2008), and the TSR was defined as the last date when the SIC fell below 10 %, prior to the observed annual sea ice minimum across the study region during summer (50-75 • N, 170-210 • E).The length of the open-water period for each pixel was defined as the number of days where the SIC was less than 10 %.Perrette et al. (2011) showed that high chl a values occur during the MIZ bloom period, which is defined as 20 days after the sea ice melt in ∼ 90 % of the Arctic Ocean.However, Niubauer et al. (1991, 1995) reported that the sea-icerelated spring bloom lasts for ∼ 2 weeks on the Bering Sea shelf unless there is a mixing event.To clarify the relationship between the TSR and phytoplankton size composition during the bloom period, we computed the 14-day averages of the variables (F L , SST, and PAR) after sea ice melted for each pixel, and defined them as the MIZ-bloom-period values.Then, a correlation analysis (Spearman's rank correlation coefficient, ρ) was conducted to evaluate the interannual variations between the MIZ bloom time variables and TSR for every one-by-one pixel.Note that the spatial statistical analysis was conducted only for the regions where sea ice was observed in at least 11 of the 16 years.
A standardized multiple regression analysis was also performed to determine the variables that contributed to the interannual variability of the APP.The APP was computed by integrating PP eu for every non-ice-covered pixel.The annual median PP eu was calculated for each pixel, and was substituted to accurately compute the APP when the PP eu was missing as a result of cloud cover.We conducted sensitivity tests and found that the substitution of the annual median PP eu to the missing values causes only a small error (∼ 10 % in maximum) to calculate the APP (data not shown).Then, the APP was regressed using the annual median F L , SST, and length of the open-water period for every one-by-one pixel: where A 1 , A 2 , and A 3 indicate the partial regression coefficients, and OWP denotes the length of the open-water period (i.e. the number of non-ice-covered days).All of the variables were standardized before conducting the multiple regression analysis by subtracting the climatological mean and dividing by the standard deviation.Thus, the contribution of the variables in Eq. ( 7) to the APP was comparable to using partial regression coefficients.Note that the multiple regression analysis was performed for areas where seasonal ice cover had been observed every year.The MATLAB statistical toolbox (The Math Works, Inc.) was used for the statistical analyses in this study.

In situ data sampling to evaluate satellite products
The in situ total and size-fractionated chl a samples were obtained during the cruises of the Bering Arctic Subarctic Integrated Survey (BASIS) program conducted in late summer to autumn (August to October) of 2005, 2006, and 2007, and the cruises of the GRENE Arctic Climate Change Research Project conducted in late summer to autumn of 2012 (September to October) and early summer of 2013 (June to July) on the Bering and Chukchi shelves (Fig. 1).For these samples, 138-525 mL of seawater was filtered onto 10, 5, and 2 µm pore-sized nuclepore polycarbonate filters and GF/F filters.All of the chlorophyll samples were frozen onboard and stored in liquid nitrogen or in a supercold freezer (−80 • C; except during 1-2 days when the samples were shipped in coolers on ice).They were analysed later (within 1-8 months) using the acidification technique (Parsons et al., 1984) with a Turner Designs 700 fluorometer for the samples from the BASIS cruises and the nonacidification technique (Welshmeyer, 1994) with a Turner Designs 10 AU fluorometer for the samples from the GRENE cruises.While different chl a measurement methods were used for the different cruises, the fluorometers were calibrated with a pure chl a standard prior to the analyses to minimize the differences between the cruises.Then, the in situ F L was calculated using Eq. ( 3) and compared with the daily matching satellite-derived F L (MODIS L2).The root- mean-square error (RMSE) between the satellite-derived F L and in situ-measured F L was calculated to evaluate the SDM performance.The in situ R rs (λ) was also measured using a PRR800/810 free-fall profiler (Biospherical Inc.) during the GRENE cruise of 2013.The in situ R rs (λ) was compared with the daily matching MODIS-observed R rs (λ) to assess the difference between them and its influence on the accuracy of F L retrieval from space.
The in situ primary productivity (PP) data, as well as the chl a data, were also used to assess the influence of the subsurface chl a maximum (SCM) on the remotely estimated PP eu .PP samples were collected during the TR/S Oshoromaru (Hokkaido University, IPY and GRENE cruise) and R/V Mirai (JAMSTEC, GRENE cruise) cruises in the same region, along with chl a samples (Fig. 1).Cruises were conducted in early summer to autumn (June to October) of 2007, 2008, 2012, and 2013.PP was measured using the stable 13 C isotope method (Hama et al., 1983) for several optical depths, and then PP eu was calculated by integrating the PP values from the surface to Z eu (see Hirawake et al., 2012 for more details of the sampling and analysis).

Evaluation of performance of satellite algorithms
The accuracy of the SDM-derived F L was evaluated by comparing the F L values from the in situ measurements and daily matched MODIS level-2 data set.Twenty-five data points were available for this examination, which were collected over a wide area of the Bering and Chukchi seas during different seasons (Fig. 1). Figure 2a compares the satellitederived F L and in situ F L .The SDM successfully retrieved the F L values for 17 of the 25 data points (68 % of the data) within a ± 20 % F L range.The RMSE was 25 %.The satellite validation was very similar to the results of Fujiwara et al. (2011), who showed that the F L derived using the in situ-measured R rs (λ) had a 69 % accuracy and an RMSE of 22.7 %.However, it should be taken into account that there was a slight overestimation in the low F L range and relatively large underestimation in the high F L range (slope = 0.48 and intercept = 0.18).We also evaluated the performance of R rs (λ) retrieval of MODIS comparing with in situ-measured R rs (λ) at 13 data points during the GRENE  2), and (c) F L derived from MODISconverted R rs (λ) vs. F L derived from in situ R rs (λ) (r 2 = 0.85, p < 0.01, N = 220).Dashed lines are 1 : 1 lines, dotted lines indicate ± 20 % from 1 : 1 (Fig. 2a and c), and coloured solid lines represent the regression line, respectively.The regression line of Fig. 2a is also shown in Fig. 2c with the blue solid line for comparison.cruise of 2013. Figure 2b shows the comparison of in situand satellite-derived R rs (λ), and we found that the MODIS significantly underestimates R rs (λ) at every wavelength (the slopes ranged from ∼ 0.34 to ∼ 0.46).The slopes and intercepts between the two R rs (λ) are listed in Table 2, which was used as the factors to convert in situ R rs (λ) to MODIS-R rs (λ).The conversion factors were applied to the same R rs (λ) data set used in Appendix B, and then, we compared the SDM performance between in situ R rs (λ) and R rs (λ) converted to MODIS-R rs (λ) to assess how the retrieval errors in the MODIS-R rs (λ) affect estimation accuracy of the F L .The relationship between the two derived F L is shown in Fig. 2c, and we found similar relationship with Fig. 2a.It can be said that the underestimations of R rs (λ) significantly caused the underestimation of F L especially in the middle to high range of F L .Although we conducted the optimization of the QAA to derive more accurate input parameters of the SDM (Appendix B), retrieval error in the R rs (λ) significantly affected the SDM accuracy in the study region.Nevertheless, the high determination coefficient indicated that the correlation was sufficient to address the spatio-temporal variability of F L .Here, we conclude that the SDM is applicable to satellite remote sensing in the Bering and Chukchi seas, and is effectively optimized for what is known as optically complex water (Matsuoka et al. 2007, Naik et al. 2013).
To confirm how the satellite-derived F L and PP eu represent a water column's phytoplankton size structure and productivity, we compared the surface and vertically integrated F L values (calculated using Eq. ( 3) with water-column-integrated chl a total and chl a > 5µm ) using an in situ data set.The vertically integrated F L showed a significant relationship with the surface F L (Fig. 3c, r 2 = 0.67,p < 0.01) in spite of the presence of the subsurface F L maximum at the many of the sampling sites, especially in low and middle range of the surface F L were observed (Fig. 3b).Similarly, PP eu was also significantly correlated with the surface PP (Fig. 3c), and ver-tical distribution of PP also represents the large contribution of surface PP to PP eu (Fig. 3e).Although the surface values of both F L and PP explain the variation of water-columnintegrated values, the depth of the maximum PP (6.1 ± 8.9 m) was significantly shallower than the depth of the maximum chl a (20.1 ± 13.7 m; Fig. 3d).

Timing of open-water bloom
In order to determine the timing of the spring phytoplankton bloom, we investigated the relationship between the TSR and date of occurrence of the yearly chl a maximum (CMAX).Figure 4a shows the climatologic median of TSR from 1998 to 2013.It indicates that the sea ice normally retreats by the summer solstice, except in the northwestern edge of the study area.The timing of CMAX showed a spatial pattern similar to that of TSR (Fig. 4b), except for coastal Alaska such as in Kotzebue Sound and Norton Sound.In these areas, CMAX occurred in late summer.However, CMAX generally occurred immediately after the sea ice retreat (fewer than ∼ 20 days after the sea ice retreat) over almost the entire shelf area (Fig. 4c).However, the CMAX date was delayed near the Alaska coast, where the Alaska Coastal Water usually flows, and in the western Bering Strait.In such regions, the timing of CMAX compared to TSR was delayed by more than 50 days.

Influence of TSR on phytoplankton size composition during MIZ blooms
Spearman's rank correlation coefficient was calculated for every one-by-one pixel between F L and TSR, SST and TSR, OHC and TSR, and PAR and TSR (Fig. 5a-d).The statistical significance proportions of ρ are listed in Table 3.We found that during the MIZ period, F L was negatively associated with TSR in 68 % of the shelf area, with a significant correlation in 15 % of the area (p < 0.05).A significantly pos-  itive (p < 0.05) relationship was found for the western side of the Bering Strait and western coast of Alaska; this area accounted for 2 % of the whole study area.During the MIZ period, SST was tightly and positively correlated with TSR over most (92 %) of the region (Fig. 5b), with a significant correlation in 65 % of the area.The spatial pattern of OHC was also similar to the SST pattern though it was inverse, with a statistically significant negative relationship for more than 61 % of the area (Fig. 5c).These results reveal that an earlier sea ice retreat is associated with a cold surface temperature and relatively large amount of heat release from the sea surface during the MIZ period.However, on the north-  western edge of the Chukchi Sea, the opposite correlations were found, with a significantly positive ρ between SST and TSR, and significantly negative ρ between OHC and TSR (Fig. 5b and c).PAR was strongly and positively correlated with TSR in the northern Bering Sea up to and through the Bering Strait, but became negative in the northern part of the Chukchi Sea shelf (Fig. 5d).The highest positive ρ between PAR and TSR was found in the area where the sea ice normally retreats before the summer solstice, at the time when the largest daily PAR would be observed (Fig. 4a).

Factors controlling annual net primary production
A standardized multiple regression analysis was used to determine the contribution of the controlling factors for the APP.Here, we compared the partial regression coefficients for the annual median values of the length of the open-water period, F L , and SST to determine the contribution to the interannual change in the APP.The spatial distributions of the partial regression coefficients for each variable are shown in Fig. 6a-c.All of the variables had positive coefficients for most of the study region.Accordingly, the APP was also successfully modelled (p < 0.05, F test) using these variables in all of the regions except the Gulf of Anadyr and coastal Alaska (Fig. 6d).This analysis revealed that longer growing seasons, larger proportions of large-sized phytoplankton assemblages, and higher SSTs generally enhanced the APP.However, the magnitude of the contribution changed regionally.Larger partial regression coefficients for the length of the open-water period (> 0.7) were found mainly in the northern shelf area of the Chukchi Sea (Fig. 6a).Larger partial regression coefficients for F L (> 0.7) were mainly found on the Bering Sea shelf and part of the central Chukchi Sea shelf (67-70 • N, 170-175 • W), where other variables showed relatively lower contributions (Fig. 6b).Conversely, the partial regression coefficients for the SST were negative on the southern edge of the Bering Sea shelf study area (Fig. 6c).
Based on the results of the standardized multiple regression analysis, in regions with seasonal ice cover, a large amount of phytoplankton positively contributed to the APP on the northern Bering Sea shelf, and a longer open-water period was a major factor for the APP in the Chukchi Sea shelf.

Evaluation of performance of satellite algorithms
The validation of the satellite-derived F L showed a sufficient correlation with the in situ F L , although the vertical distribution of the phytoplankton should be confirmed, because the SCM is commonly distributed in the study area and high Arctic Ocean (e.g.Hill et al. 2005, Ardyna et al., 2013).The omission of the SCM sometimes causes a large error in the satellite estimation of PP eu in the high Arctic (Hill et al., 2013).However, Arrigo et al. (2011b) and Ardyna et al. (2013) showed that the omission of the SCM can lead to only a small error in the PP eu retrieval over most of the Arctic Ocean.Our results supported the latter idea that surface PP explained the variation of PP eu well (Fig. 3c).Similarly, surface F L showed a very good relationship with F L retrieved from water-column-integrated chl a in spite of the presence of subsurface F L maximum (Fig. 3a and b).A similar relationship was also reported between the surface and water-column-integrated chl a in the Bering shelf (Lomas et al., 2012).Thereby, we suggest that the surface F L can also be reasonably used to predict the upper layer phytoplankton size structure.This might be due to the shallow bathymetry and SCM depth (∼ 20 m) associated with the nitracline depth (Brown et al., 2015) compared to the higher Arctic Ocean where a deeper SCM (> 40 m) is commonly found (Ardyna et al., 2013).Furthermore, the significantly shallower depth of the PP maximum compared to the depth of the SCM (Fig. 3d) causes the PP at the SCM to make a smaller contribution to PP eu .This is also consistent with the results shown by Brown et al. (2015) for the Chukchi shelf.Hence, we believe that ocean colour remote sensing is applicable to discuss the temporal and spatial relationships between the distribution of sea ice and phytoplankton variables (i.e.F L , PP eu , and chl a), at least in the study area.

Timing of MIZ bloom
We have shown that the timing of the MIZ bloom occurrence is tightly coupled with the TSR, except for the western side of the Bering Strait and coastal Alaska.Consistent with the results of Perrette et al. (2011), CMAX generally occurs within ∼ 20 days after the sea ice melts (Fig. 4c).Some previous studies have discussed the relationship between the timing of the bloom and the sea ice melt (e.g.Hunt et al., 2002Hunt et al., , 2011;;Kahru et al., 2010;Perrette et al., 2010).Sigler et al. (2014) found that an earlier sea ice retreat (earlier than mid-March) could not trigger an ice-related phytoplankton bloom, because there was insufficient light to support the pri-A.Fujiwara et al.: Influence of timing of sea ice retreat on phytoplankton size mary production as a result of the well-mixed water column in the southern Bering Sea (Hunt et al., 2002).In these years, the spring bloom occurs in late spring (May or early June).In contrast, the bloom occurs immediately after the sea ice melts in the late-ice-retreat years when it occurs later than mid-March.Brown and Arrigo (2013) evaluated the relationship between the timing of the sea ice retreat and spring phytoplankton blooms using satellite remote sensing data, and showed that the relationship described above is applicable in the southern but not in the northern Bering Sea.Since our analyses were conducted on more northward-centred regions compared to the studies of Hunt et al. (2002Hunt et al. ( , 2011)), our results revealed that the timing of CMAX generally followed the TSR, which was consistent with the results of previous studies (Kahru et al., 2011;Perrette et al., 2011;Brown and Arrigo, 2013;Ji et al., 2013).However, our results also indicated that the CMAX date is delayed in the Bering Strait and along the coast of Alaska (Fig. 4c).This suggests that the timing of the annual maximum phytoplankton biomass does not fully depend on the sea ice retreat timing and the subsequent nutrient conditions in these areas.In the Bering Strait, the CMAX delay would be because the large nutrient supply continues even in summer as a result of the inflow of nutrient-rich Anadyr Water (Springer and McRoy, 1993;Springer et al., 1996).Further, in coastal Alaska, optical contamination of the CDOM during the chl a retrieval can also affect the CMAX values, especially in the area where the Alaskan Coastal Water (ACW) flows.Matsuoka et al. (2011) suggested that the ACW contributes to the high value of CDOM absorption (a y (λ)), which can cause an overestimation of satellite chl a values, because the CDOM absorbs light at some of the same wavelengths as chl a (Matsuoka et al., 2007;Naik et al., 2013).Accordingly, in the ACW, especially in the Kotzebue and Norton sounds, CMAX should be carefully treated, because the chl a values in this area might contain large levels of uncertainty.
Our results revealed that the spring open-water bloom generally occurs during the MIZ period, at least in our study area (hereafter referred to as the MIZ bloom).Although the MIZ bloom timing is important for the growth and/or reproduction of zooplankton and more higher trophic organisms in the seasonally ice covered sea (e.g.Hunt et al., 2002Hunt et al., , 2011;;Leu et al., 2011), it is also crucial to comprehend how the phytoplankton community size structure during the MIZ bloom varies with the yearly change in the TSR, considering the energy use through predator-prey body size relationships.

Influence of TSR on phytoplankton size composition during MIZ blooms
A spatial analysis using the continuous ocean colour data from SeaWiFS and Aqua-MODIS observations revealed that the phytoplankton community size structure during the MIZ bloom responds to the yearly change in the TSR.Our results indicated that large phytoplankton tend to increase and/or maintain their high biomass proportion in early-ice-retreat years, and vice versa (Fig. 5a).The relationships between the SST and TSR, and between OHC and TSR, suggest that the development of the surface mixed layer is related to changes in the TSR.Unfortunately, the salinity stratification is difficult to infer using satellite or reanalysis data.However, Alexander and Niebauer (1981) showed that the development of both salinity and thermal stratification induces water column stability, which generally triggers the ice-edge phytoplankton bloom in the Bering Sea.Thus, OHC was used as a proxy of the development of the surface mixed layer.Because the sea ice retreat occurs before the summer solstice, even in the late-ice-retreat years on the shelf (Fig. 4a), the thermal stratification is likely to be delayed as a result of the weak solar radiation and low air temperature in early-ice-retreat years.In contrast, stronger solar radiation and warmer air temperatures can lead to faster stratification during the MIZ bloom period in late-ice-retreat years.Hence, high surface-nutrient concentrations are expected to remain for a longer period after the sea ice melts in early-ice-retreat years compared to the late years, and conversely, the nutrients can immediately be consumed by large phytoplankton in late-ice-retreat years.
In addition, under-ice phytoplankton or ice-algal blooms can also be controlling factors for the relationship between F L and the TSR.Arrigo et al. (2012) discovered that phytoplankton blooms can occur even in under-ice water columns where sufficient light for phytoplankton growth penetrates into the water as a result of thin ice or the presence of melt ponds.Lowry et al. (2014) suggested that under-ice phytoplankton blooms could frequently and widely occur on the Chukchi Sea shelf, where a negative ρ between the TSR and F L was found in the current study.Because thin first-year ice dominates in the shelf areas of the Chukchi and Bering seas (Comiso et al., 2008), sufficient light can penetrate into the water column through a melt pond or fragile ice when the solar radiation is strong enough during the melting season (Arrigo et al., 2012).We noted that the sea ice retreats before the summer solstice (i.e.before the yearly solar radiation maximum occurs) in the focus area.Accordingly, the light levels of the under-ice water column would be better for phytoplankton growth in late-ice-retreat years, and underice phytoplankton blooms may utilize nutrients in the surface layer before the sea ice retreats.As a result, the large phytoplankton taxa during the MIZ blooms cannot retain their dominance for very long in years with a late ice retreat, because as the surface nutrients become depleted, smaller phytoplankton taxa are favoured.
However, these scenarios are not the case at the western side of the Bering Strait (2 % of the study area), where the correlation between F L and the TSR showed a positive significant relationship (Fig. 5a; Table 3).It is notable that this is the same area where the CMAX date does not follow the date of the sea ice retreat (Fig. 4c).This analysis revealed that a late sea ice retreat does not limit the phytoplankton growth under post-bloom conditions in this area.One possible scenario is that the major limiting factor for phytoplankton growth is the water temperature rather than the light and nutrient availability.The Bering Strait is known to be a nutrient-rich region even in summer (Springer and McRoy, 1993), and there is sufficient light.Another scenario is that higher nutrient consumption in upstream regions can reduce the horizontal advection of nutrients into the Bering Strait in early-ice-retreat years, and vice versa.
There are several studies that referred to the importance of the local wind field to determine nutrient conditions and chl a concentration near the ice edge instead of TSR.For example, Niebauer et al. (1995) showed the effect of the wind velocity on the initiation and duration of the ice-edge bloom in the southern Bering shelf.Arrigo et al. (2014) also reported the wind direction in the northern Chukchi shelf has large impact on the upwelling of nutrients from shelf break.Although interannual variability of such wind field can also affect the phytoplankton community size structure during the MIZ bloom period in the local scale, our spatial statistical analysis revealed a general trend of the bloom time F L toward TSR.It has been reported that the sea ice retreat timing can alter the seasonal succession of phytoplankton in Svalbard (Leu et al., 2011) and the northern Chukchi Sea (Fujiwara et al., 2014).Our results also suggest that the phytoplankton community size composition can change with the timing of the sea ice retreat in the shelf region of the Bering and Chukchi seas during the MIZ bloom.These results provide important information to help understand the efficiency of energy transfer for organisms that depend on the MIZ blooms as a food source.However, it should be noted that the relationships between the TSR and phytoplankton size might not be applicable to the entire Arctic, because our region of focus is characterized by a shallow bathymetry, sea ice retreat near the summer solstice (Fig. 4a), and the dominance of thin first-year ice.These unique environmental characteristics can tightly couple the phytoplankton bloom and sea ice dynamics.

Factors controlling annual net primary production
Previous studies have reported that the variability of the annual primary production in the Arctic Ocean largely depends on the length of the open-water period and the open-water area, which contribute to the length of the growing season and the growing area, respectively (Arrigo et al., 2008(Arrigo et al., , 2011a;;Pabi et al., 2008).A similar relationship was also reported for the Bering Sea, where is seasonally covered by sea ice (Brown et al., 2011(Brown et al., , 2012)).Further, Trembley and Gagnon (2009) showed that the cumulative irradiance during the growth season is not always a primary driver of the APP in the Arctic.The initial concentration of nitrate at the onset of the growth season is very important in the high Arctic Ocean (Trembley and Gagnon, 2009).However, the nitrate concentration uniformly reaches 10-20 µM at the onset of the growth season, even in the surface water, as a result of strong vertical mixing in both the Chukchi and Bering shelves (e.g.Sambrotto et al., 1986;Codispoti et al., 2005), which is enough to support a large open-water and/or underice bloom.Thus, our findings suggest that the length of the open-water period plays a significant role in controlling the APP in the study region, especially in the northern Chukchi shelf (Fig. 6a).Nevertheless, the contribution of the length of the open-water period to the APP is relatively smaller in the Bering Sea.In contrast, the contribution of F L to the APP is larger than that of the open-water period in the Bering and central Chukchi Sea shelves (Fig. 6b).This is because, although these shelf regions have similar ranges for the interannual variation of the open-water period (1-2 months), it accounts for up to ∼ 35 % of the mean growth season on the Chukchi shelf (∼ 6 months), which is significantly larger than that on the Bering shelf (up to ∼ 20 %).It is likely that the dominance of highly productive large phytoplankton (i.e.diatoms) greatly affects the variability of the APP in regions where the open-water period is relatively longer.Therefore, some controlling factors for the phytoplankton size structure, such as the nutrient supply or grazing, are suggested to be important to determine the APP in the southern shelf region.In fact, Eisner et al. (in-press) reported that windinduced mixing during the late summer to autumn positively affected the chl a and proportion of larger-size phytoplankton in years with a lower ice extent in the southern Bering Sea.In contrast, a positive but relatively small contribution of the SST to the APP also dominated in the shelf region; its magnitude was relatively lower than the other two factors.The positive contribution of a warmer temperature to the APP was reported in Mueter et al. (2009) and Brown et al. (2011Brown et al. ( , 2013) ) in the Bering Sea.The temperature can positively and directly affect the photosynthetic activity of phytoplankton (Eppley et al., 1972), and the temperature was an indirect effect of an early sea ice retreat.Nevertheless, other controlling factors would be more important in high-latitude regions because of the small seasonal variability of the SST compared to that in mid-latitude seas.In particular, the SST contributed rather negatively to the APP along the shelf break of the Bering Sea (Fig. 6c), where the partial regression coefficient of F L was positive (Fig. 6b).This may be because a colder SST during the MIZ bloom period promotes a larger F L with high productivity (discussed in 4.2).Because the Arctic Ocean is predicted to become a more ice-free ocean by many models (Perovich and Richter-Menge, 2009), it is suggested that the contribution of the nutrient conditions and the subsequent phytoplankton size structure to the APP can be larger in the future.Therefore, understanding the relationship between the timing of the sea ice retreat and the phytoplankton community size composition during the bloom period might advance our knowledge of the biogeochemical cycles and energy transfer to higher-trophic-level organisms in this area.
To the best of our knowledge, this was the first study to capture the relationship between the phytoplankton size structure and the sea ice dynamics using satellite remote sensing data.Although satellite remote sensing data contain large observation errors compared to in situ measurements, such data represent a powerful tool to determine the spatial variability of the response of the primary producers to environmental change.In addition, the continuous observation period for ocean colour remote sensing currently exceeds 17 years over the terms of several sensors (e.g.SeaWiFS, MODIS, and MERIS).Spatio-temporal statistical analysis using a long time-series of ocean colour data can quantify the response of phytoplankton communities to environmental change, and such statistical generalization is important to predict future ecosystem responses or conduct risk assessments.
To conclude, this study revealed that the timing of the sea ice retreat can alter the phytoplankton size composition during the MIZ bloom period.Specifically, an earlier sea ice retreat was linked to larger proportions of large phytoplankton taxa in colder water during cold MIZ bloom periods, and vice versa, in the shelf regions of the Bering and Chukchi seas.Our findings could potentially be linked with the ecosystem changes reported in this region (Grebmeier, 2012 and references), considering the bottom-up energy transfer process via the food web and predator-prey body size relationships.It has been reported that algal grazers greatly depend on the spring bloom as their energy source.For example, the spring-bloom timing, sea ice, and temperature significantly affects the recruitment of copepods in the Bering Sea (Hunt et al., 2002(Hunt et al., , 2011)).The food quality during the spring bloom, which varies with the species composition, also has a significant effect on the secondary production (Søreide et al., 2010;Leu et al., 2011).Thus, comprehending the relationships between the sea ice dynamics, phytoplankton community structure, and environmental conditions is important for assessing the production, activity, or recruitment of highertrophic-level organisms in the northern Bering and Chukchi seas.
There are slight differences in the observation bands between SeaWiFS and MODIS; MODIS measures R rs (488) and R rs (667), while SeaWiFS measures R rs (490) and R rs (670).The R rs (λ) values between them also differed even though they were obtained on the same day and at the same location because of some uncertainties (i.e.atmospheric correction, sensor sensitivity, zenith angle).Therefore, bias correction between them was conducted before analysing and processing the continuous time-series data.We calculated the bias and RMSE values between SeaWiFS-R rs (λ) and MODIS-R rs (λ) (λ = 412, 443, 490, 555, and 667 nm), changing the correction factor from 0.5 to 1.5 to linearly convert SeaWiFS-R rs (λ) into MODIS-R rs (λ) (Fig. A1).The bias and RMSE calculations were conducted using level-3 daily 9 km matched pixels from 2003 to 2007.Then, the optimum correction factors for each band were determined to remove the biases (listed in Table A1).To evaluate the effect of the conversion, we calculated the biases of ocean colour products (i.e.F L , chl a, and a ph ( 443)) between MODIS and the default SeaWiFS, and converted the SeaWiFS values using the daily matched images during 2003-2007 for the study area.The frequency of the biases and their statistics are shown in Fig. A2a-c.Every product was underestimated by the default SeaWiFS relative to MODIS, although the biases were almost completely removed after the SeaWiFS-R rs (λ) conversion.Then, simple averaging between the SeaWiFS and MODIS products was conducted for overlapping observation periods to merge them and increase the number of observation points.412, 443, 490, 555, and 667, respectively.The conversion factors for each band were determined as the optimum factor when the bias was equivalent to zero (listed in Table A1).Note that both bias and RMSE were calculated for logtransformed R rs (λ).A large number of a ph (λ) values were underestimated using the default QAA (Lee et al., 2002(Lee et al., , 2009)), especially for longer λ values (i.e.λ > 490 nm).This was probably due to the uncertainty and error contained in the analytical and empirical equations for deriving a ph (λ).We attempted to overcome this problem by tuning the calculation step of QAA.QAA retrieves a ph (λ) by subtracting the absorption coefficients of detritus+CDOM (a dg (λ)) and pure sea water (a w (λ)) from the total absorption coefficient (a t (λ)): where S dg is the spectral slope of a dg (λ), and it is derived as follows, S dg = S init + 0.002 0.6 + r rs (443)/r rs (555) , (B3) where S init = 0.015 for the default QAA.The underestimation of a ph (λ) was due to over-subtracting a dg (λ) and a w (λ) from a t (λ).In other words, during the derivation of a ph (λ), the default retrieval of a dg (λ) was overestimated.To avoid such an underestimation of a ph (λ), we tuned QAA by determining the optimum S init for the study area.The optimization was performed using the in situ-collected R rs (λ) and a ph (λ).In situ data were obtained from the five cruises of R/V Mirai (Shimada, 2008;Kikuchi, 2009Kikuchi, , 2012;;Itoh, 2010;Nishino, 2013) and three cruises of TR/S Oshoro-maru conducted in the Bering Chukchi seas in 2007-2013 (N = 220).We computed the RMSE between the in situ and modelled a ph (λ) values by changing S init from 0.01 to 0.03 in 0.0005 increments (Fig. B1).Then, the optimum S init was determined to be the value that minimized the average RMSE (S init = 0.019).as a function of the initial spectral slope of a dg (S init ).Lines coloured in blue, sky blue, green, yellow, and grey denote a ph (412), a ph (443), a ph (490), a ph (555), and the average of them, respectively.The optimum value of S init was determined when the RMSE for the average was at a minimum (i.e.0.019).Note that the RMSEs were calculated for log-transformed a ph (λ).

Figure 1 .
Figure 1.Location of in situ sampling stations for the IPY cruises (red), GRENE cruises (violet), and BASIS cruises (blue).Note that the stations of the daily matches between the in situ F L measurements and MODIS level-2-derived F L values are shown in large circles, and crosses represent the sites where primary productivity measurements were conducted.Contours indicate the 100 m, 200 m, and 1000 m bathymetry lines, based on ETOPO-2.
Figure 3. (a) In situ water-column-integrated F L vs. surface F L (r 2 = 0.67, p < 0.01, N = 48), (b) vertical profiles of F L which are normalized at surface F L , (c) in situ PP eu vs. in situ surface PP (r 2 = 0.50, p < 0.01, N = 56 in log scale), (d) depth of PP maximum (Z Pmax ) vs. depth of chl a maximum (Z Cmax ), and (e) vertical profiles of PP which are normalized at surface PP, respectively.Note that the vertical axis of Fig. 3e indicates %PAR relative to surface PAR value.Vertical profiles of F L are coloured by the surface F L value (Fig. 3b).Dashed lines indicate 1 : 1 line, and red solid lines represent regression lines, respectively.

Figure 4 .
Figure 4. Distribution of 16-year median of (a) TSR, (b) date of CMAX, and (c) date difference between CMAX and TSR.Note that the colour scales of Fig. 4a and b are indicated on the upper tick marks of the colour bar, and the scale of Fig. 4c is indicated on the lower tick marks.The arrow on the colour bar denotes the date of the summer solstice (day of year: 174).Contours indicate 50 m and 100 m bathymetry lines.

Figure 5 .
Figure 5. Spatial distributions of Spearman's rank correlation (ρ) between (a) F L during MIZ period and TSR, (b) SST during MIZ period and TSR, (c) OHC during MIZ period and TSR, and (d) PAR during MIZ period and TSR.Yellow indicates a significantly positive ρ (p < 0.05); red, a positive ρ (p > 0.05); blue, a negative ρ (p > 0.05); and light blue, a significantly negative ρ (p < 0.05).The proportions of the coefficients are listed in Table 3. Contours indicate 50 m and 100 m bathymetry lines.

Figure 6 .
Figure 6.Distributions of the partial regression coefficient from Eq. (7) for (a) standardized length of open-water period, (b) standardized annual median of F L , and (c) standardized annual median of SST.The distribution of the p-value as a result of the F test for Eq.(7) is shown in Fig. 6d.Note that the colour scales in Fig. 6a-c are indicated on the upper tick marks on the colour bar, and the scale in Fig. 6d is indicated on the lower tick marks.Contours indicate 50 m and 100 m bathymetry lines.

Figure A1 .
Figure A1.Bias (left axis) and RMSE (right axis) between MODIS-R rs (λ) and converted SeaWiFS-R rs (λ) to MODIS-R rs (λ) as a function of conversion factor.Lines coloured in blue, sky blue, green, yellow, and red indicate λ = 412, 443, 490, 555, and 667, respectively.The conversion factors for each band were determined as the optimum factor when the bias was equivalent to zero (listed in TableA1).Note that both bias and RMSE were calculated for logtransformed R rs (λ).

Figure A2 .
Figure A2.Histograms and statistics of the bias between SeaWiFSand MODIS-retrieved (a) F L , (b) chl a, and (c) a ph(443).Red coloured lines and statistics denote the results from default SeaWiFS-R rs (λ), and blue coloured lines and statistics denote the results from converted SeaWiFS-R rs (λ).The products were significantly improved (p < 0.01) for means (t test) and medians (U test).Note that biases of chl a and a ph (443) were computed for the logtransformed values.

Figure B1 .
Figure B1.RMSE between in situ a ph (λ) and QAA-derived a ph (λ) as a function of the initial spectral slope of a dg (S init ).Lines coloured in blue, sky blue, green, yellow, and grey denote a ph (412), a ph (443), a ph (490), a ph (555), and the average of them, respectively.The optimum value of S init was determined when the RMSE for the average was at a minimum (i.e.0.019).Note that the RMSEs were calculated for log-transformed a ph (λ).

Table 1 .
List of symbol and units.