Accuracy Assessment of Primary Production Models with and without Photoinhibition Using Ocean-Colour Climate Change Initiative Data in the North East Atlantic Ocean

Abstract

The accuracy of three satellite models of primary production (PP) of varying complexity was assessed against 95 in situ ¹⁴C uptake measurements from the North East Atlantic Ocean (NEA). The models were run using the European Space Agency (ESA), Ocean Colour Climate Change Initiative (OC-CCI) version 3.0 data. The objectives of the study were to determine which is the most accurate PP model for the region in different provinces and seasons, what is the accuracy of the models using both high (daily) and low (eight day) temporal resolution OC-CCI data, and whether the performance of the models is improved by implementing a photoinhibition function? The Platt-Sathyendranath primary production model (PP_PSM) was the most accurate over all NEA provinces and, specifically, in the Atlantic Arctic province (ARCT) and North Atlantic Drift (NADR) provinces. The implementation of a photoinhibition function in the PP_PSM reduced its accuracy, especially at lower range PP. The Vertical Generalized Production Model-VGPM (PP_VGPM) tended to over-estimate PP, especially in summer and in the NADR. The accuracy of PP_VGPM improved with the implementation of a photoinhibition function in summer. The absorption model of primary production (PP_Aph), with and without photoinhibition, was the least accurate model for the NEA. Mapped images of each model showed that the PP_VGPM was 150% higher in the NADR compared to PP_PSM. In the North Atlantic Subtropical Gyre (NAST) province, PP_Aph was 355% higher than PP_PSM, whereas PP_VGPM was 215% higher. A sensitivity analysis indicated that chlorophyll-a (Chl a), or the absorption of phytoplankton, at 443 nm (a_ph (443)) caused the largest error in the estimation of PP, followed by the photosynthetic rate terms and then the irradiance functions used for each model.

1. Introduction

The rate of synthesis of organic matter by marine phytoplankton through the process of photosynthesis determines the energy flow through the trophic chain in the global ocean. This process, known as primary production (PP), not only fuels biological growth and fish productivity, but also regulates carbon uptake and release by the ocean [1]. Since the 1950’s, PP has been determined using radio labelled carbon [2] through in situ or simulated in situ incubations. Though these methods have provided great insight into the variability in PP, the number of measurements available is limited both spatially and temporally. The introduction of rapid fluorescent based and O₂/Ar ratio methods using membrane inlet mass spectrometry has increased the number of measurements that can be made over time and space, but still cannot provide synoptic coverage at global, basin, or regional scales. The use of ocean colour remote sensing data has revolutionised our ability to monitor PP at higher temporal and spatial scales than can be provided by in situ techniques [3,4,5,6,7,8]. Over the past two decades, there has therefore been a concerted effort to estimate PP from the satellite ocean colour using models based on the following input parameters: (1) Phytoplankton biomass expressed as Chl a, carbon concentration, or phytoplankton absorption; (2) photo-physiology or photosynthetic rates; and (3) the light field. More than twenty models have been developed that are run solely using remotely sensed data. Some of these models have been validated to ascertain their relative or absolute accuracy [9,10,11,12,13,14].

Historic studies have shown that satellite derived PP accounts for 50 to 80% of the variance in measured PP depending on the input parameters used in the model [6,7,9,15,16]. Many of these historic studies are limited by a lack of in situ data to characterize the spatial, temporal, and vertical variability in photo-physiology that cannot be measured by remotely sensed parameters alone. Some parameters used in the model, including photosynthetic rates, are input as constants or as functions of other physical parameters, and their vertical distribution is assumed to be constant with depth, which can lead to significant errors in model performance [9,10,11,17]. In addition, some models fail in regions with atypical optical and photo-physiological properties [18,19,20] and therefore require a specific regional parameterization.

Early validation studies illustrated that the main cause of the difference between modelled and in situ PP was from parameters that characterize the absorption of light by phytoplankton through Chl a, phytoplankton absorption coefficients, or carbon biomass [21,22,23]. This resulted in an over-estimate of PP in eutrophic waters and an under-estimate in mesotrophic and oligotrophic waters [15]. A series of National Aeronautics and Space Administration (NASA) PP model inter-comparisons were conducted in the 2000’s, firstly by [9], who compared twelve different PP models in nine regions of the world ocean using 89 in situ PP measurements, and found that the difference between satellite and in situ estimates of PP was highly correlated with differences between in satellite and in situ Chl a (See also [24]). Errors in the PP model can be further compounded when other input parameters are derived from Chl a. The second NASA PP inter-comparison [10] illustrated that a four-fold variation of Chl a caused a three-fold variation in PP, whereas variations in photosynthetic active radiation (PAR) and sea surface temperature (SST) had smaller effects on PP. This study also showed however, that deriving photosynthetic rates from temperature alone can cause a significant error in PP. More recently, NASA compared twenty-one PP models in ten regions of the World Ocean using 1156 in situ PP measurements, and showed that the main cause of error in PP models was in coastal, optically complex waters where satellite Chl a can be inaccurate [13]. Further to this, Milutinović and Bertino [25] conducted a comprehensive analysis of the errors associated with the parameters used in one PP model (the Vertical Generalized Production Model - VGPM [6]) and showed that up to 90% of the error is associated with the biomass specific optimum photosynthetic rate (P^B_opt). The use of regional in situ values of photosynthetic parameters in satellite derived models improves the performance of PP models [26,27,28].

In this paper, we assess the accuracy of one Wavelength-Integrated (WIM) and two Depth-Integrated (DIM) models [29] that are commonly used in the literature, in the North Atlantic. The WIM tested is the Platt and Sathyendranath model (PP_PSM) [22], which uses an exponential Photosynthesis-Irradiance (P-I) function based on the absorption of light by the photosystems [30,31]. The most commonly used DIM is VGPM (PP_VGPM), which is simple and quick to apply to satellite data. By contrast, we also assess the DIM Absorption Based Model (PP_Aph) [32,33], which uses the optical absorption signature of phytoplankton as an empirical function of phytoplankton biomass to derive PP. To date, few models include a photoinhibition term that accounts for the photosynthetic rate with increasing irradiance that can occur under high light conditions or when the phytoplankton community moves suddenly from low to high light, which can cause a decrease in the photosynthetic efficiency [34].

Model performance was assessed using 95 in situ measurements of daily water column primary production (PP_eu) from the NEA. The models were run using the ESA, OC-CCI version 3.0 data. The paper addresses the following questions: (1) which is the most accurate model for this region? (2) Does the model performance improve when implementing a photoinhibition function? (3) Does the accuracy of each model change depending on the province and season? (4) What is the accuracy of the models using both high (daily) and low (weekly) temporal resolution data? (5) Which input parameter causes the greatest error in each model?

2. Materials and Methods

2.1. Study Region

Stations of PP_eu were sampled in the North East Atlantic (20–65° N, 5–40° W) from 1998 to 2013, covering the four biological provinces according to [35]: the ARCT, NADR, NAST, and NATR (North Atlantic Tropical Gyre province) (Figure 1,

Table 1. Provinces and time periods for in situ integrated daily water column primary production (PP_eu) (mg C m⁻² day⁻¹) used for the accuracy assessment of the primary production satellite models.

Region	North East Atlantic (NEA): 20–65° N, 5–40° W
Period	1998–2013
Number of stations (n)	95
Seasons 1
summer, n = 56	June-1998, July–August 2002, June-2003, May-2004, June-2005, July–August 2007, August-2009
autumn, n = 38	September-2003, September-2004, October-2008, September-2009, October-2011, October-2012, October-2013
Provinces	ARCT 2—Atlantic Arctic Province (n = 12)
	NADR—North Atlantic Drift (n = 28)
	NAST 3—North Atlantic Subtropical Gyre (East & West) (n = 40)
	NATR 4—North Atlantic Tropical Gyre (n = 14)

¹ When all stations were considered together, additional data from April 2002 were also included. ² Biogeographical provinces were defined following [35]. ³ Analysis of stations in NAST(W) and NAST(E) were merged into NAST. ⁴ One station from CNRY (Canary Current Coastal province) was added to stations in NATR.

). A total of 95 stations were analysed during summer (May–August) and in autumn (September–October). Based on OC-CCI Chl a during the sampling periods, 54% of the stations were in oligotrophic waters (0–0.1 mg m⁻³) and 42% were mesotrophic (0.1–1 mg m⁻³) [36].

2.2. Simulated In Situ ¹⁴C Primary Production Measurements (PP_eu)

PP_eu was determined on Atlantic Meridional Transects (AMT) 6, 12, 13, 14, 15, 16, 18, 19, 21, 22, and 23 [37], Ecosystem of the Mid-Atlantic Ridge (ECOMAR) cruises in the North Atlantic JC011 and JC037, and cruise D264 in the Celtic Sea [14]. Water samples were taken from 6 to 9 depths and transferred from Niskin bottles to black carboys to prevent shock to the photosynthetic lamellae of the phytoplankton cells. Water from each sample was sub sampled into three 75 mL clear polycarbonate bottles and three black polycarbonate bottles. The bottles were pre-cleaned following Joint Global Ocean Flux Study (JGOFS) protocols to reduce trace metal contamination [38]. Each sample was inoculated with between 185 and 740 kBq (5–20 µCi) NaH¹⁴CO₃ according to the biomass of phytoplankton and transferred to an on deck (simulated in situ) incubation system using neutral density and blue filters to simulate subsurface irradiance over depth to 97%, 55%, 33%, 20%, 14%, 1%, or 0.1% of the surface value. Incubations were conducted from dawn to dusk for 10 to 16 h. The upper four light levels were maintained at the in situ temperature by pumping water from the sea surface into the incubators, and the two lower light depths were chilled with fresh water from a water bath maintained at within ±3 °C of the in situ temperature. For all cruises, after incubating the water samples, suspended material was filtered through 0.2 µm polycarbonate filters at a vacuum pressure of <20 cm Hg. The filters were then rinsed with filtered seawater, exposed to concentrated HCl fumes, and immersed in a 5 mL scintillation cocktail. ¹⁴C disintegration time per minute (DPM) was measured on board on either a Beckman LS6000 SC or a TriCarb 3900 TR liquid scintillation counter (LSC) using the external standard and the channel ratio methods to correct for quenching. The quantity of ¹⁴C added to the experimental bottles was determined by adding aliquots of the stock ¹⁴C solution to Carbosorb, a CO₂ absorbing scintillation cocktail, and counted immediately on the LSC.

2.3. Satellite Data

To estimate PP, satellite ocean colour, SST, and PAR data were used. Each model was firstly run using daily OC-CCI v.3.0 data. This was then compared using eight day composite OC-CCI data, which is freely available to end-users through the internet. From the OC-CCI database Chl a concentration [mg m⁻³], the downwelling diffuse attenuation coefficient (k_d) at 490 nm [m⁻¹] and the absorption coefficient of phytoplankton (a_ph) at 443 nm [m⁻¹] were obtained from the merged information of four ocean colour sensors (SeaWiFS, MODIS-Aqua, MERIS, and VIIRS) at a spatial resolution of 4 × 4 km. The data was firstly filtered to exclude possible outliers related to the weather conditions and edge of cloud effects. This was performed by log transforming the data and then screening for outliers, which differed by more than four standard deviations (over 21 × 21 pixels area) from the average at each match-up station. Additional filtering was conducted according to the type of water masses; only Case 1 waters, according to [39], were analysed in the study. All pixels classified as eutrophic based on a euphotic depth of <9.8 m [40], and for Case 2 waters based on k_d 490 nm > 0.47 m⁻¹, were removed.

Daily and eight-day PAR (E m⁻² day⁻¹) data with a spatial resolution of 4 × 4 km were downloaded from the NASA’s OceanColor Web. Averaged values from two ocean colour sensors, SeaWiFS and MODIS-Aqua, were used. Preliminary analysis showed that the difference between PAR values from two sensors was less than 3%. Therefore, using one sensor only, instead of the average of the two, when another sensor is not available does not cause significant errors.

Daily SST data [°C] with a spatial distribution 1 × 1 km were downloaded from the Multiscale Ultra-high Resolution Sea Surface Temperature database, which is a merged data set from AVHRR, MODIS-Aqua, and AMSR, in which the influence of clouds is reduced due to the use of both infrared and microwave satellite data. As daily SST data have a 1 × 1 km spatial distribution, we re-gridded them to 4 × 4 km, in alignment with the ocean colour data.

2.4. Satellite Models of Primary Production

The models that were assessed in this study are commonly used to estimate PP, but are architecturally different. PP_VGPM is forced by ocean colour Chl a and the variability in the photosynthetic rate is parameterized from satellite SST. PP_Aph is based on the absorption properties of phytoplankton (a_ph), which can be derived directly from the inversion of remote sensing reflectance rather than from Chl a [41]. The principal advantage of using phytoplankton absorption at 443 nm-a_ph (443) is that it minimizes the effects of coloured dissolved organic matter and non-algal absorption, which can adversely affect PP models by modifying the light field [41,42]. PP_PSM is more complex than the other two models [43], since it simulates changes in photosynthesis as a function of irradiance using a two-parameter exponential P-I function based on the theory of initial light absorption by the photosynthetic system [30,31]. By contrast, PP_VGPM and PP_Aph use the Michaelis-Menten equation [44].

PP_VGPM retrieves PP_eu. It was parameterized using Chl₉₀ [mg m⁻³], which is the average Chl a concentration in a layer of the first optical depth seen by ocean colour sensors, z_eu, which is the euphotic depth, [m], I₀ is the daily surface PAR, [E m⁻² day⁻¹], DL is the day length, [h], and an empirical coefficient, 0.66125 (see Equation (1)). P^B_opt [mg C (mg Chl)⁻¹ h⁻¹], in the VGPM, was retrieved as a function of SST according to the algorithm presented in [6]. The P-I function (given in square brackets) is based on the Michaelis-Menten equation, with the half saturated constant equal to 4.1 E m⁻² day⁻¹:

$$P{P}_{eu}=Ch{l}_{90}{z}_{eu}0.66125\left[P_{opt}^B\frac{I_0}{4.1+I_0}\right]DL,$$

(1)

PP_PSM retrieves daily PP for a specific depth, z, PP(z), based on the P-I function as an exponential curve:

$$PP\left(z\right)=Chl\left(z\right)\left[P_m^B\left(1-\mathit{exp}\left(-\frac{{\alpha }^{B}I\left(z\right)}{P_m^B}\right)\right)\right]DL,$$

(2)

where, Chl(z) is Chl a at a specific depth, z, (mg m⁻³) and I(z) is PAR at a depth, z. The P-I function is based on two independent parameters: the biomass specific maximum rate of photosynthesis (P^B_m) (mg C (mg Chl)⁻¹ h⁻¹) and the initial slope (α^B) (mg C (mg Chl)⁻¹ h⁻¹ (μ E m⁻² s⁻¹)⁻¹).

For PP_Aph, daily PP(z) is calculated as follows:

$$PP\left(z\right)=a_{ph}443\times \left[{\varphi }_m\times \frac{K_{\varphi }}{K_{\varphi }+I\left(z\right)}\right]\times I\left(z\right),$$

(3)

where a_ph (443) is the coefficient of light absorption by phytoplankton at 443 nm [m⁻¹]. In Equation (3), the expression in square brackets corresponds to the photosynthesis dependence on underwater illumination expressed in terms of quantum yield (φ), which is the number of molecules of carbon dioxide assimilated by phytoplankton per absorbed light quantum [45]. In PP_Aph, the φ-I curve is also presented as the Michaelis-Menten equation, where $\varphi =\frac{P}{I}$ [46]. For the φ-I curve, φ_m is the maximum quantum yield of photosynthesis, [mg C E⁻¹], which is observed at low illumination, and K_φ is the half-saturation constant, [E m⁻² day⁻¹]. Usually, the value of K_φ is taken equal to 10 E m⁻² day⁻¹ [28,32,33,42,47], corresponding to φ_m = 0.06 mole C E⁻¹. These values of φ_m and K_φ were obtained based on the analysis of experimental data; the value of φ_m corresponds to a theoretical maximum [44,48,49].

For PP_PSM and PP_Aph, the photosynthetic parameters were chosen to be as close as possible to those observed in the NEA. In situ values of P^B_m, α^B, and φ_m were taken during different seasons (spring, summer, and autumn) and at different depths [32,43,50,51,52,53,54,55] as average values and used in the models as constant values (

Table 2. Average and range of photosynthetic parameters used in PP_PSM and PP_Aph: α^B [mg C (mg Chl)⁻¹ h⁻¹ (μ E m⁻² s⁻¹)⁻¹], P^B_m [mg C (mg Chl)⁻¹ h⁻¹], and φ_m [mole C E⁻¹].

	αB	PBm	φm
N	12	29	13
mean	0.049	3.316	0.032
SD	0.019	2.153	0.016
min	0.017	0.947	0.010
max	0.078	9.136	0.060

). For PP_Aph, K_φ was determined from the φ-I curve for the widespread coastal phytoplankton diatom species, Skeletonema costatum, following [56]. Thus, for a value of φ_m of 0.032 mole C E⁻¹ (or 347 mg C E⁻¹ used for PP_eu retrieval), K_φ is equal to 138.6 μ E m⁻² s⁻¹. Since, the φ-I function in PP_Aph is expressed as daily values, K_φ is between 4.5–8.0 E m⁻² day⁻¹ depending on the day length.

To facilitate the comparison of PP_PSM, PP_Aph, and PP_VGPM, the P-I functions of PP_PSM and PP_Aph were integrated over the depth of the euphotic layer (presented in Appendix A,

Table A1. Mathematical equations for PP_PSM and PP_Aph after integrating their P-I functions over the euphotic depth; (a) E₁ is an exponential integral defined as: $E_1\left(x\right)={\displaystyle \underset{x}{\overset{\infty }{\int }}}\frac{e^{-t}}{t}dt$; (b) mathematical representation of photoinhibition in the models tested.

(a) Mathematical Equations for the Models
Model
PPβPSM	$P{P}_{eu}=\frac{Ch{l}_{90}\times DL\times P_m^B}{k_d}\left[\right(E_1\left(\frac{I_0{\beta }^B}{P_m^B}exp\left(-k_d{z}_{eu}\right)\right)-E_1(\frac{I_0{\beta }^B}{P_m^B}))-(E_1\left(\frac{I_0\left({\alpha }^B+{\beta }^B\right)}{P_m^B}exp\left(-k_d{z}_{eu}\right)\right)-E_1(\frac{I_0\left({\alpha }^B+{\beta }^B\right)}{P_m^B})\left)\right]$
PPPSM	$P{P}_{eu}=\frac{Ch{l}_{90}\times DL\times P_m^B}{k_d}\times \left[E_1\left(\frac{{\alpha }^B\times I_0}{P_m^B}\right)+\left(k_d\times z_{eu}\right)-E_1\left(\frac{{\alpha }^B\times I_0}{P_m^B}\times exp\left(-k_d\times z_{eu}\right)\right)\right]$
PPβAph	$P{P}_{eu}=\frac{aph443\times {\varphi }_m\times K_{\varphi }}{k_d}\times exp\left(\beta K_{\varphi }\right)\times \times \left[E_1\left(\beta \left(K_{\varphi }+I_{0}exp\left(-k_d{z}_{eu}\right)\right)\right)-E_1\left(\beta \left(K_{\varphi }+I_0\right)\right)\right]$
PPAph	$P{P}_{eu}=\frac{aph443\times {\varphi }_m\times K_{\varphi }}{k_d}\times [\ln \left(K_{\varphi }+I_0\right)-\ln \left(K_{\varphi }+I_0\times \exp \left(-k_d\times z_{eu}\right)\right)]$
(b) Mathematical Representation of Photoinhibition
Model	Photoinhibition Term	Reference
PPPSM	$exp\left(-\frac{{\beta }^B}{P_m^B}I\left(z\right)\right)$	[43]
PPAph	$exp\left(-\beta \times I\left(z\right)\right)$	[32]
PPVGPM	$\begin{array}{c}exp\left(-\beta \times \overline{I}\right)\\ \overline{I}=\frac{{\int }_0^{z_{eu}}I\left(z\right)dz}{z_{eu}}=0.215I_0\end{array}$	this paper

(a). For each model, it was assumed that the rate of photosynthesis does not depend on a wavelength of the electromagnetic spectrum and is constant during a photoperiod. The Beer–Lambert–Bouguer law, as a function of the surface PAR (I₀), was used to propagate I(z) over depth, as follows:

$$I\left(z\right)=I_{0}exp\left(-k_{d}z\right)$$

(4)

where k_d is the downwelling diffuse attenuation coefficient of solar radiation in water [m⁻¹]. The euphotic depth was calculated using the following equation, which is derived from Equation (4), assuming that the surface radiation corresponds to 100% and radiation at the euphotic depth is 1%:

$$z_{eu}=\frac{ln\left(100\right)}{k_d}=\frac{4.6}{k_d}$$

(5)

where k_d for the whole spectrum of PAR was calculated as a function of the satellite product, k_d, at 490 nm according to the algorithm given in [57].

2.5. Implementation of Photoinhibition in the Primary Production Models

PP can be modified by photoinhibition because of phytoplankton being subjected to high irradiance or from moving from low to high light conditions. A photoinhibition function was also implemented for each model to test the significance of this phenomenon in the North Atlantic. The mathematical parametrization of photoinhibition is given in Appendix A, Table A1(b), where β^B and β are the photoinhibition parameters, and the superscript, ^B, represents the photoinhibition normalized to biomass.

PP^β_PSM uses β^B [mg C (mg Chl)⁻¹ h⁻¹ (μ E m⁻² s ⁻¹)⁻¹], whereas PP^β_Aph uses β [(E m⁻² day⁻¹)⁻¹]. Photoinhibition in PP^β_VGPM is synonymous with the photoinhibition term in PP^β_Aph, where $\overline{I}$ is the weighted average irradiance in the water column [E m⁻² day⁻¹]. Values of β^B = 0.01 mg C (mg Chl)⁻¹ h⁻¹ (μ E m⁻² s⁻¹)⁻¹ [43] were used for PP^β_PSM, and β = 0.01 (E m⁻² day⁻¹)⁻¹ [32] for PP^β_Aph and PP^β_VGPM were used as constants at all stations. Figure 2 shows the P-I functions of the models used. For PP^β_Aph, the φ-I function is multiplied by I.

2.6. Validation Statistics

Satellite derived estimates of PP_eu for each 4 × 4 km pixels were compared with ship-borne in situ measurements using the statistics described in [58]. Satellite PP_eu data was extracted from 3 × 3 pixels and averaged over the area to compare with in situ point data. For the daily data, only satellite data with a coefficient of variation (CV) of <0.15 over the 3 × 3 pixels were used. Since the OC-CCI data is from the merger of four Ocean Colour sensors (SeaWiFS, MODIS-Aqua, MERIS, and VIIRS), each of these satellites have different over pass times with respect to the in situ data. The Atlantic Meridional Transect field campaigns were designed to match the overpass of these sensors within ±3 h of the different satellite overpasses. Firstly, Taylor [59] and scatter diagrams were used to give the pair correlation (r), Standard Deviation (SD), centre-pattern Root Mean Square Error (c.-p. RMSE), and linear regression (coefficients S and I) between the modelled and in situ PP_eu data. For small rows (less than 28 values), the Spearman r was calculated, while for longer ones the Pearson r was used. In addition, the bias, RMSE, and Absolute Percentage Deviation (APD) were calculated. A one-way analysis of variance (ANOVA) was used to assess whether there were significant differences between the model and in situ PP_eu means.

Statistics were calculated for all stations in individual provinces and by season except for the ARCT province, which were only available during summer.

For the accuracy assessment, results are only included in

Table A3. Validation statistics for each model for the whole region and for individual biogeographical provinces using eight day composite OC-CCI data. The numbers at the head of the table indicate statistics for one test. The most significant statistical values are highlighted in grey. r is the Pearson or Spearman pair correlation coefficient (For N > 28, Pearson was used; for N < 28 Spearman was used, but r² was not computed), r² is the coefficient of determination, SD is the Standard Deviation, RMSE is the Root Mean Square Error, S is the Slope, I is the Intercept of the linear regression coefficients, APD is the Absolute Percentage Deviation, and F and p are ANOVA (Analysis of variance) coefficients. Those models, for which the significance of r was more than 0.1, were not included into the analysis and are not presented in the Table for summer stations in NAST (N = 16) and ARCT (N = 12), autumn stations in NADR (N = 6) and NAST (N = 24).

All Stations, N = 95	1	2		3		4	5	6	7	8
	SD	r	r2	S	I	bias	RMSE	Centre-Pattern RMSE	APD, %	ANOVA, F	ANOVA, p
In situ PPeu	0.68
PPβVGPM	0.81	0.74	0.54	0.88	0.77	0.13	0.57	0.56	9	1.4	0.237
PPVGPM	0.81	0.74	0.55	0.88	0.87	0.21	0.59	0.55	9	3.9	0.051
PPβPSM	0.90	0.76	0.58	1.01	−0.85	−0.78	0.97	0.58	15	45.6	0.000
PPPSM	0.84	0.78	0.61	0.97	−0.02	−0.21	0.56	0.52	8	3.6	0.060
PPβAph	0.66	0.71	0.51	0.69	2.23	0.51	0.71	0.50	11	27.0	0.000
PPAph	0.65	0.72	0.52	0.69	2.38	0.64	0.81	0.50	13	44.2	0.000
NADR, N = 28
In situ PPeu	0.54
PPβVGPM	0.35	0.73	0.53	0.47	3.70	0.48	0.61	0.38	9	15.1	0.000
PPVGPM	0.35	0.74	0.54	0.48	3.69	0.55	0.66	0.37	10	19.5	0.000
PPβPSM	0.43	0.79	0.62	0.62	1.91	−0.40	0.52	0.34	7	9.1	0.004
PPPSM	0.43	0.82	0.67	0.64	2.25	0.10	0.33	0.32	5	0.5	0.468
PPβAph	0.37	0.80	0.63	0.54	3.18	0.41	0.53	0.34	8	10.4	0.002
PPAph	0.38	0.80	0.65	0.56	3.18	0.52	0.62	0.33	9	16.7	0.000
NAST, N = 40
In situ PPeu	0.57
PPβVGPM	0.55	0.52	0.27	0.51	2.59	0.07	0.55	0.55	9	0.32	0.572
PPVGPM	0.56	0.54	0.29	0.53	2.57	0.17	0.57	0.54	10	1.68	0.199
PPβPSM	0.49	0.51	0.26	0.44	1.88	−0.97	1.10	0.52	19	65.16	0.000
PPPSM	0.50	0.55	0.30	0.48	2.29	−0.34	0.61	0.51	10	7.76	0.008
PPβAph	0.42	0.42	0.17	0.31	4.06	0.55	0.78	0.55	14	23.80	0.000
PPAph	0.44	0.43	0.19	0.34	4.06	0.70	0.89	0.55	16	36.57	0.000
NATR, N = 14
In situ PPeu	0.33
PPβVGPM	0.66	0.51		1.49	−3.22	−0.49	0.68	0.46	11	5.78	0.024
PPVGPM	0.66	0.51		1.49	−3.13	−0.39	0.61	0.47	10	3.58	0.069
PPβPSM	0.62	0.57		1.44	−3.89	−1.44	1.50	0.42	26	54.11	0.000
PPPSM	0.63	0.54		1.44	−3.22	−0.77	0.88	0.43	14	15.14	0.001
All Summer Stations, N = 56
In situ PPeu	0.49
PPβVGPM	0.68	0.73	0.54	1.01	0.07	0.14	0.48	0.46	7	1.5	0.218
PPVGPM	0.66	0.73	0.53	0.98	0.35	0.23	0.51	0.45	7	4.3	0.041
PPβPSM	0.80	0.76	0.58	1.25	−2.20	−0.74	0.91	0.54	13	33.9	0.000
PPPSM	0.72	0.76	0.58	1.11	−0.84	−0.17	0.50	0.47	7	2.1	0.146
PPβAph	0.55	0.67	0.45	0.76	1.94	0.50	0.66	0.43	9	24.9	0.000
PPAph	0.53	0.66	0.44	0.72	2.32	0.64	0.77	0.42	11	43.0	0.000
ARCT, N = 12
In situ PPeu	0.40
PPβPSM	0.28	0.56		0.45	3.01	−0.37	0.47	0.30	7	6.3	0.020
PPPSM	0.28	0.55		0.46	3.37	0.09	0.31	0.30	4	0.4	0.547
NADR, N = 22
In situ PPeu	0.39
PPβVGPM	0.30	0.57		0.40	4.09	0.39	0.52	0.34	7	13.3	0.001
PPVGPM	0.29	0.55		0.40	4.22	0.46	0.58	0.34	8	18.9	0.000
PPβPSM	0.34	0.53		0.53	2.45	−0.47	0.57	0.33	8	16.9	0.000
PPPSM	0.32	0.59		0.50	3.14	0.05	0.32	0.31	5	0.2	0.669
PPβAph	0.28	0.58		0.46	3.68	0.33	0.45	0.30	6	9.7	0.003
PPAph	0.28	0.62		0.45	3.88	0.45	0.54	0.30	8	18.6	0.000
NATR, N = 6
In situ PPeu	0.19
PPβVGPM	0.21	0.81		0.80	0.57	−0.50	0.52	0.15	9	15.4	0.003
PPVGPM	0.21	0.81		0.76	0.93	−0.39	0.42	0.16	7	9.2	0.012
All Autumn Stations, N = 38
In situ PPeu	0.55
PPβVGPM	0.68	0.41	0.16	0.50	2.64	0.10	0.69	0.68	12	0.5	0.477
PPVGPM	0.67	0.41	0.17	0.50	2.70	0.18	0.69	0.67	12	1.6	0.213
PPβPSM	0.65	0.47	0.22	0.54	1.42	−0.87	1.07	0.62	18	38.6	0.000
PPPSM	0.61	0.50	0.25	0.55	1.97	−0.28	0.65	0.58	11	4.4	0.039
NATR, N = 8
In situ PPeu	0.38
PPβVGPM	0.84	0.69		1.67	−4.29	−0.48	0.77	0.60	12	1.9	0.190
PPVGPM	0.85	0.69		1.69	−4.29	−0.39	0.72	0.60	12	1.2	0.286
PPβPSM	0.77	0.64		1.53	−4.38	−1.38	1.48	0.53	25	18.0	0.001
PPPSM	0.78	0.69		1.58	−4.00	−0.72	0.90	0.54	13	4.8	0.046

(Appendix B) when the r between the modelled and in situ PP_eu have p-levels of <0.1, 0.05, and 0.01. For all statistics, the highest significant values are highlighted in grey to indicate the most accurate model.

2.7. Sensitivity Analysis

To estimate the contribution of each input parameter to the variability in satellite derived PP_eu, a sensitivity analysis was conducted following [14]. The average values of the satellite derived input parameters and the P-I functions were sequentially varied over their natural range.

Firstly, “reference” values of PP_eu were calculated based on averaged values of each of the parameters for the three models without photoinhibition. Then, each parameter was varied sequentially, whilst keeping the other parameters fixed at their average value. The results are presented in the form of box-and-whisker diagrams. Since, z_eu was estimated by k_d (see Equation (5)), it was not included in the sensitivity analysis, but k_d was considered instead.

All calculations and visualization of the results were carried out using Matlab R15b.

3. Results

3.1. Accuracy Assessment of Primary Production Models

Each model was firstly run using daily OC-CCI data and compared with in situ PP_eu at 46 stations in the NEA (Figure 3). Since the OC-CCI data is available as eight day composites, a secondary accuracy assessment was carried out using these data at the same stations (Appendix B,

Table A2. Validation statistics for each model using (a) daily and (b) eight day composite OC-CCI data (N = 46). The numbers at the head of the table indicate statistics for one test. The most significant statistical values are highlighted in grey. r is the Pearson pair correlation coefficient, S is the Slope and I is the Intercept of the linear regression coefficients, RMSE is the Root Mean Square Error, and APD is the Absolute Percentage Deviation.

(a) Daily, N = 46	1		2	3		4	5	6
	r	p-Level	Centre-Pattern RMSE	S	I	Bias	RMSE	APD, %
PPβVGPM	0.63	p < 0.001	0.60	0.68	1.73	−0.04	0.60	9.2
PPVGPM	0.63	p < 0.001	0.60	0.68	1.80	0.05	0.61	9.3
PPβPSM	0.69	p < 0.001	0.57	0.79	0.16	−0.99	1.15	18.8
PPPSM	0.68	p < 0.001	0.57	0.77	0.89	−0.38	0.68	9.9
PPβAph	0.71	p < 0.001	0.48	0.56	2.86	0.42	0.63	9.8
PPAph	0.70	p < 0.001	0.48	0.57	2.95	0.56	0.74	11.8
(b) Eight day composites, N = 46
PPβVGPM	0.65	p < 0.001	0.56	0.66	1.80	−0.06	0.57	8.6
PPVGPM	0.66	p < 0.001	0.56	0.67	1.86	0.03	0.56	8.6
PPβPSM	0.69	p < 0.001	0.57	0.79	0.16	−0.99	1.15	18.8
PPPSM	0.73	p < 0.001	0.49	0.75	0.98	−0.40	0.64	9.5
PPβAph	0.75	p < 0.001	0.45	0.53	3.00	0.40	0.60	9.4
PPAph	0.73	p < 0.001	0.46	0.54	3.09	0.53	0.70	11.4

) and then including more data (Appendix B, Table A3). Using N = 46, both the daily and eight day modelled data of PP_eu have almost identical results.

Using the eight day OC-CCI composite data, PP_PSM was the most accurate model over all stations (Figure 4e). When using the photoinhibition model, PP^β_VGPM tended to over-estimate in situ PP_eu at higher values and under-estimate it at lower values (Figure 4a). The net result, however, was that, compared to the other models with photoinhibition, the regression line for PP^β_VGPM was closest to the 1:1 line and the intercept was the smallest. By comparison, PP^β_Aph tended to over-estimate in situ PP_eu, especially at low values, which increased the scatter, resulting in a higher intercept (Figure 4c). For PP^β_PSM, all points were below the 1:1 line, indicating an under-estimate compared to in situ PP_eu (Figure 4b). The scatter was low, however, which increased the percentage variance explained and, though the slope was off set, it was highly correlated with the 1:1. Using no photoinhibition, the trend was the same for PP_VGPM and PP_Aph (Figure 4d,e), with similar statistical results to the models implemented with photoinhibition. For PP_PSM, there was a significant improvement in the model prediction of in situ PP_eu, with the regression line close to the 1:1 line and the intercept was reduced.

A Taylor diagram was used to illustrate the statistical relationships between modelled and in situ PP_eu (Figure 5). All models have a high correlation with in situ PP_eu, with r in the segment from 0.7–0.8. The modelled values, however, have a large c.-p. RMSE, though it did not exceed the SD of the in situ measurements. The smallest c.-p. RMSE was for PP_Aph, PP^β_Aph, and PP_PSM. PP_PSM had the largest r, but its SD exceeded that of the in situ PP_eu, whereas the SD of PP_Aph and PP^β_Aph were closer to the in situ PP_eu SD. Photoinhibition in PP_PSM significantly increased the c.-p. RMSE and created a difference in SD compared to the other models.

Table A3 (Appendix B) gives the statistical results between the modelled and in situ PP_eu. For all stations, PP_PSM was the most accurate in four out of eight statistical tests. It slightly under-estimated PP_eu (bias = −0.21) and explained 61% of the variance in in situ measurements (coefficient of determination (r²) = 0.61). There was an 8% Absolute Percentage Difference (APD) between PP_PSM and in situ PP_eu, and this model also had the smallest RMSE. PP_PSM exhibited a better performance in summer when it explained 58% of the variance in in situ PP_eu, compared to autumn when it was 25%. In addition, PP_PSM did not always accurately capture the spatial variability in in situ PP_eu since the SD was 1.5 times higher than the SD of in situ PP_eu, except for the ARCT and NADR provinces in summer.

By comparison, PP^β_VGPM was the most accurate in just two out of eight statistical tests, though it was worse in summer (in five out of eight statistical tests). Both PP^β_VGPM and PP_VGPM slightly over-estimated in situ PP_eu (bias = 0.13–0.21) except in the NATR province, where it underestimated in situ values. On average, PP^β_VGPM and PP_VGPM diffed from in situ PP_eu by 8% and 9% (APD), respectively. For all stations, r² for PP^β_VGPM and PP_VGPM explained 54 to 55% of the variance of in situ PP_eu. The PP^β_VGPM and PP_VGPM was the most accurate in the NATR province and during summer, and the statistical errors were small.

Both PP^β_Aph and PP_Aph over-estimated in situ PP_eu (bias = 0.51–0.64) compared to the other models, especially at low values (Appendix B, Table A3). APD for PP^β_Aph and PP_Aph were 11% and 13%, respectively. PP^β_Aph and PP_Aph, both captured the spatial variability in in situ PP_eu for all data, and the SD was similar to that of in situ PP_eu. In addition, the c.-p. RMSE for PP^β_Aph and PP_Aph was small compared to the other models, though they were less accurate in autumn compared to summer.

3.2. Satellite Primary Production Images

The spatial and temporal differences between the models are illustrated in satellite images during two different periods: in summer (average June–August, 1998–2011; Figure 6a–c) and autumn (average September–November, 1998–2011; Figure 6d–f). During summer in the NADR, the PP_VGPM gave the highest values (500–1500 mg C m⁻² d⁻¹), whereas the PP_PSM was 250–700 mg C m⁻² d⁻¹ and PP_Aph was 500–1000 mg C m⁻² d⁻¹. By comparison, in the NADR autumn, the range in PP_VGPM was 300–700 mg C m⁻² d⁻¹, PP_PSM was 100–300 mg C m⁻² d⁻¹, and PP_Aph was 250–500 mg C m⁻² d⁻¹. In contrast, in the NAST in summer, PP_Aph gave the highest values (250–700 mg C m⁻² d⁻¹), which were lower in autumn (200–350 mg C m⁻² d⁻¹). The PP_VGPM returned similar values in summer and autumn, which were 125–500 mg C m⁻² d⁻¹, and PP_PSM was 0–350 mg C m⁻² d⁻¹. Similarly, in the NATR, PP_Aph also returned the highest values in both summer and autumn, PP_PSM yielded the lowest (<125 mg C m⁻² d⁻¹), and PP_VGPM were 50–250 mg C m⁻² d⁻¹.

To illustrate further the spatial and temporal differences between the models, we extracted data from each model at every 4 km along the transect at 20° W from 20° to 60° N (Figure 7). In the NADR, there was a 182% and 157% Relative Percentage Deviation (RPD) between PP_VGPM and PP_PSM in summer and autumn, respectively, and 132% and 113% RPD between PP_Aph and PP_PSM, respectively. In the NAST, the differences increased to a 228% and 203% difference between PP_VGPM and PP_PSM and a 359% and 351% for PP_Aph and PP_PSM in summer and autumn, respectively. For the NATR, there was a 222% and 197% RPD between PP_VGPM and PP_PSM and a 208% and 224% for PP_Aph and PP_PSM in summer and autumn, respectively.

3.3. Model Sensitivity Analysis

For all three models, the largest contribution to the variability in PP was the remote sensing ocean colour product used, which for PP_VGPM and PP_PSM was Chl a, and for PP_Aph was a_ph (443) (Figure 8). The next most important parameter depends on the model used. For PP_VGPM, k_d for these waters caused the greatest error. For PP_PSM, both P^B_m and k_d produced the greatest error and for PP_Aph it was φ_m. The third most important parameter for the PP_VGPM was P^B_opt and DL, for PP_PSM it was I₀ and DL, and for PP_Aph it was k_d and I₀. The smallest contribution to the variability in PP_VGPM was from I₀, but for PP_PSM it was α^B, and for PP_Aph it was K_φ.

4. Discussion

4.1. Validation of Primary Production Models

Validation of PP_VGPM in the North East Atlantic Ocean has already been carried out by a number of independent studies in this and neighbouring seas [7,9,10,12,13,14,60,61]. These studies report that PP_VGPM over-estimates in situ PP_eu, which we also observed in the NEA, especially in NADR. The only exception is in the NATR province, where PP_VGPM is lower than in situ PP_eu. As previously highlighted [7,14], the use of P^B_opt to simulate model variations in photosynthesis as a function of SST may not consider the displacement of the layer of optimum light conditions for photosynthesis to lower depths, which can be caused by high surface PAR [62].

According to [10], PP_VGPM is closer to the actual global PP than PP_PSM, which tends to under-estimate global PP by approximately 13%. In our study, PP_VGPM are slightly higher than PP_PSM in summer and autumn in the Atlantic Ocean, especially in the NADR, NATR, and NAST (Appendix B, Table A3). Only in the ARCT province, where the waters are typically mesotrophic with a medium range Chl a [36], did both models converge.

Validation of PP_Aph in shelf waters off Argentina [28] showed a high correlation with in situ PP_eu (r = 0.9), but there was a tendency for the model to under-estimate PP_eu, which is due to the use of the constant, K_φ = 10 E m⁻² day⁻¹. In our study, we allowed K_φ to vary between 4.5–8.0 E m⁻² day⁻¹ depending on the day length. PP_Aph slightly overestimated in situ PP_eu, especially in autumn, which was also observed by [28] in Argentinian shelf waters.

PP_Aph in the open ocean waters of the South West Atlantic is also more accurate than PP_VGPM ([33]; PP_Aph r² = 0.74, PP_VGPM r² = 0.44). Theoretically, this is due to the advantage of deriving a_ph (443) directly from R_rs rather than using an empirical band ratio to derive Chl a [63]. However, this depends on how accurately both a_ph (443), a_CDOM (443), and a_det (443) can be determined from R_rs (443). If there is an error in deriving a_CDOM (443) and a_det (443) it will propagate to a_ph (443), which can, in turn, impact PP_Aph. Based on the sensitivity analysis, a_ph (443) causes a greater error in PP_eu compared to Chl a, which is used in PP_PSM (Figure 8). Also for PP_Aph, an accurate estimate of the photosynthetic parameters (φ_m and K_φ) is also vital in improving the performance of this model [42]. We also observed that φ_m has a large impact on PP_Aph, but K_φ has the least effect (Figure 8). In our study, PP_Aph was more accurate in the NADR province (r² = 0.65), though it still over-estimated in situ PP_eu and was worse in the ARCT province in summer (Appendix B, Table A3), though this was based on few points. Furthermore, Figure 9 shows that a_ph (443) causes a significantly greater error in PP_Aph in the ARCT compared to Chl a used in the PP_VGPM and PP_PSM.

Standard satellite validation protocols were used to assess the accuracy of the PP models. The satellite data were averaged over 3 × 3 pixels with CV < 0.15 for the majority of the match-ups (Figure 10), which illustrates that the quality of the satellite data was good. For the daily data, the CV for the three models was always <0.15. For the eight day composites, there were two data points for PP_VGPM and PP_PSM where the CV > 0.15 and six points for PP_Aph, which were mainly from the ARTC province.

We found that all models performed better in summer than in autumn. Since Chl a is the main contributor to the overall error in modelled PP_eu, it may suggest that the error in the autumn values arise from an error in the satellite Chl a. However, numerous studies have shown for the North Atlantic that satellite Chl a is accurate during this season (e.g., [14]).

According to the sensitivity analysis, P^B_m and φ_m are the second parameters that determine the overall error of PP_eu retrieval. The values of P^B_m and α^B used in the study for the NEA (Table 2) correspond with values from a global database of in situ photosynthesis parameters presented in [64]. Mean values of P^B_m and α^B averaged over four provinces (ARCT, NADR, NAST, and NATR) and over spring, summer, and autumn from this database are 3.25 mg C (mg Ch)⁻¹ h⁻¹ and 0.042 mg C (mg Chl)⁻¹ h⁻¹ (μ E m⁻² s⁻¹)⁻¹, respectively. A value of φ_m = 0.03 mol C E⁻¹ used in our study is characteristic of sunlit regions and was also used by [9] for retrieving PP_eu in different parts of the global ocean. An identical value of φ_m was obtained in waters off the southern California coast (30–34° N) [65]. Errors in PP_PSM may arise from inaccurate derivation of the photosynthetic parameters, P^B_m and α^B, which can vary over space, time, and depth depending on the dominant phytoplankton community.

To calculate PP_PSM, we used 29 values of P^B_m, and 20 values during autumn. When a P^B_m value, typical for autumn, of 2.36 mg C (mg Ch)⁻¹ h⁻¹ was used, it did not improve the correlation between in situ and modelled PP_eu (r² = 0.23). Spatial changes in the species composition of phytoplankton communities and the related changes in photosynthetic rates may also cause a weak correlation between the modelled and in situ PP_eu [66]. In addition, we input the Chl a and photosynthetic parameters as constant values over depth, which may not be representative of the vertical variability of these parameters in autumn and may therefore contribute to the error in the model. Kahru et al. [42] validated the PP_Aph algorithm in the California Current, and showed that accounting for the vertical profile of a_ph (443) within the euphotic zone improved the correlation between PP_Aph and in situ PP_eu. The coefficient of determination (r²) increased from 0.28 to 0.56. Errors in R_rs, which is used to calculate a_ph (443), could be the main cause of error in the retrieval of PP_eu. These estimates can vary widely between radiometers (15–18% for R_rs at 443 nm).

The tendency in the scatter plots is that PP_VGPM over-estimates PP_eu in summer and under- and over-estimates PP_eu in autumn, PP_Aph over-estimates PP_eu in summer, especially in autumn, and that PP_PSM shows good agreement in summer, but has a tendency to under-estimate PP_eu in autumn (Figure 4). Visual examination of the resulting mapped images provides a qualitative analysis of whether the tendencies shown in the scatter plots correspond over wider spatial and temporal scales. The extent of the differences between the models in the mapped images and data extracted from these over a transect from 20 to 60° N illustrate that, of the three models, PP_VGPM gives the highest PP_eu in summer and autumn in the NATR and NADR, whereas the PP_Aph has the highest PP_eu in the NAST, especially in autumn (Figure 7). The PP_PSM consistently has the lowest PP_eu, but is closest to the in situ PP_eu, but the differences between PP_PSM and the other models are less in autumn. The outliers in the scatter plot for the PP_Aph model (Figure 4f) may suggest that, in the higher PP waters of the NADR, this model would return the highest PP_eu. The images and transect, however, show that PP_VGPM gives the highest PP_eu, which is not reflected in the corresponding scatter plot (Figure 4d). The tendency for the PP_VGPM to over-estimate PP_eu at higher values in the scatter plot represents large areas of the NADR and NATR during summer in the corresponding mapped image (Figure 6a).

The results suggest that the PP_PSM model with the OC-CCI data is accurate for the NEA. The model is available as an OC-CCI product, but requires further testing to assess its accuracy over other basins and the entire global ocean.

4.2. Effect of Photoinhibition in Primary Production Models

Photoinhibition describes the decrease in the photosynthetic rate with increasing irradiance [67]. This phenomenon is dependent on the phytoplankton community and its adaptation to the light conditions in the water column. It is usually greatest when phytoplankton is confined to surface high light conditions or when it moves quickly from low to high light [34]. Under such high light conditions, photorespiration can be activated, which is accompanied by utilization of assimilated carbon and, therefore, a decrease in the efficiency of photosynthesis [68]. For different phytoplankton communities, it has been shown experimentally that this effect varies as a function of depth and solar radiation at the surface and photo-adaptation of phytoplankton [69,70]. For example, shade-adapted communities saturate quickly and reach photoinhibition faster than high-light adapted communities.

Since the total amount of incoming solar radiation decreases with the distance from the equator and is dependent on the time of the year, the difference in the estimates of PP_eu, both with and without photoinhibition, also decrease with increasing latitude, and has a distinctive seasonal variation associated with the variability of incoming solar radiation during the year. The maximum difference in the estimates of PP_eu is in the second half of June, when solar radiation is at a maximum. The smallest difference is at the end of December or beginning of January, when solar radiation reaches a minimum. Since both states can occur in the tropical waters of the NAST and NATR provinces where weak wind stress and density gradient can prevent the migration of the phytoplankton community and, trap it at high light conditions for a significant proportion of the day, including photoinhibition may improve model performance. We found that the PP^β_VGPM, for example, performed better in the NAST during summer.

Table 3. Relative percentage difference (%) between modelled PP_eu with and without photoinhibition. The negative sign “−” indicates how much photoinhibition reduces PP_eu.

Model	Mean	Min	Max
PPVGPM	−7	−2	−10
PPPSM	−42	−20	−50
PPAph	−11	−4	−15

illustrates the mathematical representation of photoinhibition in the models used, which does not differ significantly between each model. Photoinhibition had a greater effect on PP_PSM compared to the other two models, however. Usually, photoinhibition in the models is presented as an exponentially decreasing term [70]. Photoinhibition in PP^β_PSM is implemented as the P-I curve increases at low light [43]. After integration over the euphotic depth, its mathematical representation can cause large differences between PP^β_VGPM and PP^β_Aph, even though the photoinhibition functions are identical. According to [32], the photoinhibition term in PP^β_Aph can reduce the values of the modelled PP_eu from 10% to 55% depending on the light intensity.

5. Conclusions

Three satellite models of PP were processed using OC-CCI data, with and without photoinhibition, and compared against ship-borne in situ ¹⁴C measurements from the NEA. Validation of the models using daily OC-CCI data (N = 46) showed that the PP_PSM and PP_VGPM had a similar accuracy, whereas the PP_Aph was not suitable for the region. Using eight day OC-CCI composite data, the number of match-ups was doubled (N = 95) and PP_PSM was more accurate over all provinces. For individual provinces, the PP_PSM was more accurate in the ARCT and NADR regions whereas the PP_VGPM was more accurate in the NAST and NATR regions. The PP_Aph was the least accurate model for the NEA, which was due to errors in the a_ph (443) OC-CCI product. The use of a photoinhibition function in the PP_VGPM and PP_Aph had little effect on the model performance, whereas, in the PP_PSM, photoinhibition reduced the accuracy of the model, especially at lower range of the values. The performance of the models also varied seasonally; all models were more accurate during summer and less accurate during autumn. Over all provinces, the PP_VGPM with photoinhibition was the most accurate during summer, and the PP_PSM with no photoinhibition was more accurate during autumn. Mapped images illustrated that differences between models were far greater than indicated by the scatter plots, especially for PP_Aph and PP_VGPM in the NATR and NADR provinces. A sensitivity analysis indicated that the primary ocean colour product (Chl a for PP_PSM and PP_VGPM; a_ph (443) for PP_Aph) caused the greatest variability in model performance, followed by the photosynthetic terms, P^B_m and φ_m, in the PP_PSM and PP_Aph models, respectively, and P^B_opt and k_d in the PP_VGPM. Improvement in the accuracy of these input parameters will ultimately lead to an improvement in satellite PP algorithms for this region.