Trends computed with SST data show a significant (p < 0.001) overall warming tendency for the CanC EBUE during the 1993–2014 time period (Figure 2), although the increase in temperature is more evident in the southern part of the study area (around 0.4°C·decade⁻¹) rather than in the northern one (about 0.2°C·decade⁻¹). Nevertheless, in contrast to the open-ocean, nearly all of the water located over the shelf (<200 m depth) shows no significant change in temperature. Furthermore, some spots present slightly negative trends, coinciding with some of the major upwelling cores, i.e., Cape Blanc (21°N), Dakhla (24°N), and Cape Ghir (30°N). At subregional level (Table 1), the WPUZ and PUZ present similar non-significant warming trends over the shelf (0.09 and 0.13°C·decade⁻¹, respectively), whereas the SUZ exhibits a significant average trend of 0.3°C·decade⁻¹. However, in all three instances, the warming trends in the open-ocean and slope areas are higher than over the shelf.

Computed NPP trends show very different results depending on the model and the dataset (which comprise both the time period and sensor factors). Both the VGPM (Figure 3) and Eppley-VGPM (Figure 4) exhibit similar trends, although these are more marked in the former. For the SeaWiFS dataset two patterns can be clearly identified: an increasing trend over the shelf near Cape Blanc (around 600–800 mgC·m⁻²·decade⁻¹) and a decreasing one around it. Nevertheless, note that only in some patches the trends are significant (p < 0.05). Regarding the two other upwelling subregions, both present zones of slight increase/reduction of NPP. Average trend values for the VGPM over the shelf in the SUZ, PUZ and WPUZ are −75, 105, and 10 mgC·m⁻²·decade⁻¹, respectively, whereas the Eppley-VGPM shows increases of 175, 160, and 70 mgC·m⁻²·decade⁻¹; in neither case are they significant (Table 1). On the other hand, the MODIS datasets are dominated by decreases in NPP (Figures 3, 4). VGPM mean trends over the shelf are −530, −285 and −25 mgC·m⁻²·decade⁻¹ in the SUZ, PUZ and WPUZ, respectively, the first two being significant (p < 0.05). Similarly, the Eppley-VGPM yields non-significant changes in NPP of −355, −155 and 0 mgC·m⁻²·decade⁻¹ for the same zones (Table 1).

In general, the CbPM (Figure 5) exhibits different trends relative to the chlorophyll-based VGPM models. The SeaWiFS dataset registers a marked significant decreasing trend exceeding −400 mgC·m⁻²·decade⁻¹ (p < 0.001) over most of the WPUZ and PUZ (Figure 5 and Table 1). On the other hand, the SUZ barely shows any change. Nevertheless, it is worth noting the highly significant (p < 0.001) positive trend off Cape Blanc, with broad areas exceeding 200 mgC·m⁻²·decade⁻¹. The MODIS dataset yields different tendencies: no apparent changes are registered for most of shelf waters over the entire study area (Figure 5). However, mean trends result in significant decreases of −100 and −110 mgC·m⁻²·decade⁻¹ for the SUZ and WPUZ, respectively, and a non-significant drop of −100 mgC·m⁻²·decade⁻¹ for the PUZ (Table 1). Off the continental break, significant declines in NPP are observed widespread (Figure 5).

These results reveal that absolute values and trends in NPP greatly vary depending on the chosen dataset and, especially, model. Absolute NPP mean values show marked differences between Chl-based models and CbPM. In the most productive areas (e.g., Cape Blanc) these discrepancies are magnified, VGPM/Eppley-VGPM showing fourfold NPP values relative to CbPM. Regarding decadal change rates in NPP, a detailed comparison of the average trends in the shelf portion of each of the upwelling zones (Table 1) is shown in Figure 6. For the SeaWiFS dataset, the VGPM and Eppley-VGPM show no significant changes, while the CbPM presents remarkable reductions in the WPUZ and PUZ regions (Table 1). The three models only agree in the SUZ, where NPP estimates do not show any clear trend. Nonetheless, a common feature in the VGPM, Eppley-VGPM and CbPM is the marked seasonality of NPP in all three upwelling zones, reflecting the annual cycle in solar irradiance. On the other hand, for the MODIS dataset all three models agree fairly well in a qualitative sense, exhibiting decreases in all the upwelling zones. Notably, CbPM trends shift from marked decreases in SeaWiFS data to more gentle ones in the MODIS dataset. VGPM and Eppley-VGPM also show changes among SeaWiFS and MODIS datasets, as overall increases are registered in the former and decreases in the latter. When comparing SeaWiFS and MODIS trends in their shared period (2003–2007), both VGPM and Eppley-VGPM show trends that qualitatively match, except in some areas between the Canary Islands and the Cape Verde archipelagos (Supplementary Figure 1), where weak, non-significant trends shift from positive to negative (and viceversa). CbPM trends, on the other hand, exhibit large areas of opposing trends in the open ocean, but agree in the shelf area north of Cape Blanc (21°N). Thus, SeaWiFS and MODIS trends agree reasonably well between 2003 and 2007 for areas close to the upwelling. However, in terms of absolute NPP values it is interesting to notice that SeaWiFS always yields higher values for CbPM than MODIS, whereas the opposite is true for the VGPM and Eppley-VPGM models.

Chl-a trends exhibit no significant changes except for some areas north of Cape Ghir (30°N) and around Cap Vert (14°N) (Supplementary Figure 2). This is due to the fact that increases in the first part of the time-series are followed by a stagnant period until ~2010, when the Chl-a starts falling. Average OC-CCI Chl-a trend values over the shelf differ between upwelling zones, being 0.09, 0.17, and 0.83 mg Chl·m⁻³·decade⁻¹ in the WPUZ, PUZ, and SUZ, respectively (Table 1), although only the first one is significant at p < 0.05. However, time-series (Figure 7) seem to show an initial increasing period followed by a decreasing last one. When looking at the individual Chl-a datasets differences are observed between them. SeaWiFS Chl-a shows no significant changes, with trends of 0.06, 0.18 and −0.03 mgChl·m⁻³·decade⁻¹ in the SUZ, PUZ, and WPUZ, respectively (Table 1). On the other hand, for MODIS Chl-a significant negative trends (p < 0.01) of −1.66 and −0.38 mgChl·m⁻³·decade⁻¹ are registered in the SUZ and PUZ, respectively, whereas no significant variation is observed in the WPUZ. Overall, OC-CCI exhibits lower Chl-a concentrations than MODIS and SeaWiFS, although being very close to the latter from 2002 to 2007.

Computed SST trends reflect a non-surprising overall warming of the ocean, which has been extensively described in the literature over the last years (Levitus et al., 2000, 2005; Hansen et al., 2006, 2010; Domingues et al., 2008; Lyman et al., 2010; Gouretski et al., 2012). Warming trends calculated in the present work for this specific oceanic region are slightly weaker than those found in the literature. Demarcq (2009) used AVHRR pathfinder v5 SST data to estimate trends for the 1998–2007 period. Although the pattern is similar, his average trends were stronger than those obtained in the present work. This might be explained in part by the fact that (1) we used a longer time period for our analysis, (2) he selected different areas for his study, and (3) he used a different method to estimate the trends (Least Absolute Deviations (LAD) instead of the Theil-Sen method employed here). Nevertheless, he did obtain lower warming trends when these were computed for areas with 0–1,000 m depth, i.e., relatively close to the coast, a fact that qualitatively matches the findings presented here. Moreover, in the present work negative trends are registered over the shelf at some locations. This could mean that intensified upwelling of subsurface waters might be occurring, resulting in gentle decreases in SST as a consequence of larger amounts of cool waters reaching the ocean surface. In any case, this would be limited to certain reduced locations (Figure 2B).

McGregor et al. (2007) examined two historical sediment cores (extending back 2,500 years) from Cape Ghir, using the alkenone unsaturation index (U37K′) as a proxy for SST. For the most part of the twentieth century the sediment cores showed a SST decrease of 1.2°C. Narayan et al. (2010) also studied the Cape Ghir region within the CanC upwelling system as part of their analysis on SST and wind trends in the four major EBUEs. Comparing several wind datasets they obtained contradictory results on meridional wind trends, although the COADS (Comprehensive Ocean Atmosphere Dataset) dataset, the most reliable to their view, exhibited an equatorward wind increase in the CanC. However, caution is required when extrapolating results from Narayan et al. (2010) to the whole CanC EBUE given that the SST upwelling index was estimated exclusively in the northern portion of the WPUZ and meridional wind stress time series correspond to a small region of 3° × 5° of a non-specified location within the upwelling system.

Barton et al. (2013) questioned the results of the abovementioned works. They argued that the use of U37K′ as a proxy of SST could lead to results biased toward lower SST because of coccolithophores (from which U37K′ signal is derived) living in deeper layers; i.e., as the ocean warms and stratifies, phytoplankton will tend to occupy deeper parts of the water column and in consequence have less contact with the mixed layer temperatures. Thus, estimating SST from U37K′ as McGregor et al. (2007) did would bias measures toward lower SSTs than the real ones, as results would truly correspond to cooler subsurface waters. Like Narayan et al. (2010), and Barton et al. (2013) suggested that the evidence for increased upwelling (a greater SST difference between the coast and the open ocean) is not based on a cooling of the coastal waters but on a greater warming of the open ocean waters, as observed in the present work. Furthermore, Barton et al. (2013) analyzed several wind and SST datasets and found that there were not significant changes in meridional wind intensity off NW Africa but a significant SST increase. Consequently, they did not find any strong evidence supporting Bakun's hypothesis.

Cropper et al. (2014) estimated SST trends from HadISST (Hadley Centre Sea Ice and Sea Surface Temperature) and OISST (Optimum Interpolation Sea Surface Temperature) datasets for the summer months (June to August), i.e., when the trade winds are strongest. Their findings showed a very weak negative SST trend for the PUZ, although the trend was not significant even at the 0.1 level. They also found a significant increase in the equatorward meridional wind in several datasets but, once again, the significance limit is set at the 0.1 level, challenging the interpretation of the results.

NPP and Chl-a can be used as proxies for upwelling intensity because, although they depend on several factors, upwelling of deep waters is the primary driver of the amount of available nutrients for phytoplankton, which in turn control NPP and Chl-a concentration (Ohde and Siegel, 2010; Messié and Chavez, 2015).

Trends for the times series of the three NPP models present quite different results. Seasonal to decadal changes in the VGPM and Eppley-VGPM are nearly exclusively located around and south of Cape Blanc, and with the exception of some isolated patches they are not statistically significant (Figures 3, 4). Notably, the robustness of the data in the southern subregion of the CanC is affected by the high occurrence of cloud coverage (even on a weekly basis, Supplementary Figure 3), which is particularly persistent during summer months. Furthermore, the PUZ and WPUZ exhibit no significant trends at all. This seems to indicate that no variation in upwelling intensity is experienced for the studied period. On the contrary, the CbPM shows large areas with significant opposite trends (Figure 5). The sharp, significantly negative trends observed with the SeaWiFS dataset in the shelf portion of the PUZ and WPUZ suggest an important decadal decrease in upwelling activity. On the other hand, the MODIS dataset shows very weak trends. Discrepancies in absolute values yielded by the models and datasets also contribute to this conundrum. However, when trends are qualitatively compared for the common SeaWiFS and MODIS period overall they agree reasonably well, with differences only arising where weak, non-significant trends are registered (Supplementary Figure 1). The quantitative shift from lower NPP estimates during the SeaWiFS period to higher ones during the MODIS period observed for VGPM and Eppley-VGPM in shelf areas (Figure 6) is presumably consequence of differences in Chl-a retrievals by these two sensors. Arun Kumar et al. (2015) reported that MODIS overestimated Chl-a relative to SeaWiFS in coastal, eutrophic waters of the Arabian Sea. Demarcq and Benazzouz (2015) found similar results for the CanC region, reporting that MODIS underestimated Chl-a values relative to those estimated by SeaWiFS for low for concentrations (~<0.4 mg·m⁻³) while overestimated Chl-a for high concentrations (~<1 mg·m⁻³). The latter case is relevant when estimating Chl-a concentrations in upwelled waters of the CanC, which largely exceed that value. However, these differences occur consistently and, thus, qualitative comparisons of trends should be acceptable. As outlined by Demarcq and Benazzouz (2015), these discrepancies between SeaWiFS and MODIS estimates arise from differences in (1) calibration and atmospheric correction, (2) correction of sensor sensitivity drifts, and (3) constraints in data processing, e.g., cloud masking, which is of importance particularly in the southern part of the CanC, where cloud coverage is considerable (Supplementary Figure 3). The combination of these issues would also explain the observed differences in NPP, as estimates from the VGPM models are closely related to the Chl-a concentration. Thus, discrepancies between models seem to remain the major issue. In situ primary production (PP) measurements with which to compare model NPP estimates are scarce in the CanC system. Arístegui et al. (2004) compiled available PP data from coastal waters of the CanC upwelling system found in literature, but the bulk of the studies were carried out during the 1970s and 1980s, i.e., before the start of the NPP time-series and thus they cannot be accurately validated. However, model estimates fall between the observed in situ measurements, which ranged between ~1,000–4,000 mg C·m⁻²·day⁻¹, depending on location and season. Yearly mean in situ PP for 31.5–18°N yielded values of 2,400 ± 1,500 mg C·m⁻²·day⁻¹. If compared to mean values from this work (Table 4), Eppley-VGPM is the model that gets closest to in situ measurements, yielding 2,139 ± 610 and 2,731 ± 769 mg C·m⁻²·day⁻¹ for SeaWiFS and MODIS, respectively. On the other hand, VGPM presents higher values (~3,000–3,700 mg C·m⁻²·day⁻¹) whereas CbPM exhibits markedly lower NPP (~500–800 mg C·m⁻²·day⁻¹). This result (VGPM presenting greater NPP values than CbPM) contrasts with the findings of Westberry et al. (2008), who reported the opposite for the global ocean. Their regional results also show greater values for CbPM in gyres and subtropics, whereas VGPM exceeds CbPM estimates in high latitudes. The present work provides evidence that an independent comparison for EBUEs is needed as models show a different behavior to that exhibited in the broad subtropical domains.

Kahru et al. (2009) showed that in the California Current region the VGPM adjusted better to in situ NPP measurements (obtained from the large CalCOFI data base) than the CbPM (r² = 0.662 vs. r² = 0.389, respectively), although both models seemed to overestimate NPP. The disagreement between the model outputs, which is also observed in the present work, might arise from dissimilarities in the way the models estimate NPP, as well as in errors from the input data. As mentioned by Kahru et al. (2009), the VGPM (and Eppley-VGPM) requires three input fields while the CbPM needs five (six in the modified version), thus increasing uncertainties in derived NPP. Moreover, the authors argue that the CbPM depends on the GSM algorithm to derive Chl-a and b_bp estimates from water leaving radiance, which is a source of errors especially near the coast (mainly due to colored dissolved organic matter and suspended sediments) and when complex atmospheric conditions are present (Kostadinov et al., 2007). Additional error sources in the CanC upwelling system may arise from the phytoplanktonic community composition and the effect of dust on remote sensing measurements. In parts of the NW African coastal upwelling, where mixotrophy can certainly be a major feature of the planktonic community (Anabalón et al., 2014), Chl-a might not be the most adequate proxy to infer NPP: mixotrophic phytoplankton groups would contribute to increase carbon:Chl-a ratios, thus possibly biasing the output from NPP models based on Chl-a. However, the importance of mixotrophy regarding the entire CanC system remains poorly understood, as very few studies on phytoplanktonic community structure have been carried out. Another potential issue with Chl-a-based NPP estimates is the overestimation of Chl-a concentrations as a consequence of poor atmospheric correction. For instance, Diouf et al. (2013) argued that the SeaWiFS algorithm overestimates aerosol reflectance in the blue and green bands when weakly absorbing aerosols are present, producing lower sea surface reflectances and thus presumably yielding overestimated Chl-a concentrations. Considering the high transport rates of mineral aerosols to the study area from the Saharan desert and the Sahel (Mahowald et al., 2005), this could be a major issue. Hence, considering all the aforementioned potential error sources from input data it is difficult to contrast the performance of the NPP models themselves.

The marked seasonality exhibited by NPP in WPUZ and PUZ initially might not seem to fit with the upwelling patterns described in classical and recent literature. Wooster et al. (1976) and Van Camp et al. (1991) mentioned three similar zones regarding trade wind seasonality: (1) strong, year-round trade winds between 20 and 25°N, where permanent upwelling is present; (2) persistent trade winds during winter south of 20°N, which fade during summer; and (3) persistent trade winds during summer north of 25°N, which fade during winter. Cropper et al. (2014) followed a similar description except for the northernmost section of the upwelling system. North of 26°N they describe the aforementioned weak permanent annual upwelling zone where, despite the fact that a maximum in upwelling intensity is registered during summer months, its occurrence is year-round. Messié and Chavez (2015) concluded that upwelled macronutrients (especially Si) for the most part of the year, as well as light during winter, control NPP in the CanC upwelling system. Light availability patterns in the CanC EBUE roughly follow the same distribution of trade winds (Demarcq and Somoue, 2015). Between 20 and 26°N high light intensity is registered all year long, but out of this latitudinal band light availability decreases. North of 26°N light is reduced due to a more temperate climate, presenting a summer maximum; south of 20°N light is limited during summer months due to the increase in cloud cover consequence of the northward migration of the Intertropical Convergence Zone (ITCZ). Thus, NPP changes could be attributable mainly to the co-variability of trade winds (and, consequently, upwelled nutrients) and light: NPP peaks during summer in the PUZ and WPUZ, coinciding with the maximum of trade wind and light intensity, whereas in the SUZ highest NPP values are registered during late winter and spring, and fade in summer, following the changes in trade winds and light availability.

Demarcq (2009) estimated SeaWiFS Chl-a trends for the same period as the SeaWiFS dataset in our study. He obtained relatively similar trends, with small differences possibly due to the different areas selected for trend averaging and the distinct methods used to calculate the trends (Theil-Sen method against LAD method). In the two cases the overall decreasing trend in Chl-a (and also in NPP derived from the VGPM/Eppley-VGPM in our study) is less pronounced (or even positive) over the shelf. Demarcq (2009) suggested that this might be due to an increase in upwelling intensity or its intrinsic efficiency. However, he did not estimate the significance of trends, which probably were not significant, as in our study (Table 1). Furthermore, the MODIS dataset shows significant decreasing trends in Chl-a over the shelf in the SUZ and PUZ, which contrast with the upwelling intensification hypothesis. The OC-CCI presents positive trends, but these turned out to be non-significant. In fact, trends seem to be somehow biased by a jump in Chl-a concentrations around 2002. Afterwards, Chl-a concentrations seem to stabilize and start to decrease, matching the observed trend in the MODIS dataset. Differences in absolute Chl-a estimates between SeaWiFS/MODIS and OC-CCI are probably consequence of the merging procedure carried out to generate the latter.

In summary, NPP and Chl-a trends exhibit heterogeneous results depending on the dataset, production model and/or upwelling portion analyzed, but overall a pattern of no significant trends or significant decreases can be observed. Considering these results, one might infer either that the CanC has not been experiencing an intensification in upwelling intensity or that, in fact, a reduction is taking place, as suggested by the CbPM model.