Anomalous wind triggered the largest phytoplankton bloom in the oligotrophic North Pacific Subtropical Gyre

In summer 2010, a massive bloom appeared in the middle (16–25°N, 160–200°E) of the North Pacific Subtropical Gyre (NPSG) creating a spectacular oasis in the middle of the largest oceanic desert on Earth. Peaked in June 2010 covering over two million km2 in space, this phytoplankton bloom is the largest ever recorded by ocean color satellites in the NPSG over the period from 1997 to 2013. The initiation and mechanisms sustaining the massive bloom were due to atmospheric and oceanic anomalies. Over the north (25–30°N) of the bloom, strong anticyclonic winds warmed sea surface temperature (SST) via Ekman convergence. Subsequently, anomalous westward ocean currents were generated by SST meridional gradients between 19°N and 25°N, producing strong velocity shear that caused large number of mesoscale (100-km in order) cyclonic eddies in the bloom region. The ratio of cyclonic to anticyclonic eddies of 2.7 in summer 2010 is the highest over the 16-year study period. As a result of the large eddy-number differences, eddy-eddy interactions were strong and induced submesoscale (smaller than 100 km) vertical pumping as observed in the in-situ ocean profiles. The signature of vertical pumping was also presented in the in-situ measurements of chlorophyll and nutrients, which show higher concentrations in 2010 than other years.

In this study, we report the largest phytoplankton bloom ever observed by ocean color sensors in the NPSG. The bloom occurred in summer 2010 and appeared as a long zonal band along 20°N between Taiwan and Hawaii. Higher CHL concentrations resided within the 160°E and 200°E band (Fig. 1). Using a combination of satellite remote sensing and in-situ measurements, we showed that the eastward propagating bloom was directly related to eddies which induced vertical pumping as a result of anomalous wind-triggered westward currents.

Results
The largest summer bloom in the central NPSG from 1998 to 2013. Located in between Taiwan and Hawaii (Fig. 1a), the bloom started as a long zonal band feature with high CHL concentration waters emerging roughly within 15-25°N and 165-185°E. Chronologically, the bloom first occurred at the end of May and propagated eastward. By the third week of June, the bloom has reached its peak covering almost 40° in longitude. After that, the bloom began to decline but still propagating further east to about 200°E closer to Hawaii in July 2010 (Fig. 1c). CHL concentration within the bloom ranged from 0.1 to 0.3 mg/m 3 (Fig. 1a). At its peak, the CHL concentration can reach as high as 1 mg/m 3 (Fig. 1b). Even when the CHL concentration decreased, the zonal band feature of the bloom can still be seen clearly in between Taiwan and Hawaii (Fig. 1c).
Considering only the region of interest (i.e. black box in Fig. 1) and only for the month of June and July, the mean CHL concentration (about 0.08 mg/m 3 ) observed in 2010 is the highest in the period from 1998 to 2013 (Fig. 2a). Also, when considering only the area with CHL concentrations higher than 0.1 mg/m 3 , the bloom's area is the largest in 2010 reaching around two million km 2 . On the month-to-month time scale from May to July of 2010, Fig. 3a,b show the 8-day variability of CHL concentration and of high-CHL-concentration area averaged over the region of interest, respectively. The bloom first occurred at the end of May, peaked at mid-June and vanished by the end of July, 2010. From Fig. 3a,b, it is clear that CHL concentrations and high-CHL-concentration area are significantly higher than the climatological mean from 1998 to 2013.

The anomalies of oceanic fields and winds in 2010.
In the NPSG, anticyclonic eddies are generally more abundant than cyclonic eddies in total number in the STCC zonal band from 122°E to 170°E between 12°N and 28°N, except between 19°N and 22°N 19 . In summer, however, the number of cyclonic eddies are generally more than the anticyclonic eddies ( Fig. 2d) within the bloom region (160-200°E and 16-25°N). Particularly in 2010, the number of cyclonic eddies were 2.7 times more than the anticyclonic eddies ( Fig. 2d), which is the highest from 1998 to 2013, showing the dominance of cyclonic eddies. Moreover, the cyclonic and anticyclonic eddies have similar geostrophic-current speeds (Fig. 2e) and areas (Fig. 2f) in 2010. On the month-to-month time scale, Fig. 3c shows an increasing in the cyclonic-eddy number and a decreasing in the anticyclonic-eddy number in June 2010. The huge differences in the cyclonic-anticyclonic eddy numbers in 2010 are significant in comparison to the climatology of eddy numbers (Fig. 3c). In addition, the size of the cyclonic and anticyclonic eddies was also anomalously smaller during the bloom-peak period in June when comparing to the climatological eddy fields (Fig. 3d). The coincidence of eddy size and eddy number with the bloom occurrence suggests a bloom mechanism related to the ocean eddies in summer 2010.
Since eddy variability is mainly determined by the changes of STCC 1,3 and HLCC 6,7 , the background currents were examined to study the possible mechanism for the eddy-number difference in 2010. Figure 4 shows the background zonal currents represented by the satellite-obtained absolute geostrophic currents and relative velocity horizontal shear in the region of interest during June and July of 2010, comparing to the climatological fields from 1998 to 2013. The STCC and HLCC can be seen in Fig. 4a (dashed curve), centered along 24-25°N and 19-20°N 20 , respectively, flowing eastward at a speed reaching 2.5-5 cm/s. Between these subtropical countercurrents, a westward current flows at a maximum speed of about 2.5 cm/s, centered along 22°N. In Fig. 4b (dashed curve), current velocity shear is mostly positive except within the 20°N and 22°N region.
In contrast to the climatological fields, anomalous westward currents are found at 25°N (Fig. 4a,c) in 2010, inducing anomalous changes in current velocity shear (Fig. 4b,d) that could change potential vorticity and cause eddy polarity-asymmetry (more cyclonic eddies and less anticyclonic eddies) to the south. According to the quasi-geostrophic theory 21 , there is a necessary condition for barotropic instability or baroclinic instability with a shear flow, when the dominant sign change is mainly due to horizontal shear or vertical shear. Thus, the shear instability induced by the anomalous westward currents could contribute to the eddy polarity-asymmetry during the bloom period.
Oceanic background flow is determined by surface winds over the NPSG 3,20,22 . In spring, strong cyclonic wind anomalies in the north of the eastward-flowing STCC enhanced the SST front by decreasing SST via Ocean vertical profiles in 2010. Large number difference in cyclonic and anticyclonic eddies can generate strong eddy-eddy interactions. The strength of eddy-eddy interactions was derived using the Okubo-Weiss (OW) parameter [23][24][25] , which can provide the difference between deformation rate and vorticity. The OW-parameter averaged within the region of interest in summer 2010 is about 1 × 10 −13 s −2 , showing stronger eddy deformation than eddy rotation. The OW-parameter switched from negative to positive values at the end of May, and then was in positive values in June (Fig. 3e). This shows that the eddy deformation rate is larger than the eddy vorticity during the bloom period, suggesting strong deformation fields among the eddies. Under these mesoscale circumstances, submesoscale vertical pumping can be induced in the high-velocity regions surrounding both types of eddies 26 , via horizontal stretching done by eddy deformation fields during eddy-eddy interactions 27 .
To study the relevant ocean vertical structure, we used publicly available hydrographic profiles observed by conductivity-temperature-depth devices (CTD), Argo buoys and expendable bathythermographs (XBT) that passed through the bloom region. Locations and platform numbers/names of the profiles were superimposed on the maps of CHL concentrations and sea-surface-height anomalies (SSHAs) from June 9 to June 26 (Fig. 6), a duration that covers the whole observation periods of CTD, Argo buoys and XBT. Locations of the profiles corresponding to the regions with high CHL concentrations and eddy cores can be observed in Fig. 6. An XBT transect passed through the bloom area at about 18-20°N (Fig. 6a) during the growing phase of the bloom (Fig. 3a), whereas a meridional CTD transect at 165°E was available during the decline phase of the bloom (Fig. 6c). www.nature.com/scientificreports www.nature.com/scientificreports/ Two Argo buoys with platform numbers of 5901596 (the most southwest red dot in Fig. 6) and 2900509 (the most southeast red dot in Fig. 6) observed the background vertical structures corresponding to low CHL concentrations ( < 0.05 mg/m 3 ) observed by satellites, but the others with platform numbers of 5901538, 5901420, 5901419 and 5901539 were in the bloom region with CHL concentrations higher than >0.1 mg/m 3 .
Inside the bloom region, Fig. 7a-c,e show that the buoyancy frequencies were lower below 50 m and the mixed layer depth is shallower than outside the bloom region shown in Fig. 7d,f. Moreover, the mixed layer stratification is weaker during the bloom growing period in the first half of June than during the declining period after middle of June. These comparisons to the background fields show that the bloom occurrence was under an ocean condition that the mixed-layer stratification was weak and shallow, favorable for nutrient uplifting to the sunlit layer.   www.nature.com/scientificreports www.nature.com/scientificreports/ (TS) diagram in Fig. 8 shows that the water properties corresponding to the high and low CHL concentrations are significantly different. The TS properties where the satellite-observed CHL concentrations are high (>0.1 mg/ m 3 ) were colder and saltier than the waters with lower CHL concentrations (<0.05 mg/m 3 ) and climatology. Statistically over 50 m (roughly the mixed-layer depth), the temperature in the high-CHL regions is about 0.57 °C colder than that in the low-CHL regions; while the density and salinity in the high-CHL regions is about 0.28 kg/ m 3 and 0.13 psu larger than those in low-CHL regions, respectively. These comparisons suggest the results of vertical pumping by bringing up colder and denser waters to the mixed layer from below. Figure 7g shows the temperature profiles observed by the XBT on the WDD6033 platform passing through the bloom region from east to west during June 9 to 11 (see Fig. 6a for corresponding locations of XBT stations). Shown by the concave upward 24°C isotherm, vertical pumping/upwelling can be observed at 172.5°E, 176°E and 179.5°E. The depth of 24°C waters can reach as shallow as 80 m at 172.5°E and 179.5°E, having a horizontal length scale of submesoscale (smaller than 100 km), located at where the eddy currents are strong near boundaries between the pairs of mesoscale cyclonic and anticyclonic eddies (see Fig. 6a for corresponding locations). The depths of isotherm at 26°C, 24°C and 20°C are shallower at 176°E than the vicinity, located near the boundary of a mesoscale cyclonic eddy (see Fig. 6a for corresponding locations).   Fig. 6c for corresponding locations), indicating vertical upward motion where the eddy current speed was fast. Figure 9 shows the concentration profiles of CHL, nitrate + nitrite, phosphate and silicate along the 165°E transect in June 2010 and climatology fields from 2003 to 2012. The deep CHL maximum layer was at about 100-150 m in 2010. Higher CHL concentrations can be found at the 18° and 20°N stations when compared to the climatology profile. Within 20°N and 21°N, higher concentration of nutrients can be observed (Fig. 9) at where the isotherm outcropped (Fig. 7h). Noting that, Figs 3a and 6c show the nutrient sampling nearly outside the bloom region during the decline of bloom. Thus, the nutrient concentration is not remarkably high as shown in Fig. 9, but the occurrence of upwelling still can be seen in the vertical distribution of nutrients.

Summary
In this study, we reported an enormous phytoplankton bloom in the eddy-abundant NPSG during summer 2010. The bloom is the largest in area ever recorded by ocean color satellites in the NPSG over the period from 1997 to 2013, determined by the top-down processes of atmosphere and ocean from large scale to submesoscale.
Started from the large-scale process of air-sea interaction north of the bloom, strong anticyclonic winds warmed SST subsequently causing positive meridional SST gradients and inducing anomalous westward flows of ocean surface in the NPSG. Ocean mesoscale variability (100 km in order) responded to the velocity shear of the anomalous westward flows, by increasing the number of cyclonic eddies. This resulted in higher ratio of cyclonic to anticyclonic eddies. Subsequently, eddy deformation became more dominant than eddy rotation during the bloom period, causing strong eddy-eddy interactions. As a consequence of the strong eddy-eddy interactions submesoscale vertical upwelling (smaller than 100 km) was generated, bringing up nutrients below the sunlit layer to euphotic depth.
The submesoscale vertical upwelling and relevant ocean vertical structures favorable for nutrient uplifting were observed using available historical in-situ ocean profile observation obtained from Argo buoys, XBT and CTD sampling that passed through the bloom. The ocean stratification above 200 m was the weakest in summer 2010, by comparing to other years. At the beginning of the bloom occurrence, the mixed layer was shallower with weaker stratification inside the bloom region than outside. Above the mixed layer depth, colder and denser waters were observed at where the CHL concentrations were high during the bloom period. The submesoscale vertical upwelling was observed where eddy flows were fast near eddy boundaries. The corresponding higher nutrients and CHL observed from in-situ measurements show that a combination of atmospheric and oceanic anomalies is responsible for the largest bloom in the NPSG.

Data and Methods
We used 8-day CHL data at 1/4° resolution, available since 1997, obtained from the European Space Agency's GlobColour project (http://www.globcolour.info), which provides a merged CHL product with better coverage from sensors such as Sea-Viewing Wide Field-of-View Sensor (SeaWiFS), Medium Resolution Imaging Spectrometer (MERIS), Moderate Resolution Imaging Spectroradiometer (MODIS) onboard Aqua satellite, and Visible Infrared Imaging Radiometer Suite (VIIRS). We applied the daily altimetry SSHAs and absolute dynamic topographic (ADT) at 1/4° resolution obtained from the Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) data (http://www.aviso.altimetry.fr), to study the changes of SSHAs and background www.nature.com/scientificreports www.nature.com/scientificreports/ flows. Provided by Faghmous et al. 25 , the "daily global mesoscale ocean eddy dataset" from satellite altimetry was used to analyze eddy variability, available from 1993 to 2013.
For data consistency, the daily data of SSHAs, ADT and eddies were averaged at 8-day intervals according to the dates of CHL data. Moreover, we analyzed the data from 1998 to 2013 that covered the duration of eddy dataset and CHL data. The geostrophic velocity anomalies (u′, v′) and total geostrophic velocities (u, v) were then derived from the zonal and meridional gradients of the 8-day SSHAs and ADT 2 , respectively. The eddy variability was studied based on the given data of eddy numbers, eddy area, eddy geostrophic-current speed. Only the eddies with amplitude larger than 4 cm were considered due to the SSH accuracy of satellite altimeter of about 4 cm in the open ocean 28 . Faghmous et al. 25 defined eddies as the outermost closed-contour sea level anomalies (SLA) containing a single extreme (maximum/minimum). Bounded by the outermost closed-contour SLA, the eddy area was estimated as the surface area and the eddy geostrophic-current speed was the geostrophic-current speed averaged in the eddies. There were two uncertainties associated with the eddy-identification algorithm of Faghmous et al. 25 : (1) the eddy boundary was not necessarily associated with the eddy physical properties, and (2) small features might not be detected. However, the eddy-identification algorithm recovered about 96% of eddy features as identified by domain experts 25 .
To study eddy-eddy interactions, the Okubo-Weiss (OW) parameter [23][24][25] was calculated using where S sh is the shear-deformation rate estimated as ∂v′/∂x + ∂u′/∂y, S st is the stretch-deformation rate estimated as ∂u′/∂x − ∂v′/∂y and ξ is the vorticity estimated as ∂v′/∂x − ∂u′/∂y, using the geostrophic velocity anomalies (u′, v′). If the eddy-eddy interactions are strong, one particular circular eddy can be stretched by the deformation field of the nearby eddies 27,29 , causing the deformation rate (S sh 2 + S st 2 ) larger than the vorticity (ξ 2 ) and thus getting positive OW-parameter. Conversely, when the eddy-eddy interactions are weak, the stretched eddy recovers to a circular shape by increasing its vorticity. Thus, the deformation rate becomes smaller than the vorticity, having the negative OW-parameter.
To observe ocean vertical structure during the bloom occurrence, we applied the ocean-profile data of CTD, XBT and Argo buoys. These in-situ ocean-profile data were collected and made available by the "Coriolis project and programmes" (http://www.coriolis.eu.org). Only those Argo-buoy data with good-flagged quality were applied. For the bio-chemical parameters, we analyzed the publicly available concentration profiles of CHL, nitrate + nitrite, phosphate and silicate measured along the 165°E transect by the Japan Meteorological Agency (JMA). CHL were measured fluorometrically using a Turner Designs 10-AU-005-CE Field Fluorometer and