2.1 Data

2.2 Ocean Color Type Classification

Monthly mean reflectance values were used for defining OC types due to their improved spatiotemporal coverage and more manageable number of samples. All data available for the Northeast Atlantic Shelves area were used in the analysis (n = 1,244,874). OC types were defined via clustering, similar to the approach of Mélin and Vantrepotte (2015). Reflectance spectra were first normalized by dividing reflectance values by the area under the spectrum, resulting in an area of 1.0. Area under the spectrum was calculated with trapezoidal integration (R package caTools, Tuszynski, 2020). Clustering was conducted with k-means, using the default Hartigan-Wong algorithm (Hartigan & Wong, 1979) and four restarts to increase the stability of the cluster centers. The total number of clusters was selected according to the lowest dynamic validity index (DVI) proposed by Shen et al. (2005), which provides an objective metric for defining the optimal number of clusters. The DVI combines both a measure of compactness (minimize intra-cluster distances) with a measure of separateness (maximize between or inter-cluster distances), which allows for the identification of homogenous, well-separated clusters.

The cluster centers (mean r_n spectra) were then used to classify the OC type of two higher temporal resolution OC data sets (daily and 8-day mean reflectance). These data were used for match-ups with associated physical and biological covariates in order to further characterize the OC types (see next section). The higher temporal resolution reflectance spectra were normalized, as done with the monthly mean spectra, and prediction of OC type was based on the nearest neighbor Euclidean distance to the cluster centers (R package FNN, Beygelzimer et al., 2019).

2.3 Ocean Color Type Characterization

OC types were characterized in terms of associated covariates, both physical and biological in nature. In order to address whether differing levels of temporal aggregation affected the relationships, two temporal resolutions were used for all variables (daily and 8-day means). Since the measurement of OC data is hindered by low light conditions during winter months, only the months of March-October were used in the analyses. The predicted OC type data sets were matched to all remotely-sensed data of similar temporal aggregation, which included OC product parameters (CHL1, CHL2, POC, PIC, and TSM), and the abiotic reanalysis parameters (SST and SSS). The differences in the associated covariates among OC types were evaluated visually via boxplots displaying quantile distributions.

Canonical Correspondence Analysis (CCA) (ter Braak, 1986, 1987) was used to elucidate relationships between the zooplankton assemblages (i.e. constrained data) and environmental variables (i.e. constraining data), including OC type. CCA assumes a unimodal (e.g. Gaussian) relationship between environmental covariates and zooplankton taxa, which implies an optimal set of environmental conditions associated with the maximum abundance (or biomass) of a given taxa. CPR samples were matched to environmental covariates based on the overlap of CPR sample coordinates to cells of the 25 km resolution grid used by other data. Rare zooplankton taxa, found in < 1% of samples, were removed from the community analyses in order to focus on more prevalent groups. CPR sample abundances were log(x+1)-transformed in order to down-weight the most abundant taxonomic groups. Environmental covariates included the continuous variables log (CHL1), SST, SSS, and the categorical variables of month and OC type. The categorical factor for the month of the sample was also included to account for seasonal changes in abundance related to life history strategies, such as overwintering. Although CHL2 was used as a covariate in the characterization of OC types, it was not included as a constraining variable due to its more restricted (but subjective) association with coastal waters. Significance tests of the overall model fit, CCA axes, and marginal term effects were performed via an ANOVA-like permutation test (Legendre et al., 2011; Legendre & Legendre, 2012). The test statistic is “pseudo-F”, which is the ratio of constrained and unconstrained total inertia (Chi-square), divided by the degrees of freedom. Inertia is a measure of weighted variance among species, focusing on gradients in the overall community composition. The ratio of inertia explained by the CCA has been shown to be a poor indicator of the explained compositional variation and it is thus recommended to focus on the relative amounts of inertia explained by CCA axes and constraining variables when interpreting the model results (Økland, 1999). In the summary tables for each permutation test, relative importance is inferred from the reported Chi-square and F-values.

Similarity percentage analysis (SIMPER) (Clarke, 1993) was used to identify zooplankton taxa that differed between OC types. SIMPER performs pairwise comparisons of groups (here OC types) and quantifies the contribution of each taxa to the average between-group Bray-Curtis dissimilarity. Although abundance data were log-transformed as with the CCA, highly abundant and variable species can have high contributions to dissimilarity, and thus permutation tests are performed to determine whether these differences are significant. The R package vegan (Oksanen et al., 2019) was used for all community analyses.

2.4 Drivers of Ocean Color Province Variability

Additional focus was placed on understanding the historical dynamics of OC types in the North Sea. The fraction coverage of total North Sea area (in km²) by OC type was calculated from the monthly mean OC data. In order to summarize long-term and seasonal temporal patterns of the resulting coverage time series, a generalized additive model (GAM) was fit using the R package mgcv (Wood, 2011, 2020), which allows for the incorporation of smooth functions to describe nonlinear effects. The model assumed a quasi-binomial distribution with logit-link function, considered appropriate for proportional data, and contained two explanatory terms; 1) Year – Thin-plate regression spline applied to numeric date values for describing interannual patterns; 2) Month – Cyclic cubic regression spline applied to numeric month values for describing seasonal patterns.

For each OC type, an Empirical Orthogonal Function (EOF) analysis (also known as Principal Component Analysis) was conducted to extract the dominant spatiotemporal modes of variability; specifically, vectors describing temporal (principal components, PCs) and corresponding spatial (EOFs) modes of variability. The presence/absence (i.e. binomial) data of each OC type were logit-transformed to approximate a normal distribution. Although missing values were few in non-winter monthly-averaged fields (∼1%), gaps were filled via Data Interpolating Empirical Orthogonal Functions (DINEOF) (Beckers & Rixen, 2003) using the R package sinkr (Taylor, 2020). Following interpolation, monthly anomalies were derived by substracting the long-term monthly mean of each spatial grid, which removes the seasonal signal for analyses focused on interannual variability. An EOF was conducted on both the original and anomaly data fields; the original data were used to describe the dominant patterns in annual spatiotemporal variability, while the anomaly field was decomposed as a precursor to the Redundancy Analysis described in the next paragraph.

Redundancy analysis (RDA) (von Storch & Zwiers, 1999; Tyler, 1982; Rao, 1964) was used to elucidate physical drivers of OC type spatiotemporal dynamics on interannual scales. In particular, a field of depth-averaged current (CUR) anomalies (i.e. constraining data) was used to explain the field of each OC type's occurrence anomalies (i.e. unconstrained data). In contrast to the unimodal relationships described by CCA modes, RDA describes linear relationships between the constraining and unconstrained data sets. This assumption was deemed appropriate for evaluating the degree to which the spatiotemporal dynamics of OC provinces were related to water transport.

Monthly anomalies were derived by subtracting the long-term monthly mean of each grid for all fields. Similar to Kauker et al. (2008), we first conducted an EOF on the CUR and OC type fields in order to identify a truncated number of temporal PCs that cumulatively explain > 90% of field variance. These PCs were then used as input data to the RDA rather than the original full data fields, which has several advantages over directly using the matrices of the data fields; principally, reduced computational effort, increased stability of solution, and avoidance of problems of collinearity in both the predictor and response data sets. Unscaled PCs (i.e. with units included) were used so that their importance in terms of explained variance of original field was maintained. Temporal patterns are a direct output of the RDA, while the associated spatial components are derived by projecting the unitless EOF patterns onto the RDA loadings associated with each PC. Significance of RDA patterns were tested via an ANOVA-like permutation test (Legendre et al., 2011; Legendre & Legendre, 2012). The R package vegan (Oksanen et al., 2019) was used for RDA and permutation tests. The test statistic is “pseudo-F”, which is the ratio of constrained and unconstrained total variance, divided by the degrees of freedom. Contrary to CCA, RDA uses linear relationships to best describe variance of the unconstrained data, and thus the fraction of total explained variance is an appropriate measure of the model overall performance.

Links between other environmental drivers (NAO and SPG) and the dominant modes of OC type spatiotemporal dynamics were explored using correlation analysis of annual mean values. Mean annual PC1 loadings derived from the EOF of monthly anomalies of occurrence were used for the OC types. Yearly indices of OC types were correlated to the NAO and SPG at lag 0 and lag 1 to address possible advective delays between conditions of the North Atlantic and the North Sea. Spearman's rank correlation coefficient (ρ) was used, which is a nonparametric metric capable of describing nonlinear, but monotonic, relationships.

3.1 Ocean Color Type Classification

Clustering of the monthly mean spectra for the Northeast Atlantic Shelves area showed lowest DVI for 4 OC types (Figure 3). The variance contained within the clusters was > 80% of total variance, after which additional clusters showed only marginal improvements. Cluster mean spectra (values provided in Table S1 in Supporting Information S1) have been ordered according to decreasing normalized reflectance at 412 nm, which is a wavelength particularly associated with Chla absorbance that aided later interpretation of OC types. For the remainder of the paper, these numerical identifiers (OC1 – OC4) will be used for the OC types. An example of the resulting OC type assignment to spectra from a single year (2013) is presented in Figure 4. The spatial pattern clearly shows a more Atlantic distribution for OC1 and OC2, while OC3 and OC4 are more coastally-distributed. OC4 is particularly associated with a littoral distribution.

In general, the predicted classification of OC type was consistent between 8-day and daily match-ups (Table 1). Exact overlapping classification occurred in 76% of cases (Table 1 diagonal), while overlap with neighboring classifications (in terms of OC type order) occurred in 23% of cases (Table 1 sub- and superdiagonals). Overlap to more distant classifications occurred in 1% of cases. The high degree of overlap indicates that OC types are relatively stable over both time scales for a given grid. One exception was OC4, which was classified more frequently in daily than in 8-day mean reflectance data. The predictions of OC4 using daily mean data (n = 47,647) were split nearly equally between predictions of OC3 (n = 22,780) and OC4 (n = 24,867) using 8-day mean data, possibly indicating less persistence of this type over time.

3.2 Ocean Color Type Characterization

3.2.1 Bio-Optical and Physical Variables

23,571 of 42,563 CPR samples (56%) could be matched to the 8-day mean ocean color data (Figure 5), while daily mean data were matched in only 9,255 cases (22%) (Figure S2 in Supporting Information S1). Unmatched samples are due to missing remote OC estimates resulting from cloud coverage or low-light winter conditions. The number of matched samples was largely consistent across the considered months (March-October, Figure 5, bottom-right), although counts by year show lower values at the beginning of the series (Figure 5, top-right) due to the existence of a single ocean color sensor prior to 2003 (i.e. SeaWiFS, see Figure S1 in Supporting Information S1).

The OC types of the CPR match-up samples are associated with trends in a number of other OC products. Only 8-day mean match-ups are shown here (Figure 6), although very similar patterns were observed in the daily mean match-ups (Figure S3 in Supporting Information S1). In particular, OC types are ordered (OC1-OC4) along a gradient of increasing Chla (CHL1 and CHL2), POC, TSM, and along a decreasing gradient of SSS. As seen in Figure 4, these associations reflect spatial differences between the coastal and offshore environments. PIC increased with OC type number but shows less distinguishable distributions. SST showed no clear differences or trends among OC types.

3.2.2 Zooplankton Community

Following the filtering of CPR data to include most relevant taxa (i.e. taxa found in >1% of samples) and with match-ups to environmental parameters, 19,884 (17 taxa) and 7,821 (16 taxa) samples were included in the community analysis for 8-day and daily mean match-up data sets, respectively. Both data sets resulted in very similar results, and thus only 8-day match-up analyses are shown here. Daily-mean match-up results are available in the Supplementary Information (Figure S4 and Tables S4–S7 in Supporting Information S1).

The CCA showed that the environmental covariates significantly explained patterns in the zooplankton community (Table S2 in Supporting Information S1). The high residual inertia left unexplained by the CCA is likely due to several factors relating to a low signal to noise ratio in the zooplankton data, which may include: the large spatiotemporal coverage, the patchy nature of zooplankton distributions, and the less than optimal surface transect sampling of the CPR data set (as opposed to sampling through the water column). Despite these less than optimal characteristics, the sheer size of the data set allowed for the identification of significant community associations with the covariates.

All environmental terms were determined to be significant following a test of marginal effects (Table 2), with SST associated with the highest F-value. Variables that are more clearly seasonal in nature across the entire sampled area, e.g. SST and Month, explained a higher degree of community variance than other parameters. OC type, SSS and CHL1 all explained similar degrees of variation, however the F-value of OC type is lower due to the higher associated degrees of freedom. The biplot of the leading CCA axes shows two main directions of variability related to environmental variables (Figure 7); one associated with seasonal patterns, including Month and SST, and an orthogonal one, associated with SSS, OC type, and CHL1. Nine CCA axes were deemed significant at the p ≤ 0.05 level, but a large majority of the explained inertia is contained within in the first two (Table S3 in Supporting Information S1). The CCA coordinates of zooplankton taxa (i.e. location of labels) correspond to the environmental conditions where higher abundances occur. As the biplot only displays the first two axes, direct interpretation can be misleading; yet the results of the SIMPER analysis (see below) general support the affinity of taxa to the most closely located OC type in the CCA biplot.

Table 2. Permutation Test of CCA Term Marginal Effects (n = 999) Relating Environmental Variables (8-Day Mean Match-Ups) to Changes in Zooplankton Community

	Df	Chi-square	Pseudo-F	p-value*
OC type	3	0.009	16.651	0.001
Month	7	0.074	60.262	0.001
SST	1	0.044	251.891	0.001
SSS	1	0.018	102.948	0.001
log(CHL1)	1	0.012	69.137	0.001
Residual	19870	3.481

• Note. Chi-square is a measure of the community inertia explained by each environmental variable. Explanatory variables OC type and Month are factors and others are continuous. *All variables showed significant effects at the p ≤ 0.001 level.

A summary of the SIMPER analysis can be seen in Table 3, where only the top species contributing to 95% dissimilarity between groups, and for which permutation tests show significant differences (p < 0.05), are presented. Many taxa show clear abundance gradients correlating to the OC types (ID numbers referred to in Figure 7 and Table 3 are provided after taxa names for reference). Metridia lucens (3), Calanus finmarchicus (1), Oithona spp (14), and Clausocalanus spp (13) are more associated with Atlantic waters of OC1 and OC2, as reflected by the mean abundances by OC type in the SIMPER analysis, and their CCA scores are generally orientated toward lowest Chla concentrations. Strongest association with the most Chla-poor OC1 is seen for Metridia lucens (3) and Calanus finmarchicus (1), which show a gradient of decreasing abundance with increasing OC type number. Candacia armata (4) was also found to be significantly higher in OC1 than OC2. Acartia spp (10), Temora longicornis (9), and Centropages hamatus (12) are strongly associated with Chla-rich waters, showing a gradient of increasing abundance with increasing OC type number, and are likewise oriented towards the region of higher Chla concentrations in the CCA biplot (Figure 7). Some species show a strong association with the seasonal axis in the CCA plot, including e.g. C. finmarchicus (1), which is oriented towards the cool waters and spring months, coinciding with the season associated with their ascendance to surface waters following overwintering.

Table 3. SIMPER Analysis Summary Showing Zooplankton Groups That Significantly Contributed to Community Dissimilarity Between Ocean Color Types in the Northeast Atlantic Shelves Area (8-Day Mean Match-Ups, n = 999 Permutations for Each Comparison)

ID	Species	Size class	OC type A	OC type B	Ave. Diss.	SD	Ave. Abu. A	Ave. Abu. B	p-value
11	Centropages typicus	Small	1	2	0.115	0.140	2.819	3.017	0.001***
8	Calanus spp	Small	1	2	0.112	0.130	2.788	3.273	0.001***
2	Calanus helgolandicus	Large	1	2	0.087	0.088	2.264	2.769	0.001***
14	Oithona spp	Small	1	2	0.082	0.122	1.906	1.920	0.001***
1	Calanus finmarchicus	Large	1	2	0.055	0.082	1.171	1.062	0.001***
3	Metridia lucens	Large	1	2	0.036	0.061	0.812	0.595	0.001***
13	Clausocalanus spp	Small	1	2	0.025	0.072	0.583	0.465	0.001***
4	Candacia armata	Large	1	2	0.012	0.029	0.262	0.217	0.001***
10	Acartia spp	Small	1	3	0.137	0.161	2.527	4.020	0.001***
11	Centropages typicus	Small	1	3	0.109	0.140	2.819	2.410	0.001***
9	Temora longicornis	Small	1	3	0.106	0.140	0.679	3.497	0.001***
2	Calanus helgolandicus	Large	1	3	0.083	0.088	2.264	2.028	0.001***
14	Oithona spp	Small	1	3	0.073	0.120	1.906	1.237	0.001***
1	Calanus finmarchicus	Large	1	3	0.052	0.084	1.171	0.691	0.001***
3	Metridia lucens	Large	1	3	0.032	0.061	0.812	0.279	0.001***
13	Clausocalanus spp	Small	1	3	0.019	0.067	0.583	0.177	0.001***
10	Acartia spp	Small	1	4	0.152	0.168	2.527	4.591	0.001***
9	Temora longicornis	Small	1	4	0.124	0.150	0.679	3.991	0.001***
1	Calanus finmarchicus	Large	1	4	0.051	0.085	1.171	0.568	0.001***
3	Metridia lucens	Large	1	4	0.031	0.061	0.812	0.196	0.001***
9	Temora longicornis	Small	2	3	0.109	0.134	1.927	3.497	0.001***
12	Centropages hamatus	Small	2	3	0.020	0.061	0.338	0.579	0.001***
10	Acartia spp	Small	2	4	0.143	0.156	3.488	4.591	0.001***
9	Temora longicornis	Small	2	4	0.122	0.142	1.927	3.991	0.001***
12	Centropages hamatus	Small	2	4	0.022	0.065	0.338	0.608	0.015*
15	Corycaeus spp	Small	2	4	0.016	0.058	0.237	0.422	0.016*
10	Acartia spp	Small	3	4	0.151	0.165	4.020	4.591	0.001***
9	Temora longicornis	Small	3	4	0.139	0.155	3.497	3.991	0.001***
12	Centropages hamatus	Small	3	4	0.028	0.074	0.579	0.608	0.001***
15	Corycaeus spp	Small	3	4	0.020	0.066	0.373	0.422	0.001***

• Note. Only top contributing species to average dissimilarity (Ave. Diss.) between group comparisons are included in summary table (95% cumulative dissimilarity cut-off, SD is standard deviation of dissimilarity). Average abundances (Ave. Abu.) are log-transformed values. Group identifiers (A vs. B) provide the specific OC type numbers compared (e.g., 1 = OC1, 2 = OC2, etc.). P-value significance levels are denoted by star symbols (* ≤ 0.05, ** ≤ 0.01, *** ≤ 0.001).

3.3 Drivers of Ocean Color Province Variability

The variations in OC province coverage in the North Sea can be seen in Figure 8, along with interannual and seasonal terms identified by the GAMs. Atlantic-associated OC1 and OC2 show increased North Sea coverage during summer and late summer, with OC2 peaking around June and OC1 peaking around August. Coastally-associated OC3 and OC4 show a bimodal seasonal signal with higher spatial coverage in spring (March) and later summer (August-September). Interannual variability is most notable during the period post-2014, where a strong decrease in the spatial coverage of Atlantic-associated OC types was observed from 2014-2016, followed by a strong increase. OC1 also shows lower areal coverage at the beginning of the time series. The leading EOF of each OC type shows the main seasonal spatiotemporal pattern in distribution (Figure S5 in Supporting Information S1).

For simplicity, only the RDA1 patterns for the more Atlantic-associated OC1 and OC2 types are presented here (Figures 9 and 10), but RDA results for all OC types are available in the Supplementary Information (Figures S6-10 and Tables S8-17 in Supporting Information S1). According to the permutation tests, CUR anomalies were significant predictors of the OC type anomalies, explaining between 15-20% of total variance, depending on OC type. The leading RDA1 patterns of all OC types were significant at the p ≤ 0.05, with RDA2 being only marginally significant for OC1. It should be noted that the cumulative explained variance of the significant RDA axes is lower than that of the full RDA model. Typically, this might indicate a need to remove insignificant explanatory variables from the constraining matrix, although this is not possible in the current application given that the PCs are component parts of a single field. In this special application, the total explained variation is an appropriate measure of the model's overall predictive skill, although only the leading RDA patterns are likely associated with larger-scale features in the OC field. The RDA1 patterns for the Atlantic-associated OC1 and OC2 types were largely opposite (out of phase) to those of OC3 given that periods of increased Atlantic inflow would naturally displace the more coastally-associated OC3.

Positive temporal RDA scores, which are associated with increased occurrence and extent of OC1 and OC2 in the North Sea, are particularly associated with a strengthened Norwegian Trench Inflow, located along the western edge of the Norwegian Trench, and a weakened Fair Isle Current (FIC), flowing southward between Fair Island and mainland Shetland (Figures 9 and 10, top-left panels). A closer examination of the original current fields in both of these areas confirmed that negative anomalies reflect current weakening, with only rare cases of a reversal in direction. Current reversals of the dominant cyclonic circulation pattern in the North Sea are associated with occasional periods of easterly winds (Sündermann & Pohlmann, 2011). The increased occurrence of OC1 east of the Shetlands is associated with a lower occurrence of OC2 in that area as its distribution is shifted further to the south. The normally eastward flowing Dooley Current, found north of the 100 m depth isoline, also shows weakened flows during these positive RDA1 phases. This weakening is associated with a further southward penetration of the more coastally-associated OC2 waters along the Scottish coast and into the southern North Sea (Figure 10, upper right panel). The more Atlantic-associated OC1 (Figure 9, upper right panel) is mainly associated with an intensification of the Norwegian Trench Inflow, which shows southward oriented CUR anomalies (Figure 9, upper left panel). Decreased occurrence of OC1 and OC2 are linked with CUR anomalies in the opposite direction (i.e. flipped by 180°), which are especially associated with the period of 2014-2017 (Figures 9 and 10, lower panels). As stated above, the more coastally-associated OC3 patterns were largely opposite in their dynamics (Figure S8 in Supporting Information S1), especially when compared to a subsequent RDA using a field of combined OC1 and OC2 occurrence (OC1&2, Figure S10 in Supporting Information S1). This can be seen in the fact that the spatial patterns are opposite in sign (i.e. indicated by color). OC4, which is associated with the littoral zone, showed more unique dynamics (Figure S9 in Supporting Information S1), but with clear temporal similarities to OC3. The spatial RDA pattern for OC4 showed a greater influence of CUR anomalies associated with English Channel inflow and northward propagating continental coastal waters.

Spearman correlation coefficients (ρ) between the dominant yearly indices of OC type anomalies (PC1) and other environmental indices are presented in Table 4. Generally, correlations were higher at lag 0 than lag 1, and of similar magnitude for the NAO and SPG indices. The sign of the correlations is opposite between groups of Atlantic-associated OC types (OC1, OC2, OC1&2) and the more coastally-associated OC3. OC4 had the weakest correlations to the tested environmental indices overall, and showed the most ambiguous interannual pattern of the OC types. One exception was the correlation between OC4 at lag 1 and the NAO (ρ = −0.55).

Some of the strongest correlations were observed between the SPG index and both OC1&2 (ρ = 0.48) and OC3 (ρ = −0.46) indices, both at lag 0. The NAO index was most strongly correlated with OC2 indices at lag 0 (ρ = −0.43). A plot of the time series from the strongest correlation between the SPG index and OC1&2 shows that the strongly negative OC1&2 anomalies of 2014–2017 are related to a stronger subpolar gyre (i.e. negative SPG index) (Figure 11). The largest outliers of the relationship are from the most recent years, where OC1&2 anomalies have returned to positive values while the SPG index has recovered but remains negative.

4.1 Ocean Color Provinces and Their Characteristics

Longhurst (2007) included the North Sea with the much larger Northeast Atlantic Shelves Province ("NECS"), which encompassed the continental shelf of Western Europe, from Cape Finisterre in NW Spain to the Skagerrak north of Denmark, and the Baltic Sea. Nevertheless, he acknowledged that it could be further divided into smaller partitions, including one for the North Sea that extends from the Straits of Dover to the Shetland Islands. This definition more or less coincides with our delineation of the North Sea (Figure 2) and largely corresponds to the extent of the Greater North Sea ecoregion defined by ICES (ICES, 2004, 2020). This ecoregion has been delineated with ecosystem-based fisheries management in mind, including political, social, economic and management considerations. Several other works have defined even smaller divisions for the North Sea based on hydrographic regimes (Becker & Pauly, 1996; Laevastu, 1962; Lee, 1980; van Leeuwen et al., 2015) and fish assemblages (Ehrich et al., 2009; Frelat et al., 2017). Our results also indicate that finer province delineations are likely more appropriate and that the use of dynamic boundaries, which could be monitored through the direct observation of OC types, could provide new insight into the hydrographic and biological changes of over time.

The four described OC types are consistently ordered along gradients of several abiotic and biotic parameters (Figure 6). These differences relate largely to whether the water is coastally- (OC3 and OC4) or offshore-associated (OC1 and OC2). OC4 corresponds to waters with the highest concentrations of Chla, POC and TSM, and the most nearshore distribution (Figure 4). Its spectral shape (Figure 3) is typical of very turbid water, which often contains high concentrations of mineral particles and dissolved organic matter (Babin & Stramski, 2004; Vantrepotte et al., 2012). Of the other abiotic parameters, SSS was the most clearly differentiated among OC types. A particularly narrow distribution of values is associated with OC1 (SSS above 35 g kg⁻¹), making it the most likely choice to use as an indicator of waters with Atlantic origin. SST was not distinguishable among OC types, although was significantly related to changes in the zooplankton community (Figure 7, Table 2). This discrepancy may be due to the complexity of SST patterns in 3D (e.g. dynamic of the stratification and heat exchange with the atmosphere, currents at different depth) and may more relate to seasonal patterns, such as the vertical migration of plankton (e.g. diel vertical migration).

OC3 and OC4 show a bimodal seasonal pattern with highest spatial coverage during the spring and late summer. These periods correlate with the typical annual pattern of phytoplankton blooms, which raises the question of whether the OC type classification is registering an increase in the coverage of low-Chla waters (OC1 and OC2) in summer simply due to the diminishment of blooms rather than variability in the inflow of Atlantic waters. Support for bloom depletion lies in the fact that OC1 and OC2 are associated with the area of seasonally-stratified waters (van Leeuwen et al., 2015), which may hinder nutrient resuspension to the photic zone while simultaneously facilitating losses via increased heterotrophic grazing (sensu the “Dilution–Recoupling Hypothesis” of Behrenfeld, 2010). On the other hand, support for Atlantic origin of these waters comes from the results of the coupled RDA patterns, which show a clear overlap between current variability and pathways of inflow for these Atlantic-associated OC types. Furthermore, OC1 and OC2 are clearly associated with salinities above 35 g kg⁻¹ (Figure 6), which is a signature of the Atlantic waters in the North Sea (e.g. Lee, 1980). Despite this evidence, the two mechanisms are not necessarily mutually exclusive, and further research is required to address the stability of OC types over time.

Various studies have focused on the dynamics of zooplankton in the Northeast Atlantic, and a large majority have relied on the CPR data set given its spatiotemporal coverage (e.g. Beaugrand, 2009; Hélaouët et al., 2013). The use of zooplankton as a possible indicator of Atlantic inflow into the North Sea has previously been a subject of study (e.g. Beare et al., 2002; Corten, 1999). In particular, the copepod species Candacia armata and Metridia lucens were identified as the most likely indicators of Atlantic waters given that their trends in the northwestern North Sea aligned with known inflow events (Corten, 1999). However, there was a lack of synchronicity with abundance trends in the area to the west of Scotland, which was assumed to be the source of inflow via FIC. Given the results of our work, this apparent disparity is explained by a stronger association of these indicator species to OC1 (Table 3), which is particularly related to inflows between the Shetland Islands and the Norwegian Trench (Figure 9), whereas OC2 is more related to inflow through the FIC (Figure 10), and may be more influenced by shelf waters. The extent of OC1 is also highest in late summer and early fall (Figure 8), which coincides with the months when both species were highly abundant in the northwestern North Sea area (Corten, 1999).

OC1 contains the highest number of significantly differing zooplankton taxa (Table 3). In addition to the two aforementioned taxa, Calanus finmarchicus is also strongly associated with OC1, and its abundance clearly declines with increasing OC type number. C. finmarchicus is known to be an important species of the Atlantic Subarctic and Arctic provinces (Hélaouët et al., 2011; Hélaouët & Beaugrand, 2007) but can enter into the North Sea with intrusions of Atlantic waters and can therefore be used as an indicator species. Combined with the clearly separated optical and abiotic characteristics, these findings would again suggest that OC1 shows promise for use as a proxy for waters of Atlantic origin. Taxa associated with the more coastally-associated OC3 and OC4 include Centropages hamatus and Temora longicornis, which have been previously described as neritic indicator species (Beare et al., 2002).

Although the OC types identified in this study show clear links to biotic and physical characteristics over the study period, extrapolation of these associations over longer time scales should be done with caution. Inferring links to other bio-optically derived characteristics (e.g. Chla) is likely to be robust given the common input data used for estimation (e.g. reflectance spectra); however, the association of OC types to specific plankton species may be less reliable given that they are not estimated remotely and may experience shifts in distribution over time that are not reflected by ocean color changes. For example, C. finmarchicus has shown a significant decrease in its abundance in the North Atlantic and North Sea (Beaugrand et al., 2003; Corten, 1999; Hélaouët et al., 2013; Pitois & Fox, 2006), which is also projected to continue in the coming decades (Hélaouët et al., 2011). This is linked to increased temperatures and a regime shift in the northeast Atlantic around the mid-1980s, which has resulted in a northward shift in its preferred habitat (Beaugrand, 2004; Hélaouët et al., 2013; Hélaouët & Beaugrand, 2007). However, planktonic aggregates, based on functional types or size, may show more consistency in relation to province dynamics on the regional scale of the North Sea. Ongoing work has shown that multispectral resolution data, such as were used in this study, are effective in identifying a suite of pigments and pigment groups (Bracher et al., 2015) and have been effectively used to quantify major phytoplankton groups on a global scale (Xi et al., 2020).

4.2 Insights Into North Sea Hydrography

The ability to observe OC provinces remotely and in real-time makes their monitoring an attractive alternative over other more difficult to observe parameters. In terms of the oceanographic dynamics of the North Sea, spatiotemporal data on currents and salinity would also be obvious choices for the monitoring of water masses; however, their observation is expensive and therefore sparse in space and time. Hydrodynamic models offer valuable alternatives; yet, while their performance is steadily improving, some processes are still challenging to resolve (e.g. ocean-shelf exchange). Furthermore, the lag in new data releases may affect their utility for short-term forecasting. The links between currents and the dynamics of OC types in the North Sea potentially allow for a more straightforward and economical approach to monitoring inflow variability and the spatial distribution of water masses.

Our results indicate that the dynamics and extent of OC types in the North Sea are likely driven by a combination of large-scale circulation patterns of the eastern North Atlantic and local wind-driven processes. Interannual dynamics in OC provinces in the North Sea were shown to be correlated with the large-scale SPG index (Table 4). The strength of SPG modulates the warm, north-eastern limb of the Atlantic Meridional Overturning Circulation (AMOC), which in turn influences the flow and properties of waters onto the northwest European shelf and into the North Sea (Huthnance et al., 2016). Koul et al. (2019) found that weak phases of the SPG (i.e. positive indices) result in more saline water from the subtropical North-Atlantic dominating waters flowing into the North Sea, whereas less saline subpolar inflow water is typical during the periods of a strong SPG (i.e. negative indices). This is supported by the strongly positive correlation between the SPG index and interannual anomalies in OC1&2 (Table 4), which are associated with higher salinity (Figure 6). During the time period considered in our study, the SPG was relatively strong in the early 2000s (negative index), weakened gradually till 2010 (positive indices) and spun-up again afterwards (Figure 11). This strengthening (acceleration) of the SPG culminated during the period of 2014-2016, resulting in the lowest index measured since the existence of satellite altimetry (since 1993) (Chafik et al., 2019). The same period coincided with the lowest extent of Atlantic-associated OC types in the North Sea (Figures 8-10), whereas the peak sea level values of the SPG index in 2006 coincided with one of the highest extents of these OC types.

While variability in the SPG helps to explain variability in source waters arriving to the northwest European shelf, local-scale atmospheric processes influence the inflow and extent of these waters into the North Sea proper. In particular, the seasonal and interannual variability of OC1 and OC2 coverage (Figure 8) seem to be driven by wind- and density-driven currents in the northern and central North Sea. On seasonal scales, their coverage increases during the summer months, coinciding with periods of lowest wind (Sušelj et al., 2010) and water current velocities (Mathis et al., 2015). It has been suggested that the role of the baroclinic (i.e. density-driven) component of the Atlantic water inflow increases compared to the wind-driven component in summer (Dooley, 1983; Huthnance et al., 2016; Mathis et al., 2015).

Contrary to the positive correlation with the SPG index, a negative correlation was found between the extent of these OC types and the NAO index (Table 4). These inverse correlations are not surprising, since the NAO has been shown to drive to a large extent changes in the strength of the SPG and, thereby, the water properties in the eastern North Atlantic (Bersch et al., 2007). Curry and McCartney (2001) showed that the intensity of the SPG is a time-integrated oceanic response to persistent NAO anomalies. For the period considered in our study, the correlation between both indices was −0.47 (p < 0.05). Positive NAO indices are associated with enhanced westerly winds, increased inflow into the North Sea, and an intensified cyclonic (anti-clockwise) circulation (Emeis et al., 2015; Mathis et al., 2015). These patterns are reflected by the results of the RDA, which show that the southward penetration of Atlantic waters is associated with a weakening of the Dooley current and an intensification of the Norwegian Trench inflow (Figures 9 and 10). Using a modeling approach, Mathis et al. (2015) estimated that, on average, a large majority of the FIC and of the inflow east of the Shetland Islands (76%) turn eastward along the 100 m isobath of the Dooley Current, while 24% of inflows penetrate further southward to form the central recirculation cell north of the Dogger Bank. Of these waters, only 5% penetrate to the most southern reaches of the North Sea, where they merge with Atlantic waters entering via the English Channel. Our results suggest that the weakening of the Dooley current, in particular, allows for the penetration of OC2 into these less-commonly reached southern areas of the North Sea, around the limits of the 40 m isobath and just south of the Dogger Bank (seen in Figure 10 for OC2 and Figure S10 in Supporting Information S1 for OC1&2). The fact that this penetration increases during periods of reduced wind-forcing is further support of the importance of baroclinic forces for Atlantic inflow on the interannual scale (Núñez-Riboni & Akimova, 2017).

Based on these findings, we propose that OC provinces have predictive skill in identifying water mass characteristics of the North Sea, which are largely driven by the interaction between large-scale circulation of the eastern North Atlantic and the wind-driven and baroclinic local currents. Although the two processes operate on differing spatial scales, a degree of synergy was observed during the study period, with both the strength of the SPG and the local wind-forced currents correlating with OC province dynamics. This is supported by the existence of strong external control on North Sea circulation (Holt & Proctor, 2008; Huthnance, 1997). Recent work of Pätsch et al. (2020) also shows that a similar interplay of regional and local processes helps to explain variability in salinity and nutrient concentrations in the North Sea. Specifically, while the SPG has been shown to affect the variability in salinity of waters entering into the northern North Sea following an advective delay of up to two years (see also Koul et al., 2019), nutrient concentrations are more immediately affected. The authors propose that atmospheric drivers are more likely driving nutrient dynamics. Given that OC province dynamics also show a higher correlation with unlagged indices of the SPG and NAO indices (Table 4), as well as their association to phytoplankton concentrations (Chla, Figures 6 and 7), it is likely that their spatial variability reflects similar changes in water mass characteristics that can influence primary production, such as nutrient concentrations.

4.3 Ecological Implications and Prospectus

The impact of OC province dynamics on ecological processes is likely to be a complex mixture of trophic and physical effects. The correlation between OC provinces and associated bio-physical properties has provided evidence of a link to lower trophic levels, which, in turn, are likely to have bottom-up effects to other components of the ecosystem. Furthermore, the spatially explicit perspective sheds new light on the dynamics of Atlantic waters inflow to the North Sea in terms of coverage rather than the more typical focus on (volumetric) fluxes. Both aspects are likely to be related, albeit inversely depending on where inflow is measured; however, the areal extent of OC provinces may be of particular relevance to ecological processes that are spatial in nature. For example, positive inflow anomalies in 1988/89 and 1998 were claimed to have been the cause of increased phytoplankton biomass in the central North Sea through enhanced nutrient availability (Reid et al., 2001). Our results suggest that these increases may not be a direct consequence of the increased influx of Atlantic waters, but rather due to a larger relative coverage of the more Chla-rich provinces, OC3 and OC4, which expand their distribution further northward as currents increase and Atlantic water penetration is reduced.

From this spatial perspective, monitoring OC provinces may further elucidate other ecologically-relevant processes in higher trophic levels, such as habitat changes or the transport of eggs and larvae from broadcast spawning species to important nursery areas. The area of inflow in the northwestern North Sea is considered an important spawning area for cod in spring (González-Irusta & Wright, 2016) and herring in fall (e.g. Dickey-Collas et al., 2010), and is generally associated with highest catches for many gadoid fish species (ICES, 2019). This habitat preference may be related to the observed differences in the zooplankton community associated with the Atlantic waters. For example, cod larvae and early juvenile stages are known to prefer larger calanoid copepods like C. finmarchicus, and decreases in cod recruitment in recent decades has been in part attributed to a possible mismatch in the spatial distribution of these zooplankton (Beaugrand et al., 2003).

Such mismatches may help to explain the large degree of recruitment variability in many fish stocks, which has implications for management and the ability to provide reliable short-term advice in the form of future catch quotas. The incorporation of environmental indices in models of historical stock-recruitment relationships has demonstrated significant improvements in fit (i.e. in addition to spawning stock biomass), and the use of more spatially-focused, smaller-scale environmental signals may be of particular importance to understanding recruitment dynamics (Akimova et al., 2016). Novel statistical approaches (e.g. Machine Learning) have also shown promise in the identification of environmental drivers for recruitment from large data sets of spatially explicit data (Kühn et al., 2021). The ever-growing suite of available environmental data sets, including ocean color, will likely improve our ability to explore environmental drivers to resource dynamics in the future. The ability to monitor OC provinces in real-time provides a unique operational index for use in other modelling and forecasting applications.