Net primary productivity in the Gulf of California: spatio-temporal patterns and drivers in a warming climate (2003–2024)

The Gulf of California is a tropical–subtropical, productive and biodiverse marine ecosystem that supports nearly 70 % of Mexico's fisheries. Despite its ecological and socio-economic importance, the basin-scale patterns and drivers of seasonal and interannual spatial variability in net primary productivity (NPP) remain incompletely understood. Here, we analyse 22 years (2003–2024) of satellite-derived NPP data to identify the principal spatio-temporal patterns of NPP and their environmental controls. Monthly climatology shows that the cold season (December–May) contributes 68 % of annual NPP, whereas the warm season (June–November) accounts for only 32 %. To characterise regional variability, we applied the k-means${}^{++}$ clustering algorithm, selected the optimal number of clusters using the elbow and silhouette criteria, and validated the resulting classification through discriminant analysis. This procedure identified three distinct productivity zones: low-productivity (LPZ; < 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), mid-productivity (MPZ; 900–1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and high-productivity (HPZ; > 1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Interannually, the LPZ covered 36 % of the area but produces only 16.9 % of total NPP, whereas the HPZ contributes 61.4 % of total NPP while occupying only 30 % of the gulf's area. A pronounced NPP negative anomaly during 2013–2019 reduced NPP by approximately 9 % relative to the long-term mean (124 Tg C yr⁻¹). This decline coincided with sustained ocean warming, including multiple marine heatwaves and the strong 2015–2016 El Niño event. A partial recovery was observed during 2020–2024. Nitrate availability and water-column stratification were the dominant predictors of NPP, but their effects differed among productivity zones: nitrate was negatively associated with NPP and stratification index positively associated with NPP in the LPZ, whereas the opposite relationships were observed in the MPZ and HPZ. Wind forcing, sea-surface salinity and phosphate availability exerted additional zone-specific effects. Overall, our results indicate that anomalously warm conditions reduce the spatial extent of HPZ zones and favour the expansion of oligotrophic areas. Continued warming and enhanced stratification are therefore expected to increase the extent of LPZ, with potentially significant ecological consequences for the Gulf of California.

1 Introduction

The Gulf of California is a highly biodiverse ecosystem that supports approximately 70 % of Mexico's fisheries (Lluch-Cota et al., 2007). This gulf is influenced by multiple atmospheric and oceanic processes that complicate the characterisation of this complex marine system (Kahru et al., 2004; Bustos-Serrano and Castro-Valdez, 2006; López-Martínez et al., 2023; Marin-Enríquez et al., 2024). Among these processes, prolonged anomalous warming of the surface ocean exerts widespread negative impacts on marine ecosystems, with large-scale repercussions for ecological and biogeochemical functioning (Somavilla et al., 2017). Such events can alter particulate carbon flux (Xie et al., 2019), reduce biodiversity (Wernberg et al., 2016; Jones et al., 2018), and affect fisheries and aquaculture (Mills et al., 2013; Oliver et al., 2017). Enhanced ocean–atmosphere heat exchange during warming events strengthens water-column stratification (Jacox et al., 2016; Somavilla et al., 2017), disrupting nutrient supply and driving changes throughout marine food webs (Cavole et al., 2016). For instance, reductions in phytoplankton biomass and primary productivity have been documented during anomalous warming episodes such as El Niño events (Jacox et al., 2016).

Despite ongoing global warming, the responses of oceans and marine ecosystems to changing sea-temperature patterns remain uncertain (Somavilla et al., 2017). In recent decades, anomalous warming events, including the positive phase of El Niño-Southern Oscillation (ENSO) and the marine heatwaves, have increased in frequency across multiple ocean regions (Di Lorenzo and Mantua, 2016; Oliver et al., 2019). These events produce positive sea-surface temperature (SST) anomalies over a range of spatial and temporal scales and are generally associated with declines in primary productivity (Wernberg et al., 2016; Oliver et al., 2017). Concurrently, global ocean heat content has risen over the past five decades (NOAA, https://www.ncei.noaa.gov/access/global-ocean-heat-content/, last access: 1 October 2026), suggesting that marine ecosystems may be approaching resilience thresholds with potentially irreversible consequences (Frölicher et al., 2018).

The northeastern Pacific Ocean represents a transitional marine ecosystem (Fiedler and Talley, 2006; Cepeda-Morales et al., 2013), where warming events across diverse spatiotemporal scales have been reported in the Gulf of California (Sánchez-Cabeza et al., 2022; López-Martínez et al., 2023). Notably, the combined influence of the Blob in the North Pacific, strong equatorial El Niño event and overlapping marine heatwaves between 2013 and 2016 produced persistent warm anomalies in northwestern Mexico (Joh and Di Lorenzo, 2017; Jiménez-Quiroz et al., 2019; García-Fernández et al., 2023). These anomalies propagated through mid- and upper-trophic levels, altering species distributions, including euphausiids, copepods, pelagic crabs, and tunas, in the Eastern Tropical Pacific (Cavole et al., 2016; García-Fernández et al., 2023). As global temperatures rise, such anomalies are expected to disrupt the biological carbon pump, reducing oceanic carbon storage and increasing atmospheric carbon concentrations (Cavole et al., 2016; Xie et al., 2019). These changes will likely reshape community structures and ecosystem functioning (Cavole et al., 2016).

Figure 1 Bathymetric map of the Gulf of California, conceptually delimited from the Pacific Ocean by a red line. The gulf is a large and deep pull-apart basin formed by right-lateral strike-slip motion along faults associated with the San Andreas system. The resulting geological configuration supports oceanographic processes that sustain high phytoplankton productivity, with contrasting dynamics along the opposing coasts.

The Gulf of California (20–32° N, 105–115° W) is among the most diverse and productive marine regions in Mexico (Fig. 1), characterised by regions with high net primary productivity (NPP, defined as gross primary productivity minus respiration) and phytoplankton biomass expressed as chlorophyll a concentration (Chl a) (Lluch-Cota, 2004; O'Reilly and Sherman, 2016). Its elongate geometry, steep bathymetry, and pronounced orographic relief modulate regional pressure-gradient forcing and orographic channelling of the winds, generating seasonally coherent along-gulf winds whose wind-stress curl drives Ekman upwelling, coastal divergence, and vertical mixing that inject nutrient-rich subsurface waters into the euphotic zone (Miranda-Alvarez et al., 2020; Robles-Tamayo et al., 2020). However, anomalous warming associated with the positive phases of ENSO and Pacific Decadal Oscillation (PDO) reduce mean Chl a concentration in the gulf (Hakspiel-Segura et al., 2022; López-Martínez et al., 2023). In situ NPP measurements in the region remain scarce, particularly during winter, limiting our understanding of latitudinal gradients and seasonal cycles (Zeitzschel, 1969; Álvarez-Borrego and Lara-Lara, 1991; Mercado-Santana et al., 2017). Remote-sensing estimates of NPP and Chl a, based on Vertically Generalised Production Model (VGPM) incorporating light-use efficiency and absorbed photosynthetically active radiation, have partially addressed these gaps. However, most studies have focused on local variability rather than basin-scale-dynamics. Furthermore, while satellite-derived Chl a is often used as a proxy forphytoplankton biomass within the first optical depth, whereas NPP integrates production across the entire euphotic zone and daylight period (Balch et al., 1992; Behrenfeld and Falkowski, 1997).

In this study, we investigate the response of NPP in the Gulf of California to recent anomalous warming events, using a 2003–2024 satellite record to establish the seasonal and interannual baseline against which anomalies are assessed. We specifically focus on the 2013–2024 period, encompassing the 2014–2016 anomalously warming event, to evaluate how elevated SSTs and enhanced stratification modify the spatial distribution of productivity and its environmental controls. We hypothesize that anomalous warming events intensify water column stratification, restricting the flux of nutrient-rich water to the euphotic zone and thereby reducing NPP. By combining satellite observations with environmental reanalysis-based water-column variables, we identify the key drivers that modulate NPP under anomalously warm conditions.

2 Materials and methods

2.1 Study area

The Gulf of California (20–32° N, 105–115° W) is a marginal tropical–subtropical sea of the eastern Pacific located in northwestern Mexico (Fig. 1). It extends ∼ 1500 km from its northernmost limit to Bahía de Banderas in the south and is 150–200 km wide. South of 28° N, the gulf comprises a series of six deep basins separated by bathymetric sills (Álvarez-Borrego and Lara-Lara, 1991). The region exhibits dynamic and heterogeneous conditions driven by tidal currents, wind momentum transfer, bathymetrically controlled and wind-driven coastal upwelling, and intense solar heating (Álvarez-Borrego, 2007). In the Midriff Archipelago, bathymetric upwelling maintains low temperatures year-round. From December to May, north-westerly winds induce strong coastal upwelling along the eastern coast, transporting nutrient-rich subsurface waters into the euphotic zone, fuelling phytoplankton growth. In contrast, from July to October southeasterly winds are generally weaker and insufficient to disrupt the pronounced thermal stratification that develops during summer. June and November are considered transitional months between the cold-warm and warm-cold seasonal regimes (Santamaría-del-Ángel et al., 1994; Hidalgo-González and Álvarez-Borrego, 2004). The Baja California mountain ranges block low-lying clouds from the Pacific, resulting in predominantly clear skies over the gulf except during summer monsoon incursions (Badan-Dangon et al., 1985). For this study, the southern boundary of the Gulf of California was defined by a line connecting Cabo San Lucas on the Baja California peninsula (23° N) to Cabo Corrientes in Nayarit (∼ 20° N; Fig. 1).

2.2 Data collection

Net primary productivity (NPP) was obtained from the satellite-based VGPM product distributed through NOAA Environmental Research Division’s Data Access Program (ERDDAP). This product was selected because it provides spatially coherent NPP fields for the Gulf of California and reproduces the expected seasonal cycle of regional productivity, with enhanced winter-spring production and reduced summer productivity under stronger stratification. Reanalysis products were used only for subsurface and water-column variables, including satellite observations of mixed-layer depth, nutrients, salinity, temperature profiles, and the derived stratification index, for which long-term basin-scale observational datasets with consistent spatial and temporal coverage are not currently available in the Gulf of California. The following subsections summarize the datasets used for NPP, physical, biogeochemical, and climate variables.

2.2.1 Satellite variables

Monthly sea surface temperature (SST) and Chlorophyll a concentration (Chl a; OCI algorithm) were obtained from Aqua-MODIS Level 3 global composites (4 km × 4 km) spanning January 2003 to December 2024, downloaded from NASA Ocean Colour (https://oceancolor.gsfc.nasa.gov/l3/, last access: 1 October 2026).

Euphotic-zone depth (EZD) was derived from Chl a using the chlorophyll a dependent bio-optical parametrisation with piecewise power-law relationships for different trophic regimes (Morel and Berthon, 1989; Morel and Maritorena, 2001).

Net primary productivity (NPP) at the same nominal resolution (4 km × 4 km) was obtained from monthly Level-3 products of the Vertically Generalised Production Model (VGPM) via NOAA CoastWatch (https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1ppmday.graph, last access: 1 October 2026). Monthly 10-m wind components (u,v) were taken from Remote Sensing System (RSS) CCMP V3.1 Level-3 dataset (25 km × 25 km) for 2003–2024 (https://data.remss.com/ccmp/v03.1/, last access: 1 October 2026). To account for the gulf's orientation, the coordinate axes were rotated by 40°, aligning the total u-component with the gulf's axis and the rotated v-component perpendicular to it. The along-gulf wind u-component (URT) was used in statistical analyses.

Daily sea current components (u,v; Level 4; 12.5 km × 12.5 km) were obtained from the Copernicus Marine Environment Monitoring Service (CMEMS, https://doi.org/10.48670/moi-00148). Monthly means were computed for 2003–2024, and rotated analogously to the winds; only the along-gulf u-component (UAR) was used in statistical analyses.

Monthly global gap-free Level-4 (L4) sea surface salinity (SSS) analyses at 12.5 km spatial resolution were obtained from CMEMS (https://doi.org/10.48670/moi-00051) for the period 2003–2024 and included in the statistical analyses.

2.2.2 Reanalysis variables and derived parameters

Monthly Global Ocean Biogeochemistry Hindcast reanalysis (GLOBAL_REANALYSIS_BIO_001_029) for nitrate (NO₃), phosphate (${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$) and silicate (SiO₂) concentrations were downloaded from CMEMS (https://doi.org/10.48670/moi-00019). These Level 4 products are provided on a 25 km × 25 km grid with 75 depth levels. Mixed layer depth (MLD), potential temperature (θ), and practical salinity (PSU) were obtained from the Global Ocean Physics Reanalysis (GLOBAL_REANALYSIS_PHY_001_030; ∼ 8 km × 8 km horizontal resolution, 50 standard depth levels; https://doi.org/10.48670/moi-00021).

2.2.3 Data regridding

To ensure consistency, NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$, SiO₂, MLD, and the stratification index (Sect. 2.3.2) were bilinearly interpolated to the 4 km grid used for NPP, SST, and EZD. Only grid cells within the Gulf of California domain (Fig. 1) were retained for analyses.

2.3 Data processing

2.3.1 Identification of warm events, climate indices, and spatial trends

Anomalous warm periods were identified using monthly SST anomalies (SSTa), and marine heatwaves (MHWs) were detected using daily SST following Hobday et al. (2016) who define MHWs as periods when SST exceeds a seasonally varying 90th-percentile threshold for ≥ 5 consecutive days. To this end, daily L4 SST data (GHRSST, https://www.ncei.noaa.gov/archive/accession/GHRSST-MUR-JPL-L4-GLOB, last access: 1 October 2026) were used for the Gulf.

Figure 2 Time series for the Gulf of California of (A) monthly sea surface temperature anomaly (SSTa), (B) annual number of marine heat waves (MHW) days, and (C) multivariate ENSO index version 2 (MEI.v2) from January 2003 to December 2024.

The Multivariate ENSO Index version 2 (MEI.v2) was used to represent ENSO variability during the study period. MEI.v2 is an index for overlapping 2-month periods (December–January, January–February, etc.) based on the leading combined empirical orthogonal function of five ocean–atmosphere variables over the tropical Pacific: sea-level pressure, SST, zonal and meridional surface winds, and outgoing longwave radiation using an empirical orthogonal function (EOF) framework over the tropical Pacific (30° S to 30° N) (Wolter and Timlin, 2011; Kobayashi et al., 2015). The MEI.v2 index is calculated for 12 overlapping 2-month seasons, which allows the seasonal structure of ENSO variability to be represented more appropriately. MEI.v2 values were obtained from National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Laboratory webpage (https://psl.noaa.gov/enso/mei/, last access: 1 October 2026) for the entire study period. ENSO influence was assessed by comparing SST anomaly (SSTa) time series with MEI.v2 (Fig. 2).

In addition, the Pacific Decadal Oscillation (PDO) and the North Pacific Gyre Oscillation (NPGO) indices were included to evaluate the potential influence of large-scale Pacific climate variability on productivity patterns in the Gulf of California. The PDO index represents the dominant mode of monthly SSTa variability in the North Pacific. It is calculated as the leading EOF of monthly SSTa between 20 and 70° N, after removing the monthly climatology and the global mean SSTa. Monthly PDO values were downloaded from the NOAA Physical Sciences Laboratory climate time-series database: https://psl.noaa.gov/data/timeseries/month/PDO/ (last access: 1 October 2026).

The NPGO index represents the second dominant mode of sea surface height variability in the Northeast Pacific. It is derived from an EOF analysis of sea surface height anomalies and reflects changes in the intensity of the central and eastern branches of the North Pacific gyre circulation. The NPGO is also associated with variability in salinity, nutrients, and chlorophyll a in the Northeast Pacific. Monthly NPGO values were downloaded from the North Pacific Gyre Oscillation webpage (https://o3d.org/npgo/, last access: 1 October 2026).

Spearman's rank correlation coefficients were calculated to evaluate the relationships among hydrographic variables, climate indices, and NPP across gulf regions. To account for temporal autocorrelation, statistical significance was assessed using an effective sample size correction following the method of Pyper and Peterman (1998). Monotonic trends were assessed using the Mann–Kendall test to determine their statistical significance and the Theil–Sen slope estimator to quantify their magnitude. Both analyses were performed using the Mann–Kendall and Theil–Sen functions implemented in the Climate Data Toolbox for MATLAB (Greene et al., 2019).

2.3.2 Calculation of stratification index

Water-column stratification was quantified with a density-based index calculated from vertical profiles of potential temperature (θ) and practical salinity (PSU). The stratification index (SI; kg m⁻³) was defined as:

$$\begin{array}{}\text{(1)}& \text{SI}=\mathit{\sigma }\mathit{\theta },\text{MLD}-\mathit{\sigma }\mathit{\theta },\text{surface}\end{array}$$

where: σθ is potential density anomaly referenced to 0 dbar. Practical Salinity was converted to Absolute Salinity (g kg⁻¹) using latitude, longitude, and pressure. Potential density anomaly referenced to 0 dbar was subsequently calculated at the surface and at the mixed-layer depth (MLD), following the Thermodynamic Equation of Seawater 2010 (TEOS-10) (IOC, SCOR, and IAPSO, 2010). With this definition, positive SI values indicate stable stratification (lighter surface waters overlying denser subsurface waters), whereas values near zero denote well-mixed conditions. This metric represents the physical barrier to vertical mixing and links surface warming to reduced nutrient supply (Jacox et al., 2016).

Figure 3 Regionalisation of net primary productivity (NPP) based on (A) the k-means${}^{++}$ clustering of the January 2003–December 2024 global climatology, (B) warm seasonal (June–November), and (C) cold seasonal (December–May). NPP classes: low-productivity (LPZ < 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), mid-productivity (MPZ 900–1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and high-productivity (HPZ > 1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$).

2.3.3 Regional classification of net primary productivity

To analyse spatio-temporal variability in NPP, monthly observations for each grid cell were assigned to one of three productivity classes according to their satellite-derived NPP values. Consequently, the classes represent dynamic productivity conditions rather than fixed geographical regions, and the spatial extent of each class was allowed to vary among months, seasons, and years. Prior to clustering, outliers were removed following the interquartile range (IQR) criterion (1.5 × IQR) of Moore and McCabe (1999). The resulting filtered dataset was used in all subsequent analyses. NPP was then classified using the k-means${}^{++}$ algorithm (Kapoor and Singhal, 2017; Zhang et al., 2025). To evaluate the appropriate number of clusters, we compared alternative solutions using the elbow and silhouette methods (Fig. S1). These diagnostics indicated that a small number of clusters, particularly K = 2–4, provided the most robust classifications, whereas higher values, K≥5, resulted in lower clustering performance and increased spatial fragmentation. Because NPP varies along a continuous environmental gradient, K=3 was selected as the most parsimonious and ecologically interpretable solution, resolving three coherent NPP regimes corresponding to low-, intermediate-, and high-productivity waters.

A discriminant analysis was subsequently performed to evaluate whether the NPP-derived productivity classes corresponded to distinct environmental regimes. The productivity classes obtained from the k-means classification were used as the grouping variable, whereas the environmental variables, SST, EZD, MLD, NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$, SiO₂, URT, UAR, SI, and SSS, were included as explanatory variables (Johnson and Wichern, 2018). NPP was deliberately excluded from this analysis because it was the variable used to define the original k-means${}^{++}$ classification. Its inclusion would therefore have introduced circularity and artificially inflated the classification performance. Thus, the purpose of the discriminant analysis was not to validate the clustering procedure or the selected number of clusters, but rather to determine whether the resulting NPP-based productivity classes could be independently distinguished by hydrographic and biogeochemical characteristics.

Spatial distribution and seasonal patterns (cold season: December–May and warm season: June–November) were visualised (Fig. 3) (Santamaría-del-Ángel et al., 1994; Hidalgo-González and Álvarez-Borrego, 2004). The areal coverage for each cluster was calculated by counting grid cells within each class, converting to km², and then to percentage coverage (%) per month, season, and year. Monthly anomalies were calculated from the coverage time series (Kushnir, 1994; Ramos-Rodríguez et al., 2012). Monthly climatological NPP was estimated for each cluster, and cluster contributions were integrated to obtain total monthly productivity. NPP is reported as teragrams of carbon per month (Tg C month⁻¹) and per year (Tg C yr⁻¹). Annual distribution patterns were visualised (Figs. 5, S2 and S3).

2.3.4 Assessment of NPP drivers

To explain spatial and temporal variability in NPP, the following biophysical variables were analysed for each NPP zone: SST, EZD, MLD, NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$, SiO₂, URT, UAR, SI, and SSS. Monthly time series were computed by zone for each variable. Nutrients (NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$, SiO₂) were vertically integrated from the surface to the EZD prior to analysis.

Potential collinearity among predictors was evaluated with correlation matrices considering Spearman correlations > 0.7; highly correlated variables were excluded from the final predictor set (Table S1 in the Supplement).

To avoid potential circularity between the response variable and predictor variables, SST, euphotic zone depth (EZD), and SiO₂ were excluded from the final GLMM because these variables are directly involved in, or derived from quantities used in, the VGPM algorithm employed to estimate NPP. To identify the environmental drivers of NPP and evaluate whether their effects differed among productivity zones, we fitted a Generalised Linear Mixed Model (GLMM) using a Gamma error distribution and a log-link function, which is appropriate for continuous, positive, and right-skewed response variables such as NPP (Corrales and Cepeda-Cuervo, 2019; Miyamoto et al., 2019). Prior to modelling, all continuous predictor variables were standardised (z-scores) to facilitate comparison of effect sizes and improve model convergence. Separate models were fitted for each productivity zone; within each model, MLD, NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$, URT, UAR, SI, and SSS were included as fixed effects. To account for temporal dependence among observations and repeated measurements through time, month was included as a random intercept term (1|Month). The model structure was:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\text{NPP}}_{\text{zone}}& \sim \text{MLD}+{\mathrm{NO}}_{\mathrm{3}}+{\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}+\text{URT}+\text{UAR}+\text{SI}\\ & +\text{SSS}+\left(\mathrm{1}\mathrm{|}\text{Month}\right)\end{array}\end{array}$$

where NPP_zone indicates the response variable and the same GLMMs structure was fitted independently to each net productivity zone class (low-productivity zone (LPZ), mid-productivity zone (MPZ), and high productivity zone (HPZ) as defined in result section); therefore, zone was used only to stratify the dataset prior to model fitting and was not included as a predictor in the models. MLD = mixed-layer depth, NO₃ = nitrate concentration, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ = phosphate concentration, UAR = along-gulf (u) sea current component; URT = along-gulf (u) wind component; SI = stratification index, SSS = sea surface salinity, and 1|Month represents the random effect, and the environmental variables correspond to fixed effects.

Positive coefficients indicate that higher predictor values are associated with increased NPP, whereas negative coefficients suggest inhibitory or limiting effects. Model fitting was performed using the fitglme function in MATLAB. Model assumptions and goodness of fit were evaluated through graphical and statistical diagnostics. Because a Gamma distribution with a log-link function was employed, diagnostic evaluation focused on the distribution of deviance and Pearson residuals, the absence of systematic patterns in residual-versus-fitted value plots, and the assessment of homoscedasticity across the range of fitted values. Quantile-quantile plots of residuals were also inspected to evaluate potential departures from normality. The adequacy of the linear predictor was examined by inspecting residual patterns against each covariate, while the inclusion of month as a random effect accounted for non-independence among observations collected within the same temporal group. Model performance was further assessed using Akaike's Information Criterion (AIC) and Bayesian Information Criterion (BIC). Fixed-effect coefficients (estimate, standard error, t-statistic, and p-value) were evaluated with marginal ANOVA to identify significant predictors. Predictors with p < 0.05 were considered significant drivers of NPP.

3 Results

3.1 Anomalous warming events

Sea surface temperature anomalies (SSTa) in the Gulf of California were predominantly negative from 2003 to 2013. Beginning in 2014, anomalies transitioned to positive values, with increasing frequency and duration over time. The most intense and persistent positive anomalies occurred between late–2013 and late–2016 (Fig. 2A). During 2014–2016, the annual number of MHW days reached 118, 95, and 51 d, respectively, amounting to 50.2 % of all MHW days in the 2003–2024 record (Fig. 2B). The MEI.v2 also showed strong positive anomalies and extreme values during 2015–2016, consistent with the 2015–2016 El Niño event and was flanked by prolonged negative anomalies associated with La Niña conditions during 2010–2014 and 2016–2018 (Fig. 2C). Based on SSTa, MHW occurrence, and MEI.v2, the warmest period in the study was 2014–2016, followed by 2023–2024, which recorded the second-highest number of MHW days.

3.2 Regional classification of NPP

The k-means${}^{++}$ clustering analysis identified three groups. The discriminant analysis indicated clear environmental separation among the NPP-derived productivity classes. The overall misclassification rate was low (8.7 %), while assignment accuracy was high across groups, with values of 92.4 %, 96.1 %, and 80.1 % for Groups 1–3, respectively. Groups were categorised as low-productivity (LPZ, < 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), mid-productivity (MPZ, 900–1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and high-productivity (HPZ, > 1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) productivity zones (Fig. 3). By dividing the Gulf of California into three productivity zones, we evaluate how anomalous warming affects their temporal spatial extent, their total carbon contribution, and the physical drivers underlying these changes.

3.3 Seasonal and interannual variability in primary productivity

The 2003–2024 NPP mean seasonal cycle ranged from 642 to 3487 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The average annual spatial pattern showed a gulf-latitudinal gradient along the continental coast, extending from the mouth (∼ 20° N) to the northern gulf (∼ 31.5° N) (Fig. 3A). The most productive regions (HPZ, > 1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) included the Nayarit, Sinaloa, and Sonora coasts, the Midriff Archipelago (1200–2000 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and the Colorado River mouth extending ∼ 40 km offshore (Fig. 3A).

Mid-productivity (MPZ, 900–1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) occupied coastal waters near the gulf's mouth and southern sector, forming a coast-shore gradient from the coast to the shelf edge. The MPZ extended from ∼ 25° N to the Midriff Archipelago region and along the entire peninsular coast, with a core in the upper gulf. Low productivity waters (LPZ, < 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) characterised oligotrophic conditions from the mouth to ∼ 26° N (Fig. 3A). Seasonality was pronounced during the warm season (June–November), HPZ and MPZ coverage decreased, while LPZ expanded northward to ∼ 28° N, forming a core in the upper gulf (Fig. 3B). In contrast, during the cold season (December–May), HPZ expanded over the upper gulf and along both continental and peninsular coasts (Fig. 3C).

Figure 4 (A) Monthly mean NPP (bars) for January 2003–December 2024, showing class contribution (%) as stacked bars and total NPP (Tg C month⁻¹), (B) mean annual NPP (Tg C yr⁻¹) per zone for 2003–2024. The dashed line denotes time series mean.

Gulf-wide climatological monthly mean NPP over the 22-year period was 8.6 Tg C month⁻¹. During the cold season, NPP peaked in March at 14.9 Tg C month⁻¹, whereas during the warm season it reached a minimum in August (4.4 Tg C month⁻¹) (Fig. 4A). HPZ coverage increased from 34 % to 71 % (mean 54 %) in cold season but declined to 19 % in warm season, reaching ∼ 8 % in July–August (Figs. 3C and 4A). LPZ covered 8 %–20 % (∼ 10 %) in cold season and expanded to 49 % in warm season, peaking ∼ 66 % in July–August. MPZ accounted for ∼ 34 % annually (Figs. 3B and 4A).

Gulf-wide annual NPP per productivity zone ranged from 84 to 198 Tg C yr⁻¹ (mean of 124 ± 27.3 Tg C yr⁻¹) (Fig. 4B). Contributions over the entire period were HPZ = 61.5 %, MPZ = 21.6 %, and LPZ = 16.9 % (Fig. 4B). During 2003–2012 (predominantly negative SSTa), HPZ contributed 63.9 %, MPZ = 20.9 %, and LPZ = 15.2 %. During 2013–2024 (positive SSTa), LPZ increased to 18.4 %, while HPZ declined to 59.4 %, indicating a progressive reduction in HPZ extent, MPZ remained ∼ 21.6 %. Years 2006–2012 exceeded the historical mean (124 ± 27.3 Tg C yr⁻¹), with 2008, 2011, and 2012 being most productive (167, 198, and 164 Tg C yr⁻¹, respectively). In those years, HPZ contributed ∼ 78 % of NPP, while LPZ contributed 14 %–17 %. In contrast, the anomalously warm interval (2014–2016) recorded the lowest mean NPP (∼ 101 Tg C yr⁻¹), with LPZ contributing 21.8 % and HPZ declining to 54.7 % in 2016. MPZ exceeded 25 % of NPP in 2004, 2013, 2018, and 2021 (Fig. 4B). The longest below-average period occurred during 2013–2019, followed by recovery in 2020–2024, only 2023 (111 Tg C yr⁻¹) remained below the mean. Notably, 2024 ranked fourth (152 Tg C yr⁻¹), after 2008, 2011, and 2012.

Figure 5 (A) NPP annual maps showing low-productivity (LPZ), mid productivity (MPZ), and high productivity (HPZ) variability during 2003–2024; (B) anomalies in class coverage (%) through time for (A) low (LPZ), (B) mid (MPZ), and (C) high (HPZ). Cold (blue) represents a reduction in cover, and warm (red) represents an increase in cover relative to the historical average distribution.

3.4 Spatial and interannual variability of primary productivity

The annual spatial distribution of the LPZ, MPZ and HPZ (Fig. 5) exhibited distinct patterns closely tied to temperature regimes. Here, percentages refer to the spatial extent of each zone. During 2003–2012, when SSTa were predominantly negative, the LPZ was mostly located south of 26° N, covering on average 32.4 % of the gulf. In these years, HPZ was prominent in the upper gulf, the Midriff Archipelago, and along the continental coast, extending > 40 km offshore, and accounting for 33.3 % of the gulf area. In contrast, during 2013–2023, when SSTa were mostly positive, HPZ coverage decreased markedly in coastal areas and the upper gulf (< 26 % of the area), while MPZ expanded to ∼ 33.9 %, occupying much of the basin. The most pronounced spatial reorganisation occurred during the 2014–2016 warm anomaly period, when LPZ increased to nearly 50 %, advanced north of 26° N and reached ∼ 28° N in 2016. Years with widespread HPZ coverage (e.g., 2008, 2011, 2012, 2021 and 2024) coincided with positive NPP anomalies > 25 %, whereas 2013–2016 showed LPZ expansions > 20 % (Fig. 5).

Seasonality was most pronounced in the cold season (December–May). On average, HPZ covered 44.3 % of the gulf, followed by MPZ (33.5 %) and LPZ (22.2 %). HPZ varied from 22.9 % (2016) to 63 % (2012), MPZ from 25.7 % (2010) to 42.5 % (2013), and LPZ from 2 % (2021) to 44 % (2016). HPZ exceeded 50 % in 2008 (58 %), 2011 (61.4 %), 2012 (63 %), and 2021 (59.7 %), with LPZ nearly absent (< 10 %) in 2021. Conversely, HPZ fell below the historical average in 2014, reaching 22.9 % in 2016 and 42 % in 2017 (Fig. S2).

During the warm season (June-November), interannual variability was lower. LPZ dominated averaging 68.5 % (range 57.9 %–76.2 %) and shifting northward into the upper gulf. MPZ averaged 19.7 % (15.6 %–26.9 %), and HPZ averaged 11.7 % (7 %–22 %) concentrated near the Midriff Archipelago, Colorado River mouth, and along the continental coast, with limited offshore extent (Fig. S3). Additional latitude-time distributions of NPP and SST anomalies are shown in Figs. S4 and S5.

Table 1 GLMM results for main effects of biophysical variables on NPP by productivity zone classification (LPZ, MPZ, and HPZ). Coefficients p < 0.05 are shown in bold. Abbreviations: LPZ = low-productivity zone; MPZ = mid-productivity zone; HPZ = high-productivity zone; MLD = mixed layer depth; NO₃ = nitrate concentration; ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ = phosphate concentration; UAR = along-gulf (u) sea current component; and URT = along-gulf (u) wind component; SI = stratification index, and SSS = sea surface salinity.

3.5 Responses of biophysical drivers (GLMM analysis)

Correlation analysis indicated strong collinearity between ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ and SiO₂ (r = 0.88; p < 0.05). ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ was retained in the model given its broader relevance to primary productivity, whereas SiO₂ primarily constrains siliceous phytoplankton, like diatoms. Best-fitting Generalised Linear Mixed Models (GLMMs) for each NPP class identified the significant environmental drivers (Table 1). The model provided a good fit to the data (AIC = 1057.4, BIC = 1178.9, log-likelihood = −502.69). Month was included as a random intercept to account for temporal dependence among observations, yielding an estimated random-effect standard deviation of 0.504.

The stratification index (SI), sea surface salinity (SSS), and nitrate (NO₃) were significant predictors of NPP across the three productivity zones, although their effects differed in both magnitude and direction (Table 1). In the LPZ, SI was positively associated with NPP (β = 1.050, p < 0.001), whereas salinity (β = −0.176, p < 0.001) and nitrate (β = −1.841, p < 0.001) showed significant negative relationships. In contrast, both the MPZ and HPZ exhibited the opposite pattern. SI was negatively associated with NPP (MPZ: β = −1.116; HPZ: β = −1.117; p < 0.001), whereas salinity (MPZ: β = 0.196; HPZ: β = 0.182; p < 0.001) and nitrate (MPZ: β = 1.926; HPZ: β = 2.008; p < 0.001) were positively associated with productivity. The effect of URT also varied among productivity zones. URT was positively associated with NPP in the LPZ (β = 0.326, p < 0.001), whereas a significant negative relationship was observed in the HPZ (β = −0.336, p < 0.001). In contrast, URT was not a significant predictor in the MPZ. Phosphate (${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$) was significant only in the HPZ, where it exhibited a positive association with NPP (β = 0.294, p < 0.001). These results underscore strong regional differences in environmental controls on NPP across the Gulf of California. The complete anomaly time series for all variables (2003–2024) is provided in the Supplement (Figs. S4–S9, and Table S1).

4 Discussion

Monitoring NPP in the Gulf of California is essential for assessing the availability of phytoplankton production to upper trophic levels and for evaluating the ecological and socioeconomic impacts of climate variability and anthropogenic pressures, and their ecological and socioeconomic consequences. NPP underpins fisheries, mediates carbon cycling, and contributes to ecosystem resilience, making it a key indicator of ocean health. Despite its ecological importance, long-term NPP records, both in situ and satellite-derived, remain scarce, constraining assessments of variability from seasonal to interannual scales.

The Gulf of California's spatial organisation has been analysed through various regionalization schemes, evolving from early classifications using biological indicators and physical properties to more complex schemes leveraging satellite observations. Initial frameworks identified regions based on phytoplankton community structure (Round, 1967; Zeitzschel, 1969), thermohaline structures (Roden and Emilsson, 1979) and bathymetry (Merrifield and Winant, 1989), while later methods introduced detailed classifications based on photosynthetic pigments (Santamaría-del-Ángel et al., 1994) and multivariate variables and machine-learning techniques (Santamaría-del-Ángel et al., 1994; Kahru et al., 2004; Petatán-Ramírez, 2015; Heras-Sánchez et al., 2019; Marín-Enríquez et al., 2024). Recent approaches exhibit increased complexity, ranging from simplified divisions to schemes with up to 22 regions. Nonetheless, effective regionalisation should prioritize parsimony and mechanistic understanding to maintain interpretability and enhance insights into productivity controls.

We utilised the K-means${}^{++}$ algorithm with k-means${}^{++}$ initialisation for its probabilistic centroid selection, enhancing accuracy and convergence in classifying meteorological, hydrographic (Kapoor and Singhal, 2017; Zhang et al., 2025; Zhang and Qian, 2025) and primary productivity patterns in the Gulf of California. We established three NPP categories: HPZ (> 1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), MPZ (900–1700 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$), and LPZ (< 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Although the HPZ covers∼ 30 % and contribute ∼ 61.4 % of total NPP, LPZ covers ∼ 36 % and accounts for only ∼ 16.9 %. These findings align with previous studies indicating a latitudinal gradient in productivity linked to factors such as upwelling and tidal currents, especially in areas with strong mixing along the eastern coast and Midriff Archipelago (Valdez-Holguín and Lara-Lara, 1987; Gaxiola-Castro et al., 1995).

HPZ regions are found along the eastern coast, the Midriff Archipelago, and the upper gulf, characterised by strong mixing and upwelling (Mercado-Santana et al., 2017), particularly during winter when north-westerly winds influence the coasts of Sonora, Sinaloa, and Nayarit (Santamaría-del-Ángel et al., 1994; Mercado-Santana et al., 2017). MPZ contributes ∼ 21 % of total NPP with moderate conditions (SST ∼ 25 °C, salinity ∼ 35 ups, and a mixed layer depth ∼ 24 m), influenced by various water types that enhance productivity through complex interactions with mesoscale processes (Lavín and Marinone, 2003; Lavín et al., 2013). In contrast, LPZ is the most oligotrophic area of the southern gulf, with historically low NPP, largely due to eastern Pacific influences and warm, low-salinity waters that increase stratification and limit nutrient availability, further complicated by thermohaline structures and mesoscale dynamics (Mercado-Santana et al., 2017; Marín-Enríquez et al., 2024; Trasviña-Castro et al., 2024).

The GLMM analysis revealed that the environmental controls on NPP vary markedly among productivity regimes, highlighting the spatial heterogeneity of the mechanisms regulating phytoplankton NPP across the Gulf of California. Despite these regional differences, NO₃ availability and water-column stratification (SI) consistently emerged as the dominant predictors of NPP. The HPZ and MPZ exhibited the strongest dependence on nutrient availability, with NO₃ representing the most influential positive predictor in both spatio-temporal variable zones, whereas ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ contributed significantly only within the HPZ. In contrast, increasing stratification (SI) exerted a strong negative effect in both NPP zones, indicating that enhanced water-column stability limits the upward transport of nutrients into the euphotic layer. These results are consistent with the hydrographic dynamics of the upper gulf, the Midriff islands region, and eastern continental coast, where seasonal cooling, wind-driven mixing, and tidal processes enhance the upward transport of nutrients and sustain high rates of primary production (Santamaría-del-Ángel et al., 1994; Mercado-Santana et al., 2017; Marín-Enríquez et al., 2024). The MPZ appears to represent a transitional environment where biological production is primarily regulated by nutrient supply, and where stronger stratification restricts the vertical nutrient transport into the euphotic zone (Marín-Enríquez et al., 2024). Although nitrogen is widely recognised as the principal limiting nutrient in the Gulf of California (Álvarez-Borrego, 2012), our results indicate that ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ availability also significantly affects NPP within the HPZ. This pattern is consistent with observations from the Guaymas Basin, where intense summer stratification restricts nutrient fluxes, promoting nitrogen limitation and favouring picophytoplankton-dominated communities (Millán-Núñez et al., 2023; Martínez-López et al., 2023). In contrast, the LPZ, influenced by oligotrophic, highly stratified waters entering from the Pacific Ocean, remains persistently low in productivity and is typically dominated by picophytoplankton. Consequently, the LPZ maintains relatively low but persistent productivity, whereas the MPZ exhibits greater seasonal variability associated with winter mixing and episodic nutrient enrichment.

At the seasonal scale, NPP peaks in the cold season (December–May), with maxima in February-March (Fig. S2). This increase is driven by prevailing north-westerly winds and cooling that enhance upwelling along the eastern gulf coast, enriching the euphotic zone and stimulating phytoplankton growth (Santamaría-del-Ángel et al., 1994; Escalante et al., 2013). During this period, HPZ and MPZ together typically cover ∼ 77 % of the gulf. Hakspiel-Segura et al. (2022) reported a negative correlation of SST and meridional wind stress with surface Chl a concentrations in the central gulf, reinforcing the role of wind-driven mixing.

Conversely, during the warm season (June–November), oligotrophic conditions dominate: NPP < 900 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ covers ∼ 68 % of the basin due to weaker southeasterly winds, the intrusion of warm tropical waters, and strong stratification that limits nutrient supply (Figs. S3, S4, S8, and S9) (Lavín and Marinone, 2003; Escalante et al., 2013; Mercado-Santana et al., 2017; Trasviña-Castro et al., 2024).

The 2003–2024 NPP time series exhibited significant interannual variability. Spearman's rank correlation analyses indicated that neither the PDO nor the NPGO was significantly associated with the spatial coverage of any productivity zone (Table S1). In contrast, MEI.v2 showed a weak but statistically significant negative correlation with the extent of the HPZ (rs_[264] = −0.136, p_Corr = 0.027), suggesting that ENSO-related variability exerts a measurable, although limited, influence on regional productivity patterns. Consistent with this relationship, La Niña conditions were generally associated with lower SST, intensified north-westerly winds, a deeper mixed layer, and enhanced nutrient availability (NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ and SiO₂) within the euphotic zone, conditions that favour elevated NPP. During the La Niña events of 2008, 2011, 2012, and 2021, the HPZ contributed more than 70 % of the total gulf-wide NPP while occupying only 35 %–44 % of the gulf area. This response was most pronounced during the cold season, when the HPZ and MPZ together occupied approximately 95 % of the Gulf of California, confining the LPZ primarily to the southern entrance of the basin (Fig. S2).

Stratification index (SI) analysis indicated varying productivity regimes (Fig. S9): HPZ had low SI (0–0.5 kg m⁻³) due to mixing processes (upwelling, tidal forcing, and north-westerly winds), while LPZ exhibited high SI (up to ∼ 2.0 kg m⁻³), indicating strongly stratified waters. MPZ showed an intermediate SI (0.5–1.1 kg m⁻³). Notably, from mid–2013, warming trends led to increased stratification (Figs. S5 and S9), reduced nutrient concentrations (NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$), and negative anomalies in NPP (Fig. S8). This pattern mirrors global trends of warming-driven stratification documented in the North Atlantic, North Pacific, and elsewhere (Xie et al., 2019; Fischer et al., 2020; Honda, 2020), with cascading effects on hydrography, primary productivity, phytoplankton functional composition, and coastal carbon fluxes (Xie et al., 2019; Fischer et al., 2020; Honda, 2020; Hakspiel-Segura et al., 2022; Song et al., 2022)

The 2014–2016 period marked a significant decline in NPP in the Gulf of California, characterised by severe negative anomalies. This period was associated with strong atmospheric forcing related to ENSO (MEI.v2) and coincided with a positive phase of the PDO, a transition toward a negative NPGO, and widespread marine heatwaves. Together, these large-scale climate conditions have been previously linked to changes in the physical and biogeochemical structure of the Gulf of California (Hakspiel-Segura et al., 2022). The persistent increase in SST during this period diminished cold season HPZ coverage to below 5 % and expanded the LPZ to approximately 50 % of the gulf. SST anomalies exceeding +2 °C led to intensified thermal stratification and limited nutrient transport into the euphotic zone, conditions that are consistent with reduced productivity in the HPZ and MPZ, which depend on vertical mixing and coastal upwelling for nutrient supply (Di Lorenzo and Mantua, 2016; Hakspiel-Segura et al., 2022). Stratification increased significantly, reaching oligotrophic conditions, and caused diminished nutrient flux, correlating with declines in PO₄ and a reduction in productive habitats.

A key mechanism affecting the marine ecosystem between 2014 and 2016 was the weakening of north-westerly winds (∼ 8 m s⁻¹ versus climatological mean > 20 m s⁻¹), which reduced coastal upwelling efficiency. This led to a shift in seasonal productivity maximum and a decline in HPZ coverage. Concurrently, Tropical Surface Water intrusion reinforced stratification and limited mesoscale activity, particularly in the southern gulf and at the mouth (Di Lorenzo et al., 2008, 2013; Di Lorenzo and Mantua, 2016), resulting in the lowest NPP recorded, with both HPZ and MPZ declining by about 7 %, and areal coverages below historical averages. The HPZ remained only in the upper gulf, while the MPZ shifted north, and the LPZ expanded to ∼ 26–27° N (Fig. S3). Additionally, nutrient records indicated a significant decline in phosphate concentrations, suggesting a N-P co-limitation and highlighting phosphorus's critical role during extreme warm events (Browning and Moore, 2023).

Long-term trends derived from this study indicate declining NPP of 50–100 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{year}}^{-\mathrm{1}}$ in the upper gulf, in La Paz Bay, and off the coast of Nayarit (Fig. S10), while values at the mouth and near ∼ 28° N remained relatively stable. These patterns are consistent with global trends reported by Maishal (2024a, b), who documented ∼ 25 $\mathrm{mg}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{decade}}^{-\mathrm{1}}$ decreases in the Pacific associated with warming of 0.1–0.2 °C decade⁻¹, rates slightly lower than those reported for the Gulf of California (0.38 and 0.57 °C decade⁻¹; Hakspiel-Segura et al., 2022; Sánchez-Cabeza et al., 2022). Maishal (2024a) also noted concurrent decreases in density, salinity, wind stress, MLD, and nutrients (NO₃, ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$) during the past two decades, in agreement with our findings. Despite this, NPP rebounded strongly after the prolonged low-productivity phase of 2013–2019: 2020–2024 showed recovery, and although 2023, the warmest oceanic year on record (Cheng et al., 2024), did not yield the lowest NPP (that occurred in 2016, 85 Tg C yr⁻¹), 2024 ranked fourth (152 Tg C yr⁻¹), suggesting resilience potentially linked to the negative PDO that began in late 2019.

Globally, NPP and Chl a with ocean warming are well documented, along with strong coupling between NPP and climate variability (Behrenfeld et al., 2006; Xie et al., 2019; Fischer et al., 2020; Honda, 2020). Positive phases of climate indices such as MEI.v2 and PDO exert negative feedback on NPP via warm-water intrusion, enhanced stratification, and oligotrophication that alter pressure centres and wind regimes, reducing upwelling and nutrient supply.

This study advances understanding of regional and phenological mechanisms affecting NPP in the Gulf of California, primarily mediated by SST and its associated processes. The integration of the long-term environmental dataset demonstrated the utility of satellite-based and modelled observations for tracking NPP variability. Our parsimonious three-zone classification (LPZ, MPZ, and HPZ) captures spatial and temporal heterogeneity efficiently. Further research should elucidate the drivers of the LPZ expansion observed over the past decade and its implications for ecosystem structure and fisheries in the Gulf of California.

5 Conclusions

The average total NPP in the Gulf of California during 2003–2024 was 124 Tg C yr⁻¹. Classification of the Gulf of California into three NPP categories with the k-means${}^{++}$ algorithm provided an effective, parsimonious framework for explaining spatial and temporal variability. HPZ conditions predominated along the continental coast and in the Midriff Archipelago, contributing 61.4 % of total NPP while covering approximately 30 % of the Gulf area. In contrast, LPZ conditions dominated the oceanic region and the entrance to the Gulf, covering approximately 36 % of the area and contributing 16.9 % of total NPP.

The decline in NPP in the Gulf of California during 2014–2016 resulted from the non-linear interaction El Niño-related surface warming, weakened wind forcing, tropical-water intrusion, enhanced thermal stratification, reduced nutrient availability, and concurrent marine heatwaves. This coupling explains both the magnitude of the productivity decline and the unprecedented expansion of oligotrophic conditions observed during these years. Seasonal variability was also pronounced: the cold season, from December to May, accounted for 68 % of annual NPP, whereas the warm season, from June to November, contributed 32 %.

The GLMM results indicate that NO₃and water-column stratification (SI) were the dominant predictors of NPP overall. In the HPZ, NPP was additionally modulated by the along-gulf wind, the salinity, and the ${\mathrm{PO}}_{\mathrm{4}}^{-\mathrm{3}}$ concentrations. Warmer years were associated with reduced coverage of productive zones, including mesotrophic (MPZ) and eutrophic (HPZ), and with the expansion of oligotrophic LPZ areas. If regional warming and stratification trends persist, oligotrophic conditions and low productivity are likely to intensify, potentially extending the influence of LPZ further into the Gulf of California. These changes may have important ecological and socioeconomic implications, underscoring the need for continued monitoring and research to anticipate, understand, and mitigate the impacts of climate variability on marine phytoplankton net primary productivity.