3.1 Processing chain

Figure 2 presents a schematic diagram of the processing chain used to derive the R_rs and the associated coccolithophorid bloom map. The processor initial stages (Sect. 3.1.1–Sect. 3.1.4) are applied to each image in turn. The images are then aggregated into daily composites and monthly climatologies. Each stage of the processor is sequentially discussed below.

3.1.3 R_rs processor: atmospheric correction

The total TOA radiance measured by a satellite contains contributions from atmospheric scattering, reflections from the sea surface and the water-leaving radiance. The water-leaving component is of primary interest here and is typically much smaller than the atmospheric signal. Consequently, we must perform an atmospheric correction procedure to isolate the signal of interest.

In general terms, the TOA radiance can be written as (Franz et al., 2007)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{L}^{\mathrm{TOA}}=[L^{\mathrm{Rayl}}+L^{\mathrm{A}}+{\mathrm{td}}^{\mathrm{s}}\cdot L^{\mathrm{wcap}}+{\mathrm{td}}^{\mathrm{s}}\cdot L^{\mathrm{w}}]\cdot {\mathrm{tg}}^{\mathrm{0}}\cdot {\mathrm{tg}}^{\mathrm{s}}\cdot f_{\mathrm{p}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }\left[\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{1}}{\mathrm{nm}}^{-\mathrm{1}}\right],\end{array}$$

where L^TOA, L^Rayl, L^A, L^wcap and L^w refer to the radiance (L) at the top of atmosphere, due to Rayleigh scattering, due to aerosol scattering, due to whitecaps and due to the water-leaving components, respectively. The transmission coefficients for atmospheric gases (tg) and the associated atmospheric scaling factors (td) are superscripted according to the solar (0) and sensor (s) viewing directions. AVHRR is minimally sensitive to changes in polarisation (Zhao et al., 2004), and the polarisation factor, f_p, is set to unity.

P5.3x provides calibrated TOA bi-directional reflectance which has been normalised by $F_{\mathrm{0}}\cdot (R/R_{\mathrm{0}}{)}^{\mathrm{2}}\cdot \cos ({\mathit{\theta }}_{\mathrm{0}})$, where F₀ is the solar constant (Neckel and Labs, 1984), R∕R₀ the Earth–Sun distance ratio and θ₀ the solar zenith angle. Therefore, it is convenient to recast Eq. (1) in terms of reflectance. As remote-sensing reflectance (R_rs) is defined as the water-leaving radiance (L^w) divided by the downwelling irradiance (E_d), and E_d is proportional to the incoming solar radiance (as shown in Eq. 2), Eq. (1) can be re-written as Eq. (3).

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{d}}=F_{\mathrm{0}}\cdot (R/R_{\mathrm{0}}{)}^{\mathrm{2}}\cdot \cos ({\mathit{\theta }}_{\mathrm{0}})\cdot {\mathrm{td}}^{\mathrm{0}}[\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}],{\displaystyle }& {\displaystyle }{R}_n={\displaystyle \frac{\mathrm{1}}{{\mathrm{td}}_n^{\mathrm{0}}}}\cdot {\displaystyle \frac{\mathrm{1}}{F_{\mathrm{0}}\cdot (R/R_o{)}^{\mathrm{2}}\cdot \cos ({\mathit{\theta }}_{\mathrm{0}})}}\\ \text{(3)}& {\displaystyle }& {\displaystyle }[{\displaystyle \frac{L_n^{\mathrm{TOA}}}{{\mathrm{tg}}_n^{\mathrm{s}}\cdot {\mathrm{tg}}_n^{\mathrm{0}}\cdot {\mathrm{td}}_n^{\mathrm{s}}}}-{\displaystyle \frac{L_n^{\mathrm{Rayl}}}{{\mathrm{td}}_n^{\mathrm{s}}}}-L_n^{\mathrm{wcap}}-L_n^{\mathrm{A}}],\end{array}$$

where R_n is the corrected reflectance for the given channel, denoted by the n subscript. As the first bracketed term in Eq. (3) is proportional to the TOA reflectance provided by the P5.3x data set, Eq. (3) can be re-written as

$$\begin{array}{}\text{(4)}& {\displaystyle }{R}_n={\displaystyle \frac{\mathrm{1}}{{\mathrm{td}}_n^{\mathrm{0}}}}\cdot [{\displaystyle \frac{R_n^{\mathrm{TOA}}}{{\mathrm{tg}}_n^{\mathrm{s}}\cdot {\mathrm{tg}}_n^{\mathrm{0}}\cdot {\mathrm{td}}_n^{\mathrm{s}}}}-{\displaystyle \frac{R_n^{\mathrm{Rayl}}}{{\mathrm{td}}_n^{\mathrm{s}}}}-R_n^{\mathrm{wcap}}-R_n^{\mathrm{A}}],\end{array}$$

where $R_n^{\mathrm{Rayl}}$, $R_n^{\mathrm{wcap}}$ and $R_n^{\mathrm{A}}$ are the channel n Rayleigh, whitecap and aerosol reflectance terms, respectively. Each of these correction terms are now discussed in turn.

Figure 2Schematic diagram showing the different stages of the remote-sensing reflectance (R_rs) processing chain. The blue-shaded region generates the unfiltered R_rs product; the red-shaded region subsequently generates the filtered R_rs product. PML, NOAA and ECMWF refer to Plymouth Marine Laboratory, the National Oceanic and Atmospheric Administration and the European Centre for Medium-Range Weather Forecasts, respectively.

Download

Rayleigh scattering is a function of wavelength, and satellite and solar-viewing angles. For each wavelength, the Rayleigh reflectance, $R_n^{\mathrm{Rayl}}$, is bilinearly interpolated from a look-up table of Rayleigh radiance components, discretised by 2^∘ solar and satellite zenith angles and normalised by the extraterrestrial solar irradiance (Neckel and Labs, 1984). The look-up table was calculated using values of 0.057 and 0.02 for the Rayleigh scattering optical depth (τ) for the AVHRR visible and NIR channels, respectively (Elterman, 1970). Each entry in the table contains three sets of Rayleigh radiance components, $A_{\mathrm{0},\mathrm{1},\mathrm{2}}$, and total Rayleigh reflectance is calculated using

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{R}^{\mathrm{rayl}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathit{\pi }}{\cos \left({\mathit{\theta }}_{\mathrm{0}}\right)}}\cdot \left[A_{\mathrm{0}}\right(\mathit{\lambda })+A_{\mathrm{1}}(\mathit{\lambda })\cdot \cos (\mathit{\delta }\mathit{\varphi })+A_{\mathrm{2}}(\mathit{\lambda })\\ \text{(5)}& {\displaystyle }& {\displaystyle }\cdot \cos \left(\mathrm{2}\mathit{\delta }\mathit{\varphi }\right)],\end{array}$$

where A₀(λ), A₁(λ) and A₂(λ) are the linearly interpolated Rayleigh radiances for a given wavelength, solar zenith angle and relative azimuth ($\mathit{\delta }\mathit{\varphi }=[\mathrm{180}-({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\left)\right]$), where θ_s is the satellite zenith angle.

Reflectance due to whitecaps is a function of both wind speed and wavelength. Here, this is calculated according to the method described by Koepke (1984). The whitecap reflectance is calculated using

$$\begin{array}{}\text{(6)}& {\displaystyle }{R}^{\mathrm{wcap}}\left(\mathit{\lambda }\right)=\mathrm{2.95}\times {\mathrm{10}}^{-\mathrm{6}}\cdot U_{\mathrm{10}}^{\mathrm{3.52}}\cdot R_{\mathrm{ef}}\left(\mathit{\lambda }\right),\end{array}$$

where λ is wavelength in microns; U₁₀ is the 10 m wind speed, as provided by the ingested ERA interim fields; and R(λ)_ef is the wavelength-dependent effective reflectance. The appropriate effective reflectance value is interpolated from a look-up table derived from Koepke (1984), generalised to include NIR wavelengths.

To correct for aerosol effects, we adopt the approach used by Smyth et al. (2004) and Uz et al. (2013). The assumptions here are twofold: firstly that the aerosol reflectance for channel 1 and channel 2 is equal (Stumpf and Pennock, 1989) and secondly that R_rs in channel 2 is zero. Using these assumptions, the reflectances in the visible and NIR channels can be used to isolate and remove the aerosol signal from the reflectance calculation, using Eq. (7):

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{r}{\displaystyle }{R}_{{\mathrm{rs}}_{\mathrm{1}}}-R_{{\mathrm{rs}}_{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\cdot {\displaystyle \frac{(R_{\mathrm{1}}-R_{\mathrm{2}})}{\exp (-\mathrm{0.5}\cdot {\mathit{\tau }}_{\mathrm{0}}\cdot \mathrm{pl})}}\\ {\displaystyle }{R}_{{\mathrm{rs}}_{\mathrm{2}}}\to \mathrm{0}\left[{\mathrm{sr}}^{-\mathrm{1}}\right],\end{array}\end{array}$$

where pl is the atmospheric path length ($\mathrm{1}/\cos {\mathit{\theta }}_{\mathrm{0}}+\mathrm{1}/\cos {\mathit{\theta }}_{\mathrm{s}}$) and τ₀ is the Rayleigh optical depth for channel 1 for a path length of unity. R₁ and R₂ are the respective channel 1 and 2 TOA reflectances (R^TOA), corrected for Rayleigh scattering (R^Rayl), whitecaps (R^wcap) and atmospheric transmission (as given by Eq. 4). R_rs values are not retained where the Rayleigh reflectance calculation fails.

Ozone and water vapour absorption values for AVHRR channels 1 and 2 are provided by Liang (2005) and Tanre et al. (1992), and implemented as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathrm{tg}}^{\mathrm{s}}\left(\mathit{\lambda }\right)={\mathrm{tg}}_{{\mathrm{O}}_{\mathrm{3}}}^{\mathrm{s}}\left(\mathit{\lambda }\right)\cdot {\mathrm{tg}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}^{\mathrm{s}}\left(\mathit{\lambda }\right)\cdot {\mathrm{tg}}_{{\mathrm{NO}}_{\mathrm{2}}}^{\mathrm{s}}\left(\mathit{\lambda }\right)\cdot {\mathrm{tg}}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{s}}\left(\mathit{\lambda }\right),\end{array}$$

where the ozone and water vapour transmission components along the sensor viewing path length are calculated according to equations of the form

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{tg}}_{{\mathrm{O}}_{\mathrm{3}}}^{\mathrm{s}}\left(\mathit{\lambda }\right)=\mathrm{1.0}/(\mathrm{1.0}+B_{\mathrm{0}}(\mathit{\lambda })\cdot (\mathrm{1.0}/\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)\cdot \left[{\mathrm{O}}_{\mathrm{3}}\right]\left(\mathit{\lambda }\right){)}^{B_{\mathrm{1}}\left(\mathit{\lambda }\right)}).\end{array}$$

B₀(λ) and B₁(λ) are the wavelength-dependent absorption coefficients for each channel and gas. Concentration values for ozone [O₃], and for water vapour [H₂O], are interpolated from daily ERA-Interim fields. Equations (8) and (9) are similarly constructed for the solar-viewing path length. NO₂ and CO₂ scaling factors are set to one as their absorption is assumed to be negligible (see Liang, 2005, Fig. 2.10). Atmospheric scaling factors for each viewing path are calculated, using Eq. (10), where τ_R(λ) is the Rayleigh optical scattering depth.

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathrm{td}}^{\mathrm{s}}\left(\mathit{\lambda }\right)=\exp (-\mathrm{0.5}\cdot {\mathit{\tau }}_R(\mathit{\lambda })\cdot (\mathrm{1.0}/\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)\left)\right)\end{array}$$

Sunglint is explicitly flagged in and removed from the Px5.3 data sets, and no further correction for sunglint is applied.

3.2 Filtering, masking and identifying blooms

In previous ocean-colour-based analyses, coccolithophorid bloom maps are produced as the binary classified output of a supervised multi-spectral algorithm (e.g. Iglesias-Rodriguez et al., 2002; Brown and Yoder, 1994). In this work, the availability of only one visible channel necessitates an alternative method, and the bloom map is instead produced through temporal filtering of the R_rs product, followed by selective masking to subsequently remove false positives.

Temporal filtering of the R_rs product is performed through a comparison of each daily composite to the relevant monthly mean climatological R_rs field (produced by the total aggregator stage). R_rs signals are only classified as blooms where the per-pixel R_rs value is greater than 2 monthly standard deviations above the corresponding monthly mean value. The standard deviation in this case is calculated from the monthly mean products across the entire archive. Pixels that do not match this criterion are assumed to contribute to the background, rather than bloom signal, and are therefore set to zero. The filtered bloom product, written into the daily netCDF4 file as filtered_remote_sensing_product, is then subject to further quality controls in the masking stage, as described below.

High R_rs values, while potentially indicative of coccolithophore blooms, can also occur in regions that are subject to high concentrations of suspended sediment (e.g. estuaries), or where shallow bathymetry and clear water coincide (e.g. shelf regions in oligotrophic areas). To remove these, and other false positives, the final bloom product is derived from the filtered bloom product by subjecting the latter to a number of screening processes, as detailed below.

Firstly, to remove the effects of land contamination, the bloom map is set to zero in all points within three pixels of the land mask. Secondly, following Iglesias-Rodriguez et al. (2002), the bloom map is set to zero in areas between 47^∘ N and 47^∘ S where the bathymetry is shallower than 100 m. This removes false positives associated with the sea floor, an effect that is most noticeable in the Caribbean and Arafura seas. Thirdly, whilst flagged sea ice has been explicitly removed from the Px5.3 data (see Sect. 3.1.1), this does not comprehensively remove ice effects. As a result of missed flagging, and of glacial (Broerse et al., 2003) and river run-off, sporadic high R_rs values that are not indicative of blooms still occur at high latitudes. To correct for this, bloom map pixels are set to zero where R_rs≥0.05 (a value far above that which we would expect in water types associated with coccolithophorid blooms). Furthermore, the R_rs product is screened using sea surface temperature (SST) data obtained from contemporaneous ERA-Interim fields, and R_rs is set to zero in pixels where SST <0 ^∘C in the Northern Hemisphere, a value at which the coccolithophorid growth rate drops to near zero, even for cold-water strains such as E. huxleyi (Buitenhuis et al., 2008). Finally, bloom map pixels are set to zero where the total aggregated mean value R_rs is greater than 0.0005, removing the effects of consistent river outflows (e.g. the Amazon and around the Yellow Sea).

The final product suite is annotated with relevant metadata to ensure CF1.8 compliance, completing the processing. The contents of the data file are described fully in the following section.

5.1 Regional comparisons

Figure 3 shows a comparison between four coccolithophorid blooms detected by ocean colour sensors and the corresponding blooms in the filtered_remote_sensing_reflectance product. In all four cases where bright blooms are detected in the ocean colour sensor observations (Fig. 3a SeaWiFS; Fig. 3c, MERIS; Fig. 3e, g, MODIS) there are spatially corresponding bright patches in the AVHRR imagery. In the MERIS and MODIS cases the AVHRR imagery is from the same day (i.e. on a single overpass basis). In the SeaWiFS case, a 3-day AVHRR composite mean is used, due to differences in cloud cover at the various acquisition times, lowering the intensity of the visible bloom but preserving the spatial coverage of the scene.

Figure 3Examples of ocean-colour-derived red–green–blue (RGB) images of coccolithophore blooms matched to their filtered bloom product counterparts. Panels are as follows. (a, c, e, g) Level 2 RGB images for the North Sea and English Channel (SeaWIFS; 30 July 1999), North Atlantic and Irish Sea (MERIS; 23 May 2010), Barents Sea (MODIS; 17 August 2011) and Bering Sea (MODIS; 4 September 2014). (b, d, f, h) Matching, contemporaneous filtered bloom product composite for each location and date. Dark grey indicates land, and light grey indicates cloud, throughout. For bloom products, dark blue indicates that no bloom is present; lighter cyan colours indicate that a bloom is present.

Download

There is also evidence from in situ data in the English Channel case (Fig. 3a and b) that this was indeed a bloom of Emiliania huxleyi from cell counts and in-water radiometry (Smyth et al., 2002). The northern North Sea feature centred on 56^∘ N, 1.5^∘ E in Fig. 3a is possibly the remnants of a bloom which was the subject of an intensive field campaign during June 1999 (Burkill et al., 2002). Similarly blooms have been documented in the literature in the Barents Sea (Fig. 3e and f) which are comparable in timing and extent with some limited in situ samples (Smyth et al., 2004); Merico et al. (2003) report on blooms in the Bering Sea of similar timing and extent to those shown in Fig. 3g and h.

It is notable that the cloud (and ice) masking algorithms for the ocean colour and AVHRR sensors are in close agreement for the individual scenes shown in Fig. 3c and d, e and f, and g and h. The discrepancy in the cloud flagging for the SeaWiFS case occurs as a result of the 3-day compositing.

5.2 Global signals

Figure 4 shows the global spatial coverage of the data set, presented as decadal means of the filtered R_rs bloom product for four decadal periods. Comparing Fig. 4a with the CZCS era (1978–1986) coccolithophorid bloom map produced by Brown and Yoder (1994) suggests that the bloom signals in the North Atlantic, North Sea and Norwegian Sea and on the Argentinian coast are well captured. The presence of blooms in the Black Sea and sporadically throughout the Mediterranean is also consistent between the two. However, due to the coastal masking, the high signals around Newfoundland are not captured here. In addition, no signals are detected in the Arafura Sea or within the Indonesian archipelago, as these areas have been specifically masked due to shallow bathymetry. The signals along the coast of Greenland and in the Southern Ocean are stronger than anticipated, likely due to some remaining influence of ice and glacial river outflow.

Figure 4Mean values for the filtered remote-sensing reflectance (R_rs) bloom product by decade for the (a) 1980 to 1989, (b) 1990–1999, (c) 2000–2009 and (d) abbreviated 2010–2017 period.

Download

Mean values for the 1990–1999 period suggest the presence of coccolithophorid blooms in the North Atlantic, Norwegian Sea, Baltic Sea, Bering Sea and Southern Ocean, which is consistent with previous findings of Iglesias-Rodriguez et al. (2002). Similarly, the R_rs-based bloom product correctly attributes signals to the Benguela upwelling and in the western North Pacific Ocean. While they do not cover identical periods to the record shown here, increased R_rs in the Black, Norwegian and Baltic seas between the Brown and Yoder (1994) and Iglesias-Rodriguez et al. (2002) studies is well replicated. A reduction of R_rs along the coast of Argentina also appears to be appropriately captured (Brown and Yoder, 1994).

The change in bloom patterns between the CZCS era (Fig. 4a) and the SeaWiFS (Fig. 4b and c) mirrors the findings of Winter et al. (2013), with an intensification of bloom expression and increase in coverage in the North Atlantic. Following this, despite a substantial decrease in mid-latitude bloom expression in Fig. 4d, there is a marked increase in the intensity of blooms in the Barents Sea, consistent with Neukermans et al. (2018).

6.1 Radiometric sensitivity and grid resolution

The differences between the bloom extent and intensity observed by the ocean colour and AVHRR sensors in Fig. 3 can, in part, be attributed to a combination of lower spatial resolution and radiometric sensitivity of the AVHRR sensor. The typical spatial resolution of ocean colour sensors is between 300 m and 1.1 km, whereas the Px5.3 product is 0.1^∘. This equates to 2 orders of magnitude of coarsening and will have a particularly pronounced effect where strong spatial radiometric heterogeneities exist within blooms (see e.g. Smyth et al., 2002), resulting in lower overall bloom reflectances. This coupled with the lower radiometric sensitivity of AVHRR (3 %) will exacerbate the overall “dimming” of the bloom. It is unlikely that blooms of non-calcifying phytoplankton with low backscattering will be detectable using this approach.

6.2 Limits of detection

Given the expected “bloom dimming” (as discussed in Sect. 6.1), it is necessary to establish the concentration at which coccolithophorid blooms can be detected using this approach. Comparing Fig. 3a, c, e and g, we can see that only the densest parts of the bloom are detected by AVHRR. This will be in the parts of the bloom undergoing the closing or senescent stages where the reflectance signal will be dominated by coccoliths rather than live coccolithophore cells.

Holligan et al. (1993) in their detailed study of the biogeochemistry of coccolithophore blooms in the North Atlantic calculated that the AVHRR threshold of detectability for coccoliths is around 2×10⁴ liths mL⁻¹. A similar type of calculation using the equations found in Gordon et al. (2001) and accounting for shifts in wavelength (0.546 to 0.630 µm) yields a threshold of approximately 10⁵ liths mL⁻¹. The bloom shown in Fig. 3a is well described by Gordon et al. (2001). When comparing the in situ data from Gordon et al. (2001) Figs. 1 and 3a with the bloom pattern here, it is apparent that only the near-shore transect point is bright enough to be detected by AVHRR. However, for the same bloom Smyth et al. (2002) showed that the number of free liths enumerated via microscopy using manual counting of samples fixed in buffered formalin, as used in Gordon et al. (2001), yielded a factor of 2 greater than that using flow cytometry. This they attributed to the fixation process causing the coccolithophores to shed even more of their liths and was proved by comparing fresh and preserved samples of cultured E. huxleyi using flow cytometry. Additionally Smyth et al. (2002) achieved closer modelled optical closure when using the flow cytometric rather than the manual counts. This would therefore reduce the threshold limit described above to around 5×10⁴ liths mL⁻¹. As the Px5.3 data are cross-calibrated across all sensors, we do not expect this sensitivity threshold to be platform-specific.

6.3 Atmospheric correction

The availability of only a single visible channel greatly reduces the ability of AVHRR to spectrally determine in-water constituents. Consequently, coccolithophorid blooms are predominantly identified through the removal of likely false positives (e.g. river plumes, coastal influences, bathymetry and sea ice). Whilst coccolithophorid blooms are known to occur across the eastern (Malinverno et al., 2003), north-western (Oviedo et al., 2015) and western Mediterranean (Ignatiades et al., 2009), as shown in the bloom map of Iglesias-Rodriguez et al. (2002), they do not do so to the extent that they are expressed in the R_rs product derived here (Fig. 4b). We partially attribute the high R_rs values to the sporadic presence of Saharan dust. This requires specific atmospheric correction methods (Moulin et al., 2001) that are not implemented here and would require ancillary dust information that is of limited availability on the timescale covered by this data set.