Satellite views of the seasonal and interannual variability of phytoplankton blooms in the eastern China seas over the past 14 yr ( 1998 – 2011 )

The eastern China seas are some of the largest marginal seas in the world, where high primary productivity and phytoplankton blooms are often observed. However, little is known about their systematic variation of phytoplankton blooms on large spatial and long temporal scales due to the difficulty of monitoring bloom events by field measurement. In this study, we investigated the seasonal and interannual variability and long-term changes in phytoplankton blooms in the eastern China seas using a 14 yr (1998– 2011) time series of satellite ocean colour data. To ensure a proper satellite dataset to figure out the bloom events, we validated and corrected the satellite-derived chlorophyll concentration (chla) using extensive in situ datasets from two large cruises. The correlation coefficients between the satellite retrieval data and the in situ chl a on the logarithmic scale were 0.85 and 0.72 for the SeaWiFS and Aqua/MODIS data, respectively. Although satellites generally overestimate the chla, especially in highly turbid waters, both the in situ and satellite data show that the overestimation of satellitederived chla has an upper limit value (10 μg L −1), which can be used as a threshold for the identification of phytoplankton blooms to avoid the false blooms resulting from turbid waters. Taking 10 μg L −1 as the threshold, we present the spatial-temporal variability of phytoplankton blooms in the eastern China seas over the past 14 yr. Most blooms occur in the Changjiang Estuary and along the coasts of Zhejiang, with a maximal frequency of 20 % (about 73 days per year). The coasts of the northern Yellow Sea and the Bohai Sea also have high-frequency blooms (up to 20 %). The blooms show significant seasonal variation, with most occurring in spring (April–June) and summer (July–September). The study also revealed a doubling in bloom intensity in the Yellow Sea and Bohai Sea during the past 14 yr. The nutrient supply in the eastern China seas might be a major controlling factor in bloom variation. The time serie in situ nutrient datasets show that both the nitrate and phosphate concentrations increased more than twofold between 1998 and 2005 in the Yellow Sea. This might be the reason for the doubling of the bloom intensity index in the Yellow Sea and Bohai Sea. In contrast, there has been no significant long-term increase or decrease in the Changjiang Estuary, which might be regulated by the Changjiang River discharge. These results offer a foundation for the study of the influence of phytoplankton blooms on the carbon flux estimation and biogeochemical processes in the eastern China seas. Published by Copernicus Publications on behalf of the European Geosciences Union. 4722 X. He et al.: Seasonal and interannual variability of phytoplankton blooms


Introduction
The eastern China seas, comprising the East China Sea (ECS), Yellow Sea (YS) and Bohai Sea (BS), are some of the largest marginal seas in the world (Fig. 1a).The seas receive a tremendous amount of fresh water and suspended sediments from rivers, including two of the largest rivers in the world, the Changjiang (Yangtze River) and Huanghe (Yellow River) (Liu et al., 2010a).The Changjiang provides 90 % of the total riverine water discharge and sediment input to the ECS (Zhang et al., 2007).The annual mean discharge of the Changjiang is 29 300 m 3 s −1 , with 70 % of the total annual runoff occurring during the flood season from May to October (Wang et al., 2010).Due to the sediment carried by rivers, the eastern China seas are among the most turbid ocean regions in the world (Shi and Wang, 2010).
The eastern China seas are some of the most productive areas in the North Pacific owing to the abundant nutrients from river discharge and the incursion of the Kuroshio subsurface water on the shelf (Chen et al., 1995;Chen, 1996;Chen and Wang, 1999).Many studies have used cruise data to examine the distributions of chlorophyll concentration (chl a) and primary production in the eastern China seas (Gong et al., 1996(Gong et al., , 2000(Gong et al., , 2003(Gong et al., , 2006;;Hama et al., 1997;Chen et al., 2001Chen et al., , 2004;;Kim et al., 2009).The modelled annual mean primary production over the entire ECS shelf is up to 441 mg C m −2 d −1 (Liu et al., 2010b), and episodic phytoplankton blooms are often observed in the eastern China seas (Hyun and Kim, 2003;Furuya et al., 2003;Hu et al., 2004;Son et al., 2006;Tang et al., 2006;Isobe and Matsuno, 2008;Zhao et al., 2011).Phytoplankton blooms play a key role in the marine ecology, primary production, biogeochemical processes and marine carbon cycle, from fuelling a rich and diverse food web and associated fishery harvests (Ware and Thomson, 2005) to driving rapid fluctuations in the air-sea flux of carbon dioxide (Evans et al., 2011).Blooms also affect the severity and spatial extent of hypoxic zones (Grantham et al., 2004).However, due to the limitation of field investigations, insufficient data have been collected by traditional sampling methods to assess the spatial and temporal dynamics of phytoplankton blooms.
Satellite ocean colour data provide synoptic, long-term time series coverage, which is ideal for the identification of bloom events (McKibben et al., 2013).As a practical index of phytoplankton biomass, the spatial and temporal variability of chl a has been widely investigated in the eastern China seas using satellite-derived chl a products.Ning et al. (1998) and Tang et al. (1998) studied monthly variations in pigment concentrations measured by the Coastal Zone Color Scanner (CZCS) on board the Nimbus-7 satellite.Kiyomoto et al. (2001) examined the seasonal and spatial distributions of chl a and suspended particulate matter concentration using satellite data from the CZCS and OCTS (Ocean Color and Temperature Scanner) on board the Advanced Earth Observing Satellite (ADEOS).Recently, using 8 yr (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) of observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the Aqua satellite, Shi and Wang (2012) presented climatologically seasonal distributions of satellite-derived chl a in the eastern China seas.Yamaguchi et al. (2012) also investigated the seasonal climatology of satellite-derived chl a in the YS and ECS using 10 yr (September 1997to October 2006) of chl a data observed by the Sea-viewing Wide Field-of-view Sensor (Sea-WiFS).Although Shi and Wang (2012) and Yamaguchi et al. (2012) analysed the spatial and temporal variability of chl a in the eastern China seas, seasonal and interannual variations in phytoplankton blooms on large spatial and temporal scales have never been investigated.Moreover, little is known about the long-term changes in phytoplankton blooms in the eastern China seas, though the marine environment of China has deteriorated significantly over recent decades due to the effects of rapid economic development and increased population, and the associated consequences of habitat loss, eutrophication, pollution and overfishing (see, for example, Zhang et al., 1999;Jin, 2004;Liu and Diamond, 2005;Dong et al., 2010).
Satellite-derived chl a data can be used as a practical index for identifying phytoplankton blooms (Siegel et al., 2002;Kim et al., 2009); however, ocean colour satellite data generally overestimate chl a in turbid coastal waters, which causes false blooms to be counted.The operational satellite retrieval algorithm for chl a uses the blue-green band ratio of remote sensing reflectance, with low ratios corresponding to high chl a (O' Reilly et al., 1998).Additional terrestrial detritus and coloured dissolved organic matter (CDOM) can significantly increase the water's absorption coefficient and decrease the remote sensing reflectance in the blue band, while it enhances the particulate backscattering coefficient and increases the remote sensing reflectance in the green band.Hence, the chl a is overestimated by the blue-green band ratio algorithm (Dierssen, 2010).For this reason, Shi and Wang (2012) did not present the interannual variation of the satellite-derived chl a in the eastern China seas.Yamaguchi et al. (2012) only assessed the interannual variations of the satellite-derived chl a in summer with less water turbidity in the eastern China seas.Thus, a corrected satellite chl a product and a practical and credible chl a threshold for bloom identification are needed as the basis on which to exam the long-term bloom variation from space.
In this study, we first validate and correct the satellitederived chl a based on an abundant in situ dataset from two large cruises (up to 1255 samples) in the eastern China seas.Second, we exam and determine a credible chl a threshold for bloom identification in these turbid coastal oceans.Then, using a 14 yr (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) time series of corrected SeaW-iFS and Aqua/MODIS chl a data, we investigate the seasonal and interannual variability and long-term changes in phytoplankton blooms in the eastern China seas.Finally, we discuss the possible factors that control the long-term changes in blooms.

Study area
The study area includes the ECS, YS and BS, comprising a total area of about 1 million km 2 and an average depth of 63 m from the coastal line to the 200 m isobath (Liu et al., 2010a) (Fig. 1a).The BS is an inland shallow sea with an average depth of 18 m (Liu et al., 2003a).The YS is located between mainland China and the Korean Peninsula, with its mouth opening into the ECS.The ECS, the eleventh largest marginal sea in the world, extends from Cheju Island in the north to Taiwan Island in the south, and is bordered by the Kuroshio Current on the east and the coast of China on the west.
The circulation in the eastern China seas has been extensively investigated (see, for example, Beardsley, 1985;Chao, 1990;Fang et al., 1991;Hu, 1994;Su, 1998;Katoh et al., 2000;Naimie et al., 2001).The strong western ocean boundary current Kuroshio flows northeastward along the ECS shelf at water depths of 200 to 1000 m.The intrusion of the Kuroshio off northeastern Taiwan combines with the northward flow of the northward Taiwan Warm Current (TWWC) from the Taiwan Strait in summer (Ichikawa and Beardsley, 2002), and from the Kuroshio in winter (Chen and Sheu, 2006).In addition, as the Kuroshio turns eastward toward the northwestern Pacific, a branching current intrudes onto the shelf and penetrates into the Yellow Sea and even the Bohai Sea, forming the Yellow Sea Warm Current (YSWC, Lie et al., 2001).In summer, the Changjiang plume (or the Changjiang Diluted Water, CDW) extends northeastward to the shelf, broadly forced by the southwest monsoon (Chang and Isobe, 2003).

Satellite ocean colour datasets
The 8-day composite chl a and remote sensing reflectance at 555 nm (Rrs555) retrieved by SeaWiFS and Aqua/MODIS were acquired from the NASA ocean colour website (http://oceancolor.gsfc.nasa.gov).The spatial resolution of the datasets is 5 with global coverage.The SeaWiFS dataset spans January 1998 to December 2007, and the Aqua/MODIS dataset spans January 2003 to December 2011.Although SeaWiFS worked successfully until the end of 2010, data are unavailable for several months during 2008 to 2010 due to technical difficulties with the sensor.In this study, we needed to calculate the number of bloom events in each year.With several months of data lacking in one year, the results would be affected by the total number of valid samples.Therefore, we only used the continuous Sea-WiFS data from 1998 to 2007.We used the 8-day composite chl a datasets to obtain the statistics on phytoplankton blooms.In addition, we used the satellite-derived Rrs555 to represent the water turbidity, with larger Rrs555 corresponding to higher water turbidity (Yamaguchi et al., 2012).
For the validation and correction of the satellite-derived chl a data, we obtained the level 2 (L2) ocean colour data with native high resolution for both the SeaWiFS and Aqua/MODIS data from the NASA ocean colour website (http://oceancolor.gsfc.nasa.gov),corresponding to the in situ data period.

In situ datasets of chl a
In situ chl a datasets from two large cruises covering all of the eastern China seas (Fig. 1b-c) were collected in this study to validate and correct the satellite-derived chl a data.These cruises were initiated by the State Oceanic Administration (SOA) of China, and five vessels were used to conduct the investigations synchronously throughout the eastern China seas.The summer cruises were undertaken from 15 July to 31 August 2006, and had 662 stations.The winter cruises were from 19 December 2006 to 4 February 2007, and had 593 stations.The two cruises used almost the same sampling locations and the same measurement method for the chl a.The water samples for chl a analysis were filtered through 25 mm Whatman GF/F filters with 0.7 µm pore size.The filters were preserved in liquid nitrogen before being processed.Each filter was analysed with a Turner Designs fluorometer to obtain the chl a concentration using the non-acidification fluorometric method (Welschmeyer, 1994).In this study, we only used surface layer chl a data with a measurement depth of less than 5 m.If there were multiple samples within 5 m, we used the topmost surface layer sample.

In situ datasets of nutrients
We downloaded the time series oceanographic survey datasets measured along the coast of South Korea from the NOAA National Oceanographic Data Center (NODC).These datasets originate from the Korea Oceanography Data Center (KODC), and date from 1961 to 2005, with most of data collected between 1994 and 2005.In this study, we used surface layer (zero depth) nitrate and phosphate datasets along three latitudinal sections (35.5130 • N, 34.4300 • N and 32.0000 • N) from 1998 to 2005 (Fig. 1a), corresponding to the time period of the ocean colour satellite datasets.The north section we chose has five time series stations named N1 to N5 from west to east, with N1 in the central YS and N5 near the coast (Fig. 1a).The middle section has three time series stations (named M1, M2 and M3 from west to east), and the south section has five serial stations (named S1 to S5 from west to east).Both the nitrate and phosphate concentrations were regularly measured at these time series stations in February, April, June, August, October and December each year.
In addition, we obtained time series oceanographic survey datasets from the Japan Meteorological Agency website (http://www.data.kishou.go.jp/kaiyou/db/vessel obs/ data-report/html/ship/ship.php).The surface layer (zero depth) nitrate and phosphate datasets from 1998 to 2010 along the famous "PN" section with six time series stations (named K1 to K6 from west to east) (Fig. 1a) were used in this study.At each time series station, both nitrate and phosphate concentrations were regularly measured at least in January, April, July, and October to represent the four seasons.

Normalisation of satellite-derived chl a between SeaWiFS and Aqua/MODIS
There are inevitable differences between the chl a values retrieved by SeaWiFS and Aqua/MODIS due to the differences in sensor performance (e.g.bands, calibration accuracy) and retrieval algorithms (e.g.atmospheric correction and chl a algorithms) (Morel et al., 2007).To remove these systematic differences, we compared the annual mean chl a for each of the pixels derived from SeaWiFS with those derived from Aqua/MODIS during their overlapping periods (2003)(2004)(2005)(2006)(2007), as shown in Fig. 2. As a result, we normalised the SeaWiFSderived chl a by the Aqua/MODIS data according to a linear regression on the logarithmic scale based on the data in Fig. 2: log(Chl N ) = 0.0184 + 1.0760 × log(Chl S ), (1) with R = 0.97, SD = 0.13, P < 0.0001, where Chl S and Chl N are the pre-and post-normalised SeaWiFS-derived chl a, respectively, and R is the correlation coefficient and SD is the standard deviation.

Validation and correction of satellite-derived chl a
We followed the validation method proposed by Bailey and Werdell (2006).Specifically, the satellite (SeaWiFS or Aqua/MODIS) match-up dataset comprised the mean of the valid values found for the 5 × 5 pixels within a time window of ±3 h from the in situ observation by the following steps: 1.For each valid in situ chl a sample, extract the valid satellite chl a values (after checking the L2 flags) within a 5 × 5 box centred on the in situ location if the time difference between the satellite overpass and in situ measurement is less than 3 h.
5. The mean of the remaining valid pixels is compared with the corresponding in situ chl a value.
There are 27 and 39 valid match-ups for SeaWiFS and Aqua/MODIS, respectively, as shown in Fig. 3.It should be noted that the SeaWiFS chl a values in Fig. 3 have been normalised by Eq. ( 1).The correlation coefficients between the satellite retrieval and in situ chl a values on the logarithmic scale are 0.85 and 0.72 with SDs of 0.30 and 0.40 for the SeaWiFS and Aqua/MODIS data, respectively.Clearly, both SeaWiFS and Aqua/MODIS overestimate the chl a values in the eastern China seas, especially in highly turbid waters (corresponding to large Rrs555).To eliminate the systematic overestimations, we correct the satellite-derived chl a values with the linear regressions derived from the data in Fig. 3 for SeaWiFS and Aqua/MODIS as follows: where Chl C is the chl a value after correction, with the suffixes SWF and MOD corresponding to SeaWiFS and Aqua/MODIS, respectively; Chl N(SWF) is the SeaWiFSderived chl a value after normalisation with Eq. ( 1), and Chl (MOD) is the Aqua/MODIS-derived chl a value.

Identification of phytoplankton blooms in marginal seas with turbid waters
To date, no uniform method and threshold exists to define phytoplankton blooms due to various phytoplankton communities.For example, Henson and Thomas (2007) and Siegel et al. (2002) used a value 5 % above the annual median values of surface chl a to define the initiation times of blooms.Kim et al. (2009) used two different approaches to define phytoplankton blooms: a constant threshold, defined as two standard deviations above the long-term average, and a varying threshold, depending on each year's median value.Generally, the constant threshold approach is useful for understanding the long-term changes in phytoplankton blooms, whereas the varying threshold approach is aimed at understanding the timing and intensity of blooms on an annual timescale.
To study the long-term variations in phytoplankton blooms, a constant criterion of chl a for bloom identification needs to be defined.There are two difficulties in determining a threshold in the eastern China seas, one is that it is hard to define a uniform threshold for blooms from the literature, and the other is that the overestimated satellite-derived chl a produces "false blooms" in turbid water.Figure 4 shows the climatological annual mean distributions of satellite-derived chl a and their standard deviations after correction, and the annual mean Rrs555 derived from all of the 8-day composite datasets from SeaWiFS (1998SeaWiFS ( -2007) ) and Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).The spatial patterns of the annual mean chl a are quite similar to Rrs555 (representing the water turbidity), with gradual decreases from coast to shelf.In addition, it can be clearly seen that annual mean chl a is very high with less temporal variation (small standard deviation) along the turbid coastal waters year-round, which is probably caused by the significant overestimation by satellite-derived chl a (Yamaguchi et al., 2012;Shi and Wang, 2012).
To investigate how much overestimation of satellitederived chl a could be induced by water turbidity, we derived the maximum chl a and maximum and minimum Rrs555 from all 8-day composite datasets of SeaWiFS (1998SeaWiFS ( -2007) ) and Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), as shown in Fig. 5.In the Hangzhou Bay (HZB), Changjiang River mouth (CJM), Subei shallow (SBS) and Huanghe River mouth (HHM), where waters are extremely turbid, the maximum chl a is generally less than 10 µg L −1 for both SeaWiFS and Aqua/MODIS.This indicates that the overestimation of the satellite-derived chl a induced by water turbidity is less than 10 µg L −1 .From the validation results shown in Fig. 3a and b, we can also see that the chl a values for highly turbid waters with the highest Rrs555 are less than 10 µg L −1 , which is consistent with the result derived from the satellite data.Therefore, 10 µg L −1 is a conservative critical threshold for the identification of phytoplankton blooms in marginal seas with turbid waters.In the following analysis, we took 10 µg L −1 as a threshold to identify blooms in the eastern China seas.This threshold is within the threshold values of major blooms (17.9 µg L −1 ) and minor blooms (5.7 µg L −1 ) used by Kim et al. (2009).
Note that in this study, we only focused on the large blooms in the eastern China seas with satellite-derived chl a higher than 10 µg L −1 .In Sect.4.4, we also test this criterion by comparing the statistical results with thresholds of 8 and 12 µg L −1 to examine the influence of threshold values.

Statistical indices of bloom events
We binned the whole study area into 0.25  defined for various time periods (e.g.seasonal, annual) as follows: (1) The percentage of valid chl a in a grid cell: where N is the total number of projected 8-day composite satellite chl a datapoints in a certain period, e.g.N = 46 per year for one satellite; and I Chl = 1 if there is at least one projected pixel with a valid value (chl a > 0) within the cell, otherwise I Chl = 0.For satellite ocean colour remote sensing, PE is usually determined by the degree of cloud coverage.If the PE of an area is too small (e.g.PE < 30 %), the bloom statistical result could have large uncertainty due to insufficient valid data.
(2) Bloom occurrence frequency (defined as chl a larger than 10 µg L −1 ) in a grid cell: where I bloom = 1 for a bloom event occurring in a grid cell, otherwise I bloom = 0.A bloom event was assumed to have occurred in a grid cell if at least one projected pixel within the cell had chl a greater than 10 µg L −1 .
(3) Bloom intensity index (defined as chl a larger than 10 µg L −1 ) in a specific region: where M is the total number of cells in the specific region, S bloom is the averaged chl a of all bloom events occurring in a cell, and A is the area of a cell.BI is an area-integrated parameter representing the chl a inventory of phytoplankton blooms in the 1 m surface layer in a specific region in units of kg m −1 .

Uncertainties in the statistical results
The greatest challenge for the analysis of interannual variability using satellite ocean colour data is the amount of missing valid data due to cloud coverage.Merging multiple satellite datasets may reduce the effect of cloud coverage, but for long-duration clouds, such improvement is still limited.In this study, to increase the proportion of valid data, we used 8day composite chl a data for bloom analysis instead of daily data.Figure 6a shows the annual climatology of PE (the percentage of valid chl a in a grid cell) derived from all 8-day composite chl a datasets from SeaWiFS (1998SeaWiFS ( -2007) ) and Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).Among the eastern China seas, the BS and northern Yellow Sea (NYS) have the lowest cloud coverage, with PEs larger than 75 %.The PE is also greater than 50 % in the ECS, although the cloud coverage is relatively heavy, especially along the Kuroshio and along the inner shelf.The low PE areas do not affect our result much, because the Kuroshio is oligotrophic water and the inner shelf lacks sufficient radiance for phytoplankton growth.In extremely turbid waters, the failure of atmospheric correction or oversaturation of the radiances measured by Aqua/MODIS will result in invalid satellite-derived chl a (He et al., 2012).
In this study, we limit the bloom statistics to the annual scale, because it is hard to determine the phenological statistics of phytoplankton blooms (e.g.initiation date and duration) in the eastern China seas with the current satellite dataset.The low PE in the ECS during wintertime also makes it difficult to calculate seasonal bloom statistics for each year.Although using the 8-day composite satellite chl a data may increase the proportion of valid data compared with daily data, it cannot resolve the short-term bloom events with durations shorter than 8 days.Therefore, the 8-day temporal resolution may systematically influence the bloom occurrence frequency.However, such influence is expected to be limited in the climatological or annual statistics in this study.

Climatological annual distribution of blooms
Figure 6b shows the annual climatology of F bloom derived from all 8-day composite chl a data from SeaWiFS (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).Blooms mostly occur in the Changjiang Estuary (CJE) and on the coasts of Zhejiang (CZJ), where the maximum F bloom is larger than 20 % (about 73 days per year).Meanwhile, along some coastal ar- eas of the NYS and BS, F bloom is also as high as 20 %, though the area with high F bloom is much smaller than that in the CJE and CZJ.In the middle of the southern Yellow Sea, F bloom is as high as 5 % for large blooms with chl a higher than 10 µg L −1 , which are mainly the spring blooms.In contrast to the coastal waters, F bloom in the middle and outer shelves of the ECS is less than 3 %.Almost no strong blooms occur in the oligotrophic Kuroshio region.In addition, no strong blooms occur in the extremely turbid waters of the HZB, CJM and SBS, which is further evidence for the suitability of taking 10 µg L −1 as the threshold (to avoid false blooms) for identifying a large bloom event in marginal seas with turbid waters.

Climatological seasonal distribution of blooms
Figure 7 shows the seasonal climatology of F bloom and PE derived from all 8-day composite chl a data from SeaWiFS (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).The spring, summer, autumn and winter are defined as April to June, July to September, October to December and January to March, respectively.Both F bloom and PE have significant seasonal variations, with most of the blooms occurring in spring and summer.In spring, the highest F bloom occurs along the CJE and CZJ.In some coastal areas of the BS and NYS, there are also large F bloom values in spring.In addition, spring phytoplankton blooms occur in the central YS with F bloom values up to 8 %.In summer, the region along the CZJ with high F bloom shrinks to the northern part.Meanwhile, the high F bloom area in the CJE extends northeastward along the CDW during summer.The spring blooms occurring in the central YS disappear in summer.Similar to the CJE, the west coast of the Korean Peninsula also has high F bloom values.F bloom reaches its highest in summer in the BS.In autumn, F bloom is much smaller than in spring and summer, and the high F bloom values are located in the NYS, CJE and CZJ.In addition, autumn blooms are found in the central YS.In winter, the bloom frequencies decrease further compared to those in autumn, with high F bloom values still occurring in the NYS.In addition, short-term blooms occur in winter in the middle shelf of the southern ECS.
Regarding the seasonal variations in PE (Fig. 7b), summer has the highest value, indicating the minimal cloud coverage year-round in the eastern China seas.The spatial distribution of PE is quite uniform, with PE greater than 75 % in summer, except in the extremely turbid waters in the HZB, CJM and SBS.Winter has the lowest PE, with values of less than 50 % in the ECS and a minimum of less than 25 %.The low PE may affect the statistics in winter, but with the combination of the 14 yr of data and two satellite datasets, the influence of low PE may not be significant, as evidenced by the low correlations between the spatial distributions of F bloom and PE in winter.

Interannual variation between 1998 and 2011
Figure 8 presents the annual PEs derived from SeaWiFS and Aqua/MODIS during the past 14 yr.The spatial distributions of annual PE for all years from both satellites are quite similar, with high values in the BS and YS, and low values in the ECS.The average annual PEs for the whole of the eastern China seas derived from SeaWiFS and Aqua/MODIS are highly consistent in both quantity and interannual variation (Fig. 9).Both SeaWiFS and Aqua/MODIS show that the annual PE averaged over the whole of the eastern China seas fluctuated between 55 and 65 % over the past 14 yr, with the minimum value in 2006.Although the annual PE val-ues from both datasets show slightly decreasing long-term trends, the rate of decrease is quite small, with −0.204 % per year for SeaWiFS (1998-2007) and −0.297 % per year for Aqua/MODIS (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011).These results show that the valid data for the 14 yr study period (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) are quite stable for the following analysis of interannual bloom variation.Figure 10 presents the annual F bloom derived from SeaW-iFS and Aqua/MODIS for the past 14 yr.Similar to the annual climatology, high F bloom occurs in the CJE and CZJ, but with significant interannual variation.Although we normalised and corrected the satellite-derived chl a, the bloom occurrence frequency derived from SeaWiFS is still systematically lower than that derived from Aqua/MODIS in the YS and BS.The difference may be caused by multiple factors, including the different sensor bands, atmospheric correction, chl a retrieval algorithms, etc., and further efforts should be taken to eliminate this difference in the future.
To quantitatively examine the interannual variation in F bloom , we selected two regions with a high frequency of bloom occurrence in the eastern China seas to obtain the annual bloom intensity index (BI) defined by Eq. ( 5), as shown in Fig. 11a.Region A located at the CJE corresponds to the area strongly influenced by the Changjiang River runoff.Region B covers the BS and YS, corresponding to the area weakly influenced by the Changjiang plume.These two regions cover almost all of the areas with high F bloom .Figure 11b and c show the annual BIs derived by SeaWiFS and Aqua/MODIS from 1998 to 2011 in these two regions.Although there are systematic differences between the annual BIs derived by SeaWiFS and Aqua/MODIS, the interannual variations are quite consistent.Region A in the CJE has significant interannual variation, but there is no significant increasing or decreasing trend over the 14 yr period, as revealed by Fig. 11b.However, in Region B, covering the YS and BS, there is a significant increase after 2002 and the BI more than doubles over the 14 yr (Fig. 11c).
We carried out a sensitivity test to examine the effect of threshold values (8 µg L −1 , 10 µg L −1 , 12 µg L −1 ) on the bloom results, as shown in Fig. 11b and c.In comparison with the BI derived using the 10 µg L −1 threshold, the mean differences in the BIs derived with 8 µg L −1 and 12 µg L −1 thresholds are 9.8 % and 6.4 % for Region A, and 22.7 % and 14.7 % for Region B. Although the BI varies slightly with the threshold values between 8 µg L −1 and 12 µg L −1 , the in-terannual variation and long-term trends are consistent.The results prove that a slight variation of our chosen 10 µg L −1 threshold for large bloom identification has little effect on the long-term trends.
We further examined the relationship between the longterm changes in the BI and the El Niño-Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO).Although it is difficult to assess the effects of ENSO and PDO on the blooms based on only 14 yr time series data, we believe it is still useful to analyse their relationships.We compared the interannual variations in the annual BI with the monthly Multivariate ENSO Index (MEI) and PDO index, as shown in Fig. 12.In Region A in the CJE, the interannual variations in the BI and MEI/PDO match relatively well, with high BI generally occurring in warm phases and low BI in cool phases, as revealed by Fig. 12a and c, indicating that the ENSO and PDO may modulate the bloom intensity in the CJE.In Region B, covering the YS and BS, there are less significant variations compared with those in the CJE, as revealed by Fig. 12b and d.To quantitatively analyse these relationships, we generated scatter plots to compare the 12-month averaged MEI and PDO indices with the annual BI for Region A, as shown in Fig. 12e and f, with 0, 3, 6 and 12 months delay, meaning that the 12-month averaged MEI and PDO indices are for the 0, 3, 6 and 12 months preceding the annual BI.The correlations between the annual BI and the 12-month averaged MEI and PDO indices are all positive, with the highest correlation coefficients for the 6-month delay.In addition, the MEI has higher correlation coefficients than the PDO index.It should be noted that we only present the relationship between bloom intensity and the ENSO/PDO indices over a relatively short timescale (14 yr).Longer time series data are needed to reveal the effects of climate change on the variations in phytoplankton blooms and community structures in the future, when more data become available.

Discussion
The necessary prerequisites for a phytoplankton bloom to start are the availability of inorganic nutrients and sufficient solar radiation (Fennel et al., 1999).The ECS and YS feature spring phytoplankton blooms, while the highest ocean primary production occurs in late summer and early autumn in the BS (Gong et al., 2003;Shi and Wang, 2012).Our results on the seasonal climatology of bloom occurrence frequency in Fig. 7a also reveal that phytoplankton blooms in the eastern China seas mainly occur in spring and summer.During spring and summer, the upper mixed-layer depths are usually less than the euphotic layer depths.In addition, the classic Sverdrup's critical depth is generally larger than the euphotic layer depth (corresponding to the compensation depth in the Sverdrup's critical depth model) (Sverdrup, 1953).Thus, light may not be a limiting factor for phytoplankton blooms in the eastern China seas, except in highly turbid coastal waters.The nutrient supply is expected to be the major factor controlling the interannual variation of phytoplankton blooms in the eastern China seas.Thus, in the following sections we focus on the discussion of the nutrient supply in two typical areas, Region A and Region B in Fig. 11a.

Nutrient supply in Region A
The Changjiang is the world's fifth largest river in terms of volume discharge, and the fourth largest in terms of sediment load (Milliman and Farnsworth, 2011).The huge discharge of the Changjiang transports abundant terrestrial material to the estuary and further to large plume areas (Bai  et al., 2013).On the one hand, the Changjiang inputs large amounts of nutrients into the CJE, benefiting phytoplankton blooms.The annual nutrient flux to the CJE from the Changjiang is markedly higher than that to other estuaries in China (Shen et al., 1992).The annual fluxes of total inorganic nitrogen, phosphate, silicate and nitrate are about 8.9 × 10 6 t, 1.4 × 10 4 t, 2.0 × 10 6 t and 6.4 × 10 6 t, respectively, which are an order of magnitude higher than the inputs from the Huanghe (Gao and Song, 2005).On the other hand, the annual sediment load of the Changjiang is about 4.9 × 10 8 t, with an average suspended sediment concentration of about 500 mg L −1 (Gao and Song, 2005).Therefore, as phytoplankton growth is light-limited in the water column in the turbid estuary, phytoplankton grows slowly relative to dilution and vertical mixing and cannot accumulate in the euphotic zone in the estuary in spite of high nutrient concentrations (Yin et al., 2004;Wang et al., 2010).Thus, phytoplankton blooms mostly occur in the plume area outside the turbidity maximum in regions rich with terrestrial nutrients and low water turbidity (Fig. 8) (Zhu et al., 2009).The strong phytoplankton blooms, combined with column stratification and the inflow of bottom water from the TWS, lead to the development of hypoxia off the CJE during summer (Wang et al., 2012).
As the Changjiang is a major source of nutrients for the inshore area of the ECS (Zhang, 1996;Chai et al., 2009), interannual variation in the Changjiang discharge may control bloom intensity in the CJE (Chen et al., 2007).  in the CJE were positively correlated with the Changjiang discharge with a time lag of 0-2 months, suggesting that the interannual variation in chl a is controlled by interannual variation in the Changjiang discharge.Yet, from Fig. 13b, we can see that there are anomalies in several years, especially in 2006, suggesting that other factors might affect the phytoplankton blooms besides river discharge, such as the direction and area of the plume extension, the stability of the water column, tidal cycles, wind forcing and the grazing rates of zooplankton (Cloern, 1996;Yin et al., 2004).At the same time, increases in both nitrogen and phosphorus and decreases in silicate from the outflow of the Changjiang may have disrupted the ecosystem in the ECS, although the exact effect is yet unknown (Zhou et al., 2008;Zhu et al., 2009).In addition, the construction of the Three Gorges Dam (TGD), the world's largest hydropower project, impounded from June 2003, reduced the sediment load to 147 Mt yr −1 in 2004, which was 35 % of the average for the previous five decades (Yang et al., 2006;Xu and Milliman, 2009).Meanwhile, the annual maximal discharge from the Changjiang covaries very well with the ENSO and PDO indices, as revealed by Fig. 12a and c.This is consistent with the notion that large Changjiang discharges generally occur in warm phases, whereas low discharges occur in cool phases (Jiang et al., 2006;Zhang et al., 2007).Therefore, the ENSO and PDO may modulate the Changjiang discharge, and the interannual variation in the bloom intensity is hence affected in the CJE.Our results also indicate that the blooms might be affected by ENSO and PDO events six months later (Fig. 12e and  f).Similarly, Volpe et al. (2012) found that the phytoplankton biomass anomalies associated with extreme atmospheric anomalies (such as the cold 1998-1999 winter and the summer 2003 heatwave) showed a significant correlation with a ∼ 5 month time lag in the Mediterranean Sea.
There has been no significant long-term increase or decrease during the past 14 yr in the CJE in spite of the www.biogeosciences.net/10/4721/2013/Biogeosciences, 10, 4721-4739, 2013 significant interannual variation, according to the long-term changes in the BI in Fig. 11b.Despite the large reduction in sediment since the construction of the TGD, there has been no corresponding significant decrease in water discharge (Chai et al., 2009;Xu and Milliman, 2009).The TGD may have limited effects on annual discharge and associated nutrients, as indicated by the annual mean and maximal discharge in Fig. 13b.Nutrients added by the Changjiang have increased dramatically since the 1960s, mainly due to the use of chemical fertilizers and industrial/domestic sewage discharges (Zhang, 1996;Gao and Wang, 2008;Zhou et al., 2008), and such nutrient supply would compensate the reduction in sediments due to the TGD.Furthermore, the significant decrease in suspended sediment supply by the Changjiang may have increased water transparency in the estuary, thus favouring phytoplankton blooms, and could also compensate the effect of nutrient reduction due to the slight decrease in the river discharge.

Nutrient supply in Region B
The BS and YS (Region B) make up a semi-enclosed marginal sea between mainland China and the Korean Peninsula, with its mouth opening into the ECS.Both SeaWiFS and Aqua/MODIS reveal a doubling of the bloom intensity during the past 14 yr in the YS and BS (Fig. 11c).The doubling of blooms is very important for biogeochemical processes and the marine carbon cycle because of the high background biomass in these regions (Fig. 4a and d).Using 10 yr (Sept 1997 to Oct 2006) of monthly mean chl a data observed by SeaWiFS, Yamaguchi et al. (2012) also found that the chl a during summer had gradually increased over 10 yr in the YS, indicating possible eutrophication.Therefore, a question is whether nutrients have significantly increased in these regions or not? Figure 14 shows the temporal variations in the surface layer nitrate and phosphate concentrations along two latitudinal sections in the southeastern YS (Fig. 1a) between 1998 and 2005.The maximal concentrations of nitrate and phosphate generally occur in December or February due to strong vertical mixing in the winter.As expected, both the nitrate and phosphate concentrations along the two sections increased significantly from 1998 to 2005.At the central YS represented by station N1, the maximum nitrate concentration increased over the same period from 6.84 to 14.61 µmol L −1 , and the maximal phosphate concentration increased from 0.41 to 1.01 µmol L −1 .Similarly, at station M1 located in the central YS, the maximum nitrate concentration increased from 5.58 to 11.17 µmol L −1 , and the maximal phosphate concentration increased from 0.58 to 1.13 µmol L −1 .The average rates of increase for all of these serial stations were 0.37 ± 0.08 µmol L −1 and 0.05 ± 0.01 µmol L −1 per year for the nitrate and phosphate concentrations, respectively.Therefore, the more than twofold increase in nutrients may have supported the doubling of the bloom intensity in the YS and BS during the past 14 yr.
Many factors could have caused the nutrient increase, such as the incursion of the Kuroshio, and riverine inputs.From the nutrient data along two sections of the Kuroshio and ECS shelf (Fig. 15), we can see that no significant increases or decreases in nitrate and phosphate concentrations have occurred in the surface layer over the past 14 yr, indicating that the incursion of Kuroshio might not be solely responsible for the twofold increase in nutrients in the YS.The riverine nutrients in the eastern China seas have increased  dramatically during the past decades mainly due to the use of chemical fertilizers and industrial/domestic sewage discharges (Zhang, 1996;Gao and Wang, 2008;Zhou et al., 2008;Qu and Kroeze, 2010).Moreover, the water flush time in the YS and BS is about 2.2 yr due to the semi-enclosed boundary, which is much longer than that in the ECS (Liu et al., 2003b).The dramatic increase in riverine nutrients and long residence time might be responsible for the systematic increase in nutrients in the YS and BS during the past 14 yr.Climate change and dust deposition may also affect the level of nutrients.Although the twofold increase in nutrients might dominate the doubling of bloom intensity in the YS and BS, the ENSO and PDO can also modulate the interannual variation of bloom intensity in this region, with bloom intensity X.He et al.: Seasonal and interannual variability of phytoplankton blooms increasing in warm phases and decreasing in cool phases, as revealed by Fig. 12b and d.However, the effects of ENSO and PDO are much smaller than the influences in the CJE.At the Scripps Pier in the Southern California Bight, Kim et al. (2009) found that the annual mean chl a exhibited a similar increasing trend, although the interannual variation in chl a did not correlate with the ENSO index.Using satellite wind data, Chen et al. (2012) found that annual wind speeds increased slightly during 2000-2011 in all of the eastern China seas.The increase in wind speed might enhance the vertical mixing and favour the transportation of nutrient-rich bottom water to the euphotic layer.Stramska (2005) argued that wind forcing was a main factor influencing the interannual variability of phytoplankton blooms in the open ocean.In addition, the annual dust input to the entire ECS is estimated to be 17 Mt, which appears to be a significant source of sediments and dissolved Fe compared with riverine discharge (Hsu et al., 2009).Tan et al. (2011) also found that chl a and primary production increased and eventually became a bloom 1-21 days after the occurrence of dust storms in the southern YS.It would be worth investigating the longterm trends of dust input when such data become available in the future.

Conclusions
Using the time series satellite ocean colour data provided by SeaWiFS and Aqua/MODIS, we studied the interannual variability in phytoplankton blooms over a 14 yr period from 1998 to 2011.In the first stage, we dealt with the issue of bloom identification in turbid water to ensure a proper dataset to investigate the bloom variation.Then, we examined the annual and seasonal climatological distributions of phytoplankton blooms in the eastern China seas and analysed their spatial and temporal variations.The findings revealed a doubling of the bloom intensity index in the YS and BS during the past 14 yr, which we attribute to the more than twofold increases in nitrate and phosphate concentrations.Despite significant interannual variation, there was no long-term increase or decrease in the bloom intensity in the CJE, which is mainly regulated by the Changjiang River discharge.
The results of the present study offer a foundation for investigating the effects of phytoplankton blooms on biogeochemical processes and the marine carbon cycle in the eastern China seas, such as the air-sea flux of carbon dioxide and hypoxic zones induced by blooms.In future studies, it would be worth investigating the sources of the twofold nutrient increase in the YS and the long-term variation in nutrients in the ECS.We have presented a preliminary report on the relationship between bloom intensity and ENSO/PDO indices, but more evidence and deeper analyses are needed in the future when more data become available.

Fig. 2 .
Fig. 2. Comparison of annual mean chl a retrieved by SeaWiFS and Aqua/MODIS during 2003-2007.The black dashed line is the 1 : 1 line, and the red solid line is the linear regression on the logarithmic scale.

Fig. 3 .
Fig. 3. Comparisons of matched chl a from the satellite and in situ data: (a) comparison for SeaWiFS, and (b) comparison for Aqua/MODIS.

Fig. 9 .
Fig. 9. Interannual variation in PE averaged for all of the eastern China seas derived from SeaWiFS and Aqua/MODIS.

Fig. 11 .
Fig. 11.The locations of the two selected regions in the eastern China Seas (a), and the interannual variations in the annual bloom intensity index (kg m −1 ) in Region A (b) and Region B (c) derived from SeaWiFS (SWF) and Aqua/MODIS (MOD) with the threshold values of 8 µg L −1 , 10 µg L −1 and 12 µg L −1 .
Figure 13a shows the monthly discharges of the Changjiang from 1998 to 2011.The discharges were measured at the Datong Gauge Station, which has no tidal influence and is 624 km upstream from the river mouth.The maximum discharge usually occurs in June or July, inducing the highest bloom frequency during summer in the CJE, as revealed by Fig. 7a.Over the past 14 yr, 1998 and 2010 had the largest annual mean discharges, while 2006 and 2011 had the lowest.The maximum monthly discharges for 1998 and 2010 were 63 426 and 61 310 m 3 s −1 , compared with 38 534.48 and 36 916.13 m 3 s −1 for 2006 and 2011, respectively.Generally, the annual BI covaries with the annual mean discharge and the maximal discharge, as revealed by Fig.13b, indicating that the Changjiang discharge is related to blooms in the CJE.Compared with the annual mean discharge, the maximal discharge has a larger correlation coefficient (R = 0.62) between the annual BI and discharge because most blooms occur in spring and summer with high discharge (Fig.13c).Using 10 yr (1998 to 2006) of SeaWiFS-derived chl a data,Yamaguchi et al. (2012) also found that summer chl a maxima between 1998 and 2006

Fig. 12 .
Fig. 12. Comparisons of the interannual variations among the bloom intensity index (BI in kg m −1 ), the annual maximal discharge of the Changjiang (m 3 s −1 ) and the MEI and PDO indices.(a) and (b) are the comparisons between the BI and MEI in regions A and B as shown in Fig. 11a.(c) and (d) are the comparisons between the BI and PDO indices in Regions A and B. (e) and (f) are the scatter plots between the 12-month averaged MEI and PDO indices with the BI for Region A with 0, 3, 6 and 12 months of delay, meaning that the 12-month averaged MEI and PDO indices are for the 0, 3, 6 and 12 months prior to the annual BI.

Fig. 13 .
Fig. 13.(a) The monthly averaged Changjiang discharge measured at the Datong Gauge Station from 1998 to 2011 (from http: //xxfb.hydroinfo.gov.cn).(b) The correlation between the annual bloom intensity index (kg m −1 ) of Region A and the annual mean discharge and maximal discharge from 1998 to 2011.(c) The scatter plots between the annual mean discharge and maximal discharge with the bloom intensity index for Region A.

Fig. 14 .
Fig. 14.Temporal variations in the nitrate and phosphate concentrations at time-serial stations in the southeastern YS from 1998 to 2005.The upper and lower rows are sections N1-N5 and M1-M3, respectively, as shown in Fig. 1a.The dashed line is the linear regression versus time (year) for each serial station.

Fig. 15 .
Fig. 15.Temporal variations in the nitrate and phosphate concentrations at time-serial stations in the ECS from 1998 to 2005.The upper and lower rows are sections S1-S5 and K1-K6, respectively, as shown in Fig. 1a.